Clear Filters
Clear Filters

Different outcome by using self time iteration and ode15s.

2 views (last 30 days)
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 Comments
Sara
Sara on 9 Jul 2014
Can you attach a code that works? There are some missing inputs in the code you posted, x for instance.
Bendo Voon
Bendo Voon on 9 Jul 2014
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!

Sign in to comment.

Accepted Answer

Sara
Sara on 9 Jul 2014
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 Comments
Sara
Sara on 10 Jul 2014
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
Bendo Voon on 11 Jul 2014
Edited: Bendo Voon on 11 Jul 2014
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!