% articulo Eskinat AIChE Journal february 1991, vol. 37, No. 2
% elaboro eqm, 23oct'17
% adaptado de programa 5.7, Luyben, pag. 130. (eqm) 18 nov 98.
% "binary distillation column dynamics" 
%-------------------------------------------------------------
% 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=25; nf=12; mdo=0.5; mbo=0.5; mo=0.5; ro=1.477;
vo=1.977; f=1; beta=0.1; alpha=2; r=ro; v=vo;
xb=0.005; xd=0.995; Ts=10; nTs=0;
x=[8.9448e-03, 1.5161e-02, 2.4848e-02, 3.9711e-02, 6.1978e-02,...
 9.4169e-02, 1.3839e-01, 1.9507e-01, 2.6156e-01, 3.3197e-01,...
 3.9885e-01, 4.5616e-01, 5.0178e-01, 5.5763e-01, 6.2155e-01,...
 6.8930e-01, 7.5550e-01, 8.1527e-01, 8.6548e-01, 9.0517e-01,...
 9.3506e-01, 9.5677e-01, 9.7212e-01, 9.8277e-01, 9.9006e-01];
kcd=1000; kcb=1000; taud=5; taub=1.25; delta=0.005;
time=0; tprint=0; erintd=0; erintb=0; xx=[];
a1=-0.796; b1=0.204; g1=1.01; g2=-7.13; g3=-7.43;
% disturbance
z=0.5;
% 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<=9000)
   % tray liquid hydraulics and VLE
%    if(time>0 && time<100), r=ro-0.05; end
%    if(time>100 && time<200), r=ro-0.1; end
%    if(time>200 && time<300), r=ro; end
%    if(time>300 && time<400), r=ro-0.05; end
%    if(time>400 && time<500), r=ro+0.1; end
%    if(time>500 && time<600), r=ro; end
%    if(time>600 && time<700), r=ro+0.05; end
%    if(time>700 && time<800), r=ro-0.15; end
%    if(time>800 && time<900), r=ro+0.1; end
%    if(time>900 && time<1000), r=ro+0.15; end
   for n=1:nt
      y(n)=alpha*x(n)/(1+(alpha-1)*x(n));
      l(n)=r+f;
      if(n>nf) l(n)=r; end
   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)*0;
   %r=ro+kcd*(errd+erintd/taud)*0;
   % 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;
   [ r,v ] = RV_time( time,ro,vo );
   if(time>=tprint)
      %[ r,v ] = RV_time( time,ro,vo );
      xx=[xx;nTs,xb,x(10),xd,r,v];
      tprint=tprint+Ts; nTs=nTs+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
u=xx(:,5)-ro; yt(1)=0; dR=u; dXd=xx(:,4)-0.995;
for i=1:nTs
    xt(i)=g1*u(i)+g2*u(i)^2+g3*u(i)^3;
    yt(i+1)=b1*xt(i)-a1*yt(i);
end
dV=xx(:,6)-vo;
[xx1,yy1]=stairs(xx(:,1),dR);
[xx2,yy2]=stairs(xx(:,1),dV);
figure
subplot(2,1,1),plot(xx1,yy1),ylabel('Delta(R)'),grid
subplot(2,1,2),plot(xx2,yy2),ylabel('Delta(V)'),grid
%subplot(2,1,2),plot(xx(:,1),xx(:,4)-0.995,u,yt(1:length(u)),'--')
figure
plot(xx(:,1),dXd);%,xx(:,1),yt(1:length(xx(:,1))),'--')
ylabel('Delta(Xd)'),grid

