Help with this code (fixed pivot technique)

2 visualizaciones (últimos 30 días)
Cesar García
Cesar García el 2 de Nov. de 2015
Comentada: Walter Roberson el 2 de Nov. de 2015
Hi, Regards, i need help with this code, in order to get the molecular weight distribution of the polymer. This is using the fixed pivot technique.
here, i define the optimal stimated parameters:
clear all
clc
Param=[4.3276103255301335 0.10215543561231799 12.346130949956718 0.04565511507929096 0.0074744247798137625 0.0022142903501092513 1.6957452439712215 0.008166212885461337 0.028435528710812962 6.156650857937292 2.0566837008694536 0.9088238699999999 0.8669780915448081 8.072715437027492 0.010453310496367572 5.260491048 5.325728554979126 6.837140753063702 0.03947867100809781];
Optim=[1.4989542634814161 1.4335418277418217 1.5744495703184183 0.2559519340674219 0.7975221812123433 1.3319125502728477 1.1153538963112686];
Sens=[0.5 0.973342635 2 8.827704916 0.87519757 1.2];
% Parameters for Fed-Batch Models
um=0.54*Param(1)*0.6*0.6*0.5*Optim(1);
Ksr=0.15*Param(2);
Sm=0.3*Param(3)*Optim(2);
Csx=1.7640*Param(4)*0.59*25*Optim(3);
Rcsx=0.022*Param(5)*Optim(4);
Csp=1.7460*Param(6);
k1=0.14*Param(7)*1.5*Optim(5);
k2=0.74*Param(8);
Cnx=0.336*Param(9);
KL=0.1*Param(10);
k3=93*Param(11);
k4=85*Param(12);
O2Leq=7.6*10^-3*Param(13);
nk=1.22*Param(14)*Optim(6);
Rcnx=0.16*Param(15)*Optim(7);
alfa1=0.143*Param(16);
alfa2=0.0000004*Param(17);
alfa3=0.0001*Param(18);
Kox=0.02*Param(19);
Now i define the state variables and feeding strategy for fed-batch fermentation solved by the Euler´s method
%Parameters for polymerization model
ki=0.62*10^4;
kp=0.46*10^5;
kt=0.14*10^1;
km1=0.11*10^-3;
km2=0.86*10^7;
kd=0.83*10^2;
%Initial Conditions validation data
X(1)=1.16;
S(1)=20.05;
N(1)=2;
P(1)=1.01;
O2L(1)=0.002432;
CO2(1)=0.01;
V(1)=4;
V1(1)=0;
V2(1)=0;
%Initial Conditions for polymerization
M(1)=0;
ESHM(1)=0;
dp1dt(1)=0;
dp2dt(1)=0;
dDdt(1)=0;
F3=0.1;
Sin=300*Sens(3);
Nin=7*Sens(4);
O2in=0.002432;
Yms=3.48*10^-3;
Jf=0.2818;
Jm=Yms*Jf;
%Specific growth rate
u(1)=um*((N(1)/S(1))/((N(1)/S(1))+Ksr))*(1-((N(1)/S(1))/Sm)^nk)*((O2L(1)/(Kox*X(1)+O2L(1))));
tsim=50;
t(1)=0;
dt=0.001;
i=1;
Nt=2;
%Euler Method for ODE´s solution
while t(i)<tsim
if t(i)<6
F1=0;
F2=0;
elseif t(i)>6 && t(i)<10
F1=0;
F2=70*(1/1000)*Sens(5);
Fobj31=1.3*F2*4;
elseif t(i)>10 && t(i)<16
F1=80*(1/1000)*Sens(6);
F2=70*(1/1000)*Sens(5);
Fobj32=3.2*F1*6+1.3*F2*6;
elseif t(i)>16 && t(i)<20
F1=80*(1/1000)*Sens(6);
F2=0;
Fobj33=3.2*F1*4;
Fobj3=Fobj31+Fobj32+Fobj33;
else
F1=0;
F2=0;
end
V(i+1)=V(i)+(F1+F2)*dt;
V1(i+1)=F1;
V2(i+1)=F2;
% Fed-Batch States
X(i+1)=X(i)+(u(i)*X(i)*dt-((F1+F2)/V(i))*X(i)*dt);%Biomass
S(i+1)=S(i)-((Csx*u(i)*X(i))+(Rcsx*X(i))+Csp*((k1*u(i)*X(i))+(k2*X(i))))*dt+(F1/V(i))*Sin*dt-((F1+F2)/V(i))*S(i)*dt;%Carbon Source
P(i+1)=P(i)+((k1*u(i)*X(i))+(k2*X(i)))*dt-((F1+F2)/V(i))*P(i)*dt;%Polymer
N(i+1)=N(i)-((Cnx*u(i)*X(i))+(Rcnx*X(i)))*dt+(F2/V(i))*Nin*dt-((F1+F2)/V(i))*N(i)*dt;%Nitrogen Source
if N(i+1)<0.15
N(i+1)=0.15;
end
O2L(i+1)=O2L(i)+((KL*(O2Leq-O2L(i)))-(k3*u(i)*X(i))-((k4*k1*u(i)*X(i))+(k4*k2*X(i))))*dt+(F3/V(i))*O2in*dt-((F1+F2)/V(i))*O2L(i)*dt;%Dissolved Oxygen
if O2L (i+1)<0.002432
O2L(i+1)=0.002432;
end
CO2(i+1)=CO2(i)+((alfa1*u(i)+alfa2)*X(i)*dt)+(alfa3*dt)-((F1+F2)/V(i))*CO2(i)*dt;%Dissolved CO2
u(i+1)=um*((N(i+1)/S(i+1))/((N(i+1)/S(i+1))+Ksr))*(1-((N(i+1)/S(i+1))/Sm)^nk)*((O2L(i)/(Kox*X(i)+O2L(i))));
if u(i+1)<0
u(i+1)=0;
end
Now, i am getting problem with the next part
%Monomer and Complex definition
dp1dt(1)=1;%
for n=1:Nt
M(i+1)=M(i)+ (-((F1+F2)/V(i))*M(i)*dt+(-((Csx*u(i)*X(i))+(Rcsx*X(i))+Csp*((k1*u(i)*X(i))+(k2*X(i))))*Yms)-(km1*M(i)*dt)-(km2*M(i)*dp1dt(n))*dt);
ESHM(i+1)=ESHM(i)+(-((F1+F2)/V(i))*ESHM(i)*dt+(km1*M(i)*dt)-(ki*(ESHM(i)*dt)));
for j=1:Nt
for k=1:j
dp1dt(j)=(-((F1+F2)/V(j))*dp1dt(j)+(ki(ESHM))-(km2*(dp1dt(j)*M))+kp*(dp2dt(k)*A(j,k)-kt*dp1dt(j)));
dp2dt(j)=(-((F1+F2)/V(j))*dp2dt(j)+(km2*(dp1dt(j)*M))-(kp(dp2dt(j))));
dDdt(j)=(-((F1+F2)/V(j))*dDdt(j)+(kt*dp1dt(j))-(kd*dDdt(j))+(kd*(dDdt(k)*B(j,k))));
end
end
end
t(i+1)=t(i)+dt;
i=i+1;
end
besides, A(j,k) and B(j,k) need to be solved by the fixed pivot technique, i´ll appreciate all your help.
  1 comentario
Walter Roberson
Walter Roberson el 2 de Nov. de 2015
"i am getting problem with the next part" is not very specific. Is the code meowing and demanding to be fed Purina Code Chow? Do vandals keep coming by and defacing the code? Is the code refusing to pay its taxes?

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Chemistry en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by