% programa 5.7, W.L.Luyben, pag. 130
% "binary distillation column dynamics" 
% modif. E.Quintero.M 18nov'98
%------------------------------------------------------
% assumptions: constant relative volatility, equimolal
%     overflow, theoretical trays, simple liquid tray 
%     hydraulics
% feedback controllers manipulate R and V to control XD
% and XB, disturbance is a feed composition change from
% 0.50 to 0.55 at time equal zero
clear
format short e
% initial conditions and parameter values
nt=20; nf=10; mdo=100; mbo=100; mo=10; ro=128.01;
vo=178.01; f=100; beta=0.1; alpha=2;
xb=0.02; xd=0.98;
x=[.035,.05719,.08885,.1318,.18622,.24951,.31618,...
 .37948,.43391,.47688,.51526,.56295,.61896,.68052,...
 .74345,.80319,.85603,.89995,.93458,.96079];
kcd=1000; kcb=1000; taud=5; taub=1.25; delta=0.005;
time=0; tprint=0; erintd=0; erintb=0; xx=[];
% disturbance
z=0.55;
% condiciones iniciales
for n=1:nt
   m(n)=mo;
   mx(n)=m(n)*x(n);
   lo(n)=ro+f;
   if(n>nf) lo(n)=ro; end
end
while(time<=50)
   % tray liquid hydraulics and VLE
   for n=1:nt
      y(n)=alpha*x(n)/(1+(alpha-1)*x(n));
      l(n)=lo(n)+(m(n)-mo)/beta;
   end
   yb=alpha*xb/(1+(alpha-1)*xb);
   % two PI feedback controllers
   errb=0.02-xb;
   errd=0.98-xd;
   v=vo-kcb*(errb+erintb/taub);
   r=ro+kcd*(errd+erintd/taud);
   % perfect level controllers in reflux drum and column base
   d=v-r;
   b=l(1)-v;
   if(r<0 | v<0 | d<0 | b<0)
      error('level too low or composition unreal');
   end
   % evaluate derivatives
   xbdot=(l(1)*x(1)-v*yb-b*xb)/mbo;
   mdot(1)=l(2)-l(1);
   mxdot(1)=v*(yb-y(1))+l(2)*x(2)-l(1)*x(1);
   for n=2:nf-1
      mdot(n)=l(n+1)-l(n);
      mxdot(n)=v*(y(n-1)-y(n))+l(n+1)*x(n+1)-l(n)*x(n);
   end
   % feed plate
   mdot(nf)=l(nf+1)-l(nf)+f;
   mxdot(nf)=v*(y(nf-1)-y(nf))+l(nf+1)*x(nf+1)-l(nf)*x(nf)+f*z;
   for n=nf+1:nt-1
      mdot(n)=l(n+1)-l(n);
      mxdot(n)=v*(y(n-1)-y(n))+l(n+1)*x(n+1)-l(n)*x(n);
   end
   mdot(nt)=r-l(nt);
   mxdot(nt)=v*(y(nt-1)-y(nt))+r*xd-l(nt)*x(nt);
   xddot=v*(y(nt)-xd)/mdo;
   if(time>=tprint)
      xx=[xx;time,xb,x(10),xd,r,v];
      tprint=tprint+0.1;
   end   
   % euler integration
   time=time+delta;
   xb=xb+delta*xbdot;
   for n=1:nt
      m(n)=m(n)+mdot(n)*delta;
      mx(n)=mx(n)+mxdot(n)*delta;
      x(n)=mx(n)/m(n);
      if(x(n)<0 | x(n)>1 )
         error('level too low or composition unreal');
      end
   end
   xd=xd+xddot*delta;
   erintd=erintd+errd*delta;
   erintb=erintb+errb*delta;
end
subplot(2,1,1),plot(xx(:,1),xx(:,6)),ylabel('XB'),grid
subplot(2,1,2),plot(xx(:,1),xx(:,5)),ylabel('XD'),grid
