$*************************************************************
$ A. Aggarwal and C.A. Floudas, "Synthesis of General
$   Distillation Sequences: Nonsharp Separation", Comput.
$   Chem. Eng., 1990, 14 (6), 631--653.
$
$ Example 2
$
$ Optimal Solution:  62596.7
$*************************************************************

OPTIONS{{
  FTOL = 1e-5;
  OTOL = 1e-5;
  RFIX; BALLR = "1.0"; BALLS = "1.0"; RFIXBOTH; BALLT = "1e+10";
#  RETRY = 3;
}}

VARIABLES{{
  INDEX {i,j,k,l,m,n,p,q,s};

  SET J =  |1:2|;   #columns
  SET I =  |1:3|;   #components
  SET N =  |0:20|;  #streams
  SET S =  |0:4|;   #splitters
  SET Mc = |1:2|;   #mixer prior to each column
  SET Mf = |1:2|;   #mixer prior to each product
  SET P =  |1:2|;   #products

#This is to reduce the number of variables defined
  SET Ncf = |5,13|; #feed streams to the columns

  SET Sin(S) = {|0|,|6|, |9|, |14|,|18|};
  SET Sout(S) = {|1,2,3,4|, |7,8|, |10,11,12|, |15,16,17|, |19,20|};
  SET Mcin(Mc) = {|1,15|, |2,10|};
  SET Mcout(Mc) = {|5|, |13|};
  SET Mfin(Mf) = {|3,7,11,16,19|, |4,8,12,17,20|};
  SET SU(J) = {|5|,|13|};
  SET SUtop(J) = {|6|, |14|};
  SET SUbot(J) = {|9|, |18|};
  SET LK(J) = {|1|, |2|};
  SET HK(J) = {|2|, |3|};
  SET ND(N) = {|.|, |.|, |.|, |.|, |.|, |.|, |3|, |.|, |.|, |.|, |.|,
                    |.|, |.|, |.|, |.|, |.|, |.|, |.|, |1|, |.|, |.|};

  PARA a0(J) = {23947.0,75835.0};
  PARA a1(J) = {-1399.04,-6615.88};
  PARA a2(J) = {935.14,3381.47};
  PARA a3(J) = {773.08,3733.49};
  PARA b(I,J) = {-57.19,163.71,
                 426.56,2889.96,
                   0   ,   0   };
  PARA C(I,P) = {80,20,
                 30,70,
                 20,80};

  PARA U = 300;

  XVAR {x(I,N), rlk(I,J), rhk(I,J)};
  YVAR {y(J), F(N), f(I,Ncf)};
  BINA {y(J)};

  LBDS x(I,N) = <i E I| <n E N| 0> >;
  UBDS x(I,N) = <i E I| <n E N| 1> >;
  LBDS x(I,0) = <i E I| 1/3>;
  UBDS x(I,0) = <i E I| 1/3>;


  LBDS rhk(I,J)  = {0.85, 0.85, 0.85, 0.85, 0.85, 0.85};
  UBDS rhk(I,J)  = {1, 1, 1, 1, 1, 1};

  LBDS rlk(I,J)  = {0.85, 0.85, 0.85, 0.85, 0.85, 0.85};
  UBDS rlk(I,J)  = {1, 1, 1, 1, 1, 1};

  STP y(J)  = {1, 1};

  LBDS F(N) = {300,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0};
  UBDS F(N) = {300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300};

  LBDS f(I,Ncf) = {0,  0,
                   0,  0,
                   0,  0};
  UBDS f(I,Ncf) = {100,100,
                   100,100,
                   100,100};

#Columns 1 and 2 in series
#  STP  F(N) = {300,300,  0,0,0,300,100,80,20,200,200, 0, 0,200,100,  0,30,70,100,20,80};
#  STP  f(I,Ncf) = {100,  0,
#                   100,100,
#                   100,100};

#Columns 2 and 1 in series
  STP  F(N) = {300,0,300,0,0,200,100,80,20,100,  0,30,70,300,200,200, 0, 0,100,20,80};
  STP  f(I,Ncf) = {100,100,
                   100,100,
                     0,100};
}}

MODEL {{
  MIN: <<j E J| a0(j)*y(j) +
               <<k E SU(j)| (a1(j) + <<i E LK(j)| a2(j) *rlk(i,j) >> 
                                   + <<i E HK(j)| a3(j) *rhk(i,j) >> 
                                   + <<i E I    | b(i,j)*x(i,k) >> )*F(k) >> >>;

  split(s E S): <<l E Sin(s)| F(l)>> =e= <<l E Sout(s)| F(l)>>;
  keyl(j E J, i E LK(j), n E SU(j), p E SUtop(j)): F(p)*x(i,p) - rlk(i,j)*f(i,n) =e= 0;
  keyh(j E J, i E HK(j), n E SU(j), q E SUbot(j)): F(q)*x(i,q) - rhk(i,j)*f(i,n) =e= 0;
  flow(j E J, i E I, n E SU(j)): f(i,n) =e= F(n)*x(i,n);
  comp(i E I, j E J, n E SU(j), p E SUtop(j), q E SUbot(j)): f(i,n) - F(p)*x(i,p) - F(q)*x(i,q) =e= 0;
  mixc(i E I, m E Mc): <<l E Mcin(m)| F(l)*x(i,l)>> - <<l E Mcout(m)| f(i,l)>> =e= 0;
  nonk(k E N, i E ND(k)): x(i,k) =e= 0;
  mixf(i E I, p E P): <<l E Mfin(p)| F(l)*x(i,l)>> =e= C(i,p);
  sumx(s E S, k E Sin(s)): <<i E I| x(i,k) >> =e= 1;
  sumx(m E Mc, k E Mcout(m)): <<i E I| x(i,k) >> =e= 1;
  logi(j E J, k E SU(j)): F(k) - U*y(j) =l= 0;
  matx(s E S, n E Sout(s), i E I): x(i,n) =e= <<m E Sin(s)| x(i,m)>>;
}}