$*************************************************************
$ Multiperiod MILP Model for Planning of Chemical Complexes
$ N. V. Sahinidis and I. E. Grossmann
$
$ CACHE Process Design Case Studies
$ M. Morari and I.E. Grossmann
$
$ Optimal Solution:  -1697.61
$*************************************************************

DECLARATION {{
  INDEX {i,j,k,l,t,tt,m};

  SET P = |1:3|; # Processes
  SET C = |1:3|; # Chemicals
  SET K = |1:1|; # Markets
  SET T = |1:3|; # Time Periods

  SET TT = |1:100|; #years

  SET L(P) = {|1,2|, |2,3|, |2,3|}; #inputs and outputs of the processes
  SET M(P) = {|2|,|3|,|3|}; #main products of the processes
  SET I(C) = {|1|, |2,3|, |.|}; #Plants which consume each chemical
  SET O(C) = {|.|, |1|, |2,3|}; #Plants which produce each chemical


  PARA np = {3}; #Number of processes
  PARA nc = {3}; #Number of chemicals
  PARA nm = {9}; #Number of markets
  PARA nt = {3}; #Number of time periods

  PARA lenp(T)   = {2,3,5};          #length of the time periods [years]
  PARA ci(T)     = {200,300,400};    #capital available for investment [1e5 $/yr]
  PARA qelb(P,T) = <i E P|<t E T|   0 > >; #lower bound on expansion
  PARA qeub(P,T) = <i E P|<t E T| 100 > >; #upper bound on expansion
  PARA nexp(P)   = {3,3,3};          #maximum number of expansions
  PARA excap(P)  = {0,50,0};         #existing capacities [kton/yr]
  PARA svalf(P)  = {0.1,0.1,0.1};    #salvage value factor
  PARA wcapf(P)  = {0.15,0.15,0.15}; #working capital factor

  PARA int = 0.10;
  PARA tax = 0.45;

#variable investment coefficient [1e2 $/ton]
  PARA alpha(P,T) = {1.38,1.56,1.78,
                     2.72,3.22,4.6 ,
                     1.76,2.34,2.84};

#fixed investment coefficient [1e5 $/ton]
  PARA beta(P,T) = { 85, 95,112,
                     73, 82,102,
                    110,125,148};

#operating expenses per unit of main product
  PARA delta(P,T) = {0.4,0.5,0.6,
                     0.6,0.7,0.8,
                     0.5,0.6,0.7};

#cost for purchase of one unit of chemical [1e2 $/ton]
  PARA lam(C,K,T) = <j E C|<k E K|<t E T| 0 > > >; 
  PARA lam(1,1,T) = {4  , 5.24, 7.32};
  PARA lam(2,1,T) = {9.6,11.52,13.52};

#purchase lower bounds
  PARA plb(C,K,T) = <j E C|<k E K|<t E T| 0 > > >;

#purchase upper bounds
  PARA pub(C,K,T) = <j E C|<k E K|<t E T| 0 > > >;
  PARA pub(1,1,T) = { 6, 7.5, 8.6};
  PARA pub(2,1,T) = {20,25.5,30  };

#sale prices
  PARA gamma(C,K,T) = <j E C|<k E K|<t E T| 0 > > >;
  PARA gamma(3,1,T) = {26.20,29.20,35.20};

#sales lower bounds
  PARA slb(C,K,T) = <j E C|<k E K|<t E T| 0 > > >;

#sales upper bounds
  PARA sub(C,K,T) = <j E C|<k E K|<t E T| 0 > > >;
  PARA sub(3,1,T) = {65,75,90};

#mass balance coefficients
  PARA mu(P,C) = {1.11,1   ,0,
                  0   ,1.22,1,
                  0   ,1.05,1};


#Determine the discounted values of the parameters
  PARA x(T) = <t E T| 0 + <<k E T & k<t| lenp(k)>> >;

  PARA fx = <<t E T| lenp(t) >>;

  PARA int1 = 1/(1+int);

  PARA c2(T) = <t E T| <<tt E TT & (tt >= x(t)+1) && (tt <= fx)| int1^tt >> >;

  PARA t1(T) = <t E T| -int1^x(t) >;

  PARA t2(P,T) = <i E P| <t E T| wcapf(i)*(t1(t)+int1^fx) > >;

  PARA t3(P,T) = <i E P| <t E T| (1-svalf(i))*tax*c2(t)/(fx-x(t)) > >;

  PARA t4(P) = <i E P| svalf(i)*int1^fx >;

  PARA c1(P,T) = <i E P| <t E T| t1(t)+t2(i,t)+t3(i,t)+t4(i) > >;

  PARA alphar(P,T) = <i E P| <t E T| -alpha(i,t)*c1(i,t) > >;
  PARA betar(P,T)  = <i E P| <t E T| -beta(i,t) *c1(i,t) > >;

  SET T2 = |1:2|;
  PARA scoeff(T) = <t E T| 0>;
  PARA scoeff(T2) = <t E T & t<card[T] | <<tt E TT & (tt>=(x(t)+1)) && (tt<=x(t+1))|
                             int1^tt >>*(1-tax) >;
  PARA scoeff(3) = <<tt E TT & (tt>=x(3)+1) && (tt<=fx)|
                             int1^tt >>*(1-tax);

  PARA lam(C,K,T) = <j E C| <k E K| <t E T| scoeff(t)*lam(j,k,t) > > >;
  PARA gamma(C,K,T) = <j E C| <k E K| <t E T| scoeff(t)*gamma(j,k,t) > > >;
  PARA delta(P,T) = <i E P| <t E T| scoeff(t)*delta(i,t) > >;

  XVAR {q(P,T),   #capacities
        qe(P,T),  #capacity expansions
        y(P,T),   #integer decision variables
        w(P,C,T), #flow rates
        s(C,K,T), #Sales
        p(C,K,T)  #Purchases
       };

  POSI {q(P,T),qe(P,T),w(P,C,T),s(C,K,T),p(C,K,T)};
  BINA {y(P,T)};

  LBDS p(C,K,T) = {plb(C,K,T)};
  UBDS p(C,K,T) = {pub(C,K,T)};

  LBDS s(C,K,T) = {slb(C,K,T)};
  UBDS s(C,K,T) = {sub(C,K,T)};
}}

