$*************************************************************
$ Optimum Feed Plate Location
$ J. Viswanathan and Ignacio E. Grossmann
$
$ CACHE Process Design Case Studies
$ M. Morari and I.E. Grossmann
$
$ Optimal Solution:  -13.406
$*************************************************************

OPTIONS {{
  AUTOINIT;
  SNOPT = "Iterations Limit 1000000";
  SNOPT = "Major Iterations Limit 10000";
}}

DECLARATION {{
  INDEX {i, j, k};
  SET I = |1:9|; #Trays including condenser and reboiler
  SET R = |1|;   #Reboiler
  SET S = |2:8|; #Trays (for feed location)
  SET C = |9|;   #Condenser
  SET J = |1:2|; #components (Benzene-Toluene)

  SET K = |1:5|;
  SET K2 = |1:4|;

#Constants
#Gas constant
  PARA rg = {8.314};
#molar mass
  PARA mm(J) = {78.114, 92.141};
#Normal boiling point [K]
  PARA Tb(J) = {353.2, 383.8};
#Critical Termperature [K]
  PARA Tc(J) = {562.2, 591.8};
#Critical Pressure [bar]
  PARA Pc(J) = {48.9, 41.0};
#Omega
  PARA omega(J) = {0.212, 0.263};
#Liquid density
  PARA rho(J) = {0.885, 0.867};
#
  PARA Tden(J) = {289, 293};

#Constants for the vapor pressure
  PARA vpcon(J,K) = {-6.98273, 1.33213, -2.62863, -3.33399, 288,
                     -7.28607, 1.38091, -2.83433, -2.79168, 309};

#Constants for the isobaric heat capacity equation
  PARA Cp(J,K2) = {-3.392e1, 4.739e-1, -3.017e-4, 7.130e-8,
                   -2.435e1, 5.125e-1, -2.765e-4, 4.911e-8};

#
  PARA Tcon = {355};
#Integration constant for enthalpy
  PARA h0(J) = <j E J| <<k E K2| Tcon^k*(Cp(j,k)/k) >> >;

#Feed specifications
  PARA F = {100};    #[moles]
  PARA pf = {1.12};  #[bar]
  PARA Tf = {359.6}; #[k]
  PARA Tf = {363.368}; #K
  PARA zf(J) = {0.70, 0.30};

#Pressure on each tray [bar]
  PARA preb = 1.20;
  PARA pbot = 1.12;
  PARA ptop = 1.08;
  PARA pcon = 1.05;
  PARA p(I) = <i E I| (i==1)*preb + (i>1 && i<9)*(pbot - ((pbot-ptop)/6*(i-2))) + (i==9)*pcon >;

  PARA hf = <<j E J| zf(j)*( <<k E K2| Tf^k*Cp(j,k)/k >> - h0(j)
                + rg*Tc(j)*( vpcon(j,1)*(1-Tf/Tc(j))
                           + vpcon(j,2)*(1-Tf/Tc(j))^1.5
                           + vpcon(j,3)*(1-Tf/Tc(j))^3
                           + vpcon(j,4)*(1-Tf/Tc(j))^6 )
                + rg*Tf*(   vpcon(j,1)
                        + 1.5*vpcon(j,2)*(1-Tf/Tc(j))^0.5
                        +   3*vpcon(j,3)*(1-Tf/Tc(j))^2
                        +   6*vpcon(j,4)*(1-Tf/Tc(j))^5 ) ) >> ;

  XVAR {x(I,J),  #mole fraction of the jth component in liquid on ith tray
        y(I,J),  #mole fraction of the jth component in vapor on ith tray
        L(I),    #molar flow rate of liquid leaving tray i
        V(I),    #molar flow rate of vapor leaving tray i
        T(I),    #temperature of tray i
        feed(I), #molar flow rate of feed entering tray i
        r,       #reflux ratio
        P1,      #molar flow rate of top product
        P2,      #molar flow rate of bottom product
        hl(I),   #molar enthalpy of liquid on tray i
        hv(I)    #molar enthalpy of vapor on tray i
       };

  YVAR {q(I)}; #existence of tray i
  BINA {q(I)};

  LBDS x(I,J) = <i E I|< j E J| 0 > >;
  LBDS x(C,1) = {0.95}; #Product purity
  STP  x(I,J) = {0.400000, 0.6,
                 0.462500, 0.537500,
                 0.525000, 0.475000,
                 0.587500, 0.412500,
                 0.650000, 0.350000,
                 0.712500, 0.287500,
                 0.775000, 0.225000,
                 0.837500, 0.162500,
                 0.900000, 0.10000};
  UBDS x(I,J) = <i E I| <j E J| 1 > >;

  LBDS y(I,J) = <i E I| <j E J| 0 > >;
  STP  y(I,J) = {0.2, 0.8,
                 0.3, 0.7,
                 0.4, 0.6,
                 0.5, 0.5,
                 0.6, 0.4,
                 0.7, 0.3,
                 0.8, 0.2,
                 0.9, 0.1,
                 1.0, 0.0};
  UBDS y(I,J) = <i E I| <j E J| 1 > >;

  LBDS L(I) = <i E I| 0 >;
  STP  L(I) = {40,127,127,127,127,127,27,27,27};
  UBDS L(I) = <i E I| 1000 >;

  LBDS V(I) = <i E I| 0 >;
  STP  V(I) = {87,87,87,87,87,87,87,87,87};
  UBDS V(I) = <i E I| 1000 >;

  LBDS T(I) = <i E I| 345 >;
  STP  T(I) = <i E I| 375 >;
  UBDS T(I) = <i E I| 390 >;

  LBDS feed(I) = <i E I| 0 >;
  STP  feed(I) = {0,0,0,0,0,100,0,0,0};
  UBDS feed(I) = <i E I| F >;

  LBDS r = {0.1};
  STP  r = {0.45};
  UBDS r = {0.95};

  LBDS P1 = {50};
  STP  P1 = {60};
  UBDS P1 = {80};

  LBDS P2 = {20};
  STP  P2 = {40};
  UBDS P2 = {50};

  LBDS hl(I) = <i E I| -60000 >;
  STP  hl(I) = <i E I| -30000 >;
  UBDS hl(I) = <i E I| -10000 >;

  LBDS hv(I) = <i E I| 0 >;
  STP  hv(I) = <i E I| 2000 >;
  UBDS hv(I) = <i E I| 6000 >;

}}

