Help with this code (fixed pivot technique)

Question

Cesar García 2015-11-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/252223-help-with-this-code-fixed-pivot-technique

评论： Walter Roberson 2015-11-2

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2015-11-2

"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?

请先登录，再进行评论。

请先登录，再回答此问题。