MODEL {{
  MIN:  <<i E P| <<t E T| alphar(i,t)*qe(i,t) + betar(i,t)*y(i,t) >> >>
      + <<i E P| <<t E T| delta(i,t)*<<m E M(i)| w(i,m,t)>> >> >>
      - <<j E C| <<l E K| <<t E T| gamma(j,l,t)*s(j,l,t) - lam(j,l,t)*p(j,l,t) >> >> >>;

  con1(i E P, t E T): qe(i,t) =g= y(i,t)*qelb(i,t);
  con2(i E P, t E T): qe(i,t) =l= y(i,t)*qeub(i,t);

  cap(i E P, t E T & t==1): q(i,t) =e= excap(i) + qe(i,t);
  cap(i E P, t E T & t>1 ): q(i,t) =e= q(i,t-1) + qe(i,t);

  prcp(i E P, t E T): <<m E M(i)| w(i,m,t) >> =l= q(i,t);

  pbal(i E P, j E L(i), t E T): w(i,j,t) =e= mu(i,j)*<<m E M(i)| w(i,m,t) >>;

  cbal(j E C, t E T): <<l E K| p(j,l,t)>> + <<k E O(j)| w(k,j,t)>> =e=
                      <<l E K| s(j,l,t)>> + <<k E I(j)| w(k,j,t)>>;

  limit(i E P):  <<t E T| y(i,t)>> =l= nexp(i);

  limit(t E T):  <<i E P| alpha(i,t)*qe(i,t) + beta(i,t)*y(i,t) >> =l= ci(t);
}}