MODEL {{
  MIN: 50*r - P1;

#Phase equilibrium
  phase(i E I, j E J): y(i,j)*p(i) - x(i,j)*Pc(j)*exp[Tc(j)/T(i)* 
                      ( vpcon(j,1)*(1-T(i)/Tc(j))
                      + vpcon(j,2)*(1-T(i)/Tc(j))^1.5
                      + vpcon(j,3)*(1-T(i)/Tc(j))^3 
                      + vpcon(j,4)*(1-T(i)/Tc(j))^6   )] =e= 0;

  errk(i E I): <<j E J|y(i,j)>> =e= <<j E J|x(i,j)>>;

#Component material balances
  cmbr(i E R, j E J): L(i)*x(i,j) + V(i)*y(i,j) =e= L(i+1)*x(i+1,j);
  cmbt(i E S, j E J): L(i)*x(i,j) + V(i)*y(i,j) =e= L(i+1)*x(i+1,j) + V(i-1)*y(i-1,j)
                                                     + feed(i)*zf(j);
  cmbc(i E C, j E J): L(i)*x(i,j) + P1*x(i,j) =e= V(i-1)*y(i-1,j);

#total material balances
  tmbr(i E R): L(i) + V(i) =e= L(i+1);
  tmbt(i E S): L(i) + V(i) =e= L(i+1) + V(i-1) + feed(i);
  tmbc(i E C): L(i) + P1   =e=          V(i-1);

  efln(i E C):  L(i) =e= r*P1;
  efp2(i E R):  L(i) =e=   P2;

#Equilibrium expressions
  efhl(i E I):  hl(i) - <<j E J| x(i,j)*( <<k E K2| T(i)^k*Cp(j,k)/k >> - h0(j)
                + rg*Tc(j)*( vpcon(j,1)*(1-T(i)/Tc(j))
                           + vpcon(j,2)*(1-T(i)/Tc(j))^1.5
                           + vpcon(j,3)*(1-T(i)/Tc(j))^3
                           + vpcon(j,4)*(1-T(i)/Tc(j))^6 )
                + rg*T(i)*(   vpcon(j,1)
                        + 1.5*vpcon(j,2)*(1-T(i)/Tc(j))^0.5
                        +   3*vpcon(j,3)*(1-T(i)/Tc(j))^2
                        +   6*vpcon(j,4)*(1-T(i)/Tc(j))^5 ) ) >> =e= 0;

  efhv(i E I): hv(i) - <<j E J| y(i,j)*(<<k E K2|T(i)^k*Cp(j,k)/k >> - h0(j) ) >> =e= 0;

#Enthalpy balance
  eb(i E S): L(i)*hl(i) + V(i)*hv(i) =e= L(i+1)*hl(i+1) + V(i-1)*hv(i-1)
                                          + feed(i)*hf;

  sumf:<<i E S| feed(i) >> =e= F;

#Feed location constraints
  sumb:<<i E S| q(i) >> =e= 1;
  lfeed(i E S): feed(i) =l= F*q(i);
  lex1:q(1) =e= 0;
  lex2:q(9) =e= 0;
}}