# 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

### Accepted Answer

Sara
on 9 Jul 2014

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

### Categories

### Community Treasure Hunt

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

Start Hunting!