$*************************************************************
$ Dynamic Optimization of a Binary Distillation Column
$ C. A. Schweiger and C. A. Floudas
$
$ ~Step disturbance in feed composition
$
$ Optimal Solution:
$*************************************************************

OPTION {{
  INTCUT;
  FTOL = "1e-6";
  OTOL = "1e-6";
}}

DECLARATION {{
  INDEX {i};
  SET  T  = |1:30|; #Trays
  SET  T1 = |0:31|; #Reboiler, Trays, and Condenser

  PARA alpha = 2.5;
  PARA N = 30; #possible number of trays
  PARA m = 0.175;
  PARA xdset = {0.98};
  PARA xbset = {0.02};

  XVAR {f(T), #feed flowrate to each tray
        r(T), #reflux flowrate to each tray
        v,    #vapor boilup,
        mu
       };

  YVAR {yf(T), #existence of feed to each tray
        yr(T)  #existence of reflux to each tray
       };

  ZVAR {l(T), #liquid flowrate from each tray
        x(T1), #liquid composition on each tray
        y(T1), #vapor composition on each tray
        d,     #distillate flowrate
        b,     #bottoms flowrate
        obj    #objective integral
       };

  BINA {yf(T),yr(T)};
  POSI {f(T),r(T)};

  LBDS f(T) = <i E T| 0>;
  STP  f(T) = <i E T| 0 + (i==15)*1 >;
  UBDS f(T) = <i E T| 1>;

  LBDS r(T) = <i E T| 0>;
  STP  r(T) = <i E T| 0 + (i==30)*9.5 >;
  UBDS r(T) = <i E T| 100>;

  LBDS v = {0.1};
  STP  v = {10};
  UBDS v = {100};

  STP yf(T) = <i E T| 0 + (i==15)*1 >;
  STP yr(T) = <i E T| 0 + (i==30)*1 >;

  ISPE {obj, x'(T1)};

  TIME {0,400};
}}

MODEL {{
  MIN: mu;

#Objective (Integral of the Squared Error)
  obj:   obj' - (x(31)-xdset)^2 - (x(0)-xbset)^2 =e= 0;

#Reboiler Component Balance
  dae1:  10*m*x'(0) + v*y(0) + b*x(0) =e= l(1)*x(1);

#Tray component balances
  dae2(i E T & i<N ):  m*x'(i) + l(i)*x(i) + v*y(i)
                     =e= l(i+1)*x(i+1) + v*y(i-1) + f(i)*(0.54 - 0.09*exp[-10*t]) + r(i)*x(31);
  dae2(i E T & i==N):  m*x'(i) + l(i)*x(i) + v*y(i)
                     =e=                 v*y(i-1) + f(i)*(0.54 - 0.09*exp[-10*t]) + r(i)*x(31);

#Condenser Component balance
  dae3:  10*m*x'(31) + v*x(31) =e= v*y(30);

#Reboiler Total Balance
  dae4:  v + b =e= l(1);

#Tray total balances
  dae5(i E T & i<N ): l(i) =e= l(i+1) + f(i) + r(i);
  dae5(i E T & i==N): l(i) =e=          f(i) + r(i);

#Condenser Total balance
  dae8:  <<i E T| r(i)>> + d =e= v;

#Reboiler Equilibrium
  dae9:  alpha*x(0) =e= y(0)*(1+x(0)*(alpha - 1));

#Tray equilibrium
  dae10(i E T):  alpha*x(i) =e= y(i)*(1+x(i)*(alpha - 1));

  dae11: y(31) =e= 0;

#Point (end-time) Constraints
  pcon1[1]: mu - obj =g= 0;
  pcon2[1]: x(0) =l= xbset;
  pcon3[1]: x(31) =g= xdset;

#Logical constraints
  l1:  <<i E T| yf(i) >> =e= 1;
  l2:  <<i E T| yr(i) >> =e= 1;
  l3:  <<i E T| i*yr(i) >> =g= 2;
  l4:  <<i E T| i*yf(i) >> =g= 1;
  l5:  <<i E T| i*yr(i) >> =g= <<i E T| i*yf(i)>>+2;
#  l6:  v - <<it|r(it+1)>> =e= 5;
  l7(i E T):  f(i) =e=   1*yf(i);
  l8(i E T):  r(i) =l= 50*yr(i);
}}