Different outcome by using self time iteration and ode15s.
2 views (last 30 days)
Show older comments
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
on 9 Jul 2014
Can you attach a code that works? There are some missing inputs in the code you posted, x for instance.
Accepted Answer
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
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.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!