#Isothermal Reactor Network Synthesis
#Constant Density Assumption
#just CSTRs for algebraic formulation
#
# C. A. Schweiger and C. A. Floudas, "Optimization
#   Framework for the Synthesis of Chemical Reactor
#   Networks", Ind. Eng. Chem. Res. 1998, submitted
#

OPTIONS {{
  FTOL = "1e-6";
  OTOL = "1e-6";
  SNOPT = "Iterations Limit 100000";
  SNOPT = "Major Iterations Limit 10000";
  MINOS = "Iterations Limit 100000";
  MINOS = "Major Iterations Limit 10000";
}}

DECLARATIONS {{
  INDEX {i,j,k,k',l,l',r,p};
#Sets
  SET I = |1:5|; #reaction components
  SET J = |1:7|; #reaction paths
  SET L = |1:5|; #units (reactors, etc.)
  SET R = |1|; #feeds
  SET P = |1|; #products

#Parameters
#reaction rate constants [1/s,1/s,1/s,L/gmol-s,L/gmol-s,1/s,1/s]
  PARAMETER k0(J) = {0.33384,0.26687,0.14940,0.18957,0.009598,0.29425,0.011932};
#stoichiometric coefficients
  PARAMETER nu(I,J) = {-1, -1,  0,  0, -2,  0,  0,
                        0,  0,  1,  0,  0, -1,  0,
                        0,  0,  0,  1,  1,  0, -1,
                        0,  1, -1, -2,  0,  1,  2,
                        1,  0,  0,  0,  0,  0,  0};

#feed flowrate [L/s]
  PARAMETER fa(R) = {100};
#feed composition [gmol/L]
  PARAMETER ca(R,I) = {1,0,0,0,0};

  XVAR {fad(R,L), #flowrate from feed splitters to CSTR mixers
        fd(L),    #flowrate of CSTR inlets
        cd(L,I),  #composition of CSTR inlets
        fg(L),    #flowrate of CSTR outlets
        cg(L,I),  #composition of CSTR outlets
        fgd(L,L), #flowrate from CSTR splitters to CSTR mixers
        fgh(L,P), #flowrate from CSTR splitters to product mixers
        fh(P),    #flowrate of product
        ch(P,I),  #composition of product

        vm(L),    #CSTR reactor volumes
        rm(L,J)   #CSTR reaction rates
       };

#Bounds and starting points
  LBDS fad(R,L) = {0};
  STP  fad(R,L) = {100,0,0,0,0};
  UBDS fad(R,L) = {1000};

  LBDS fd(L) = {0};
  STP  fd(L) = {100};
  UBDS fd(L) = {1000};

  LBDS cd(L,I) = {0};
  STP  cd(L,I) = {ca(R,I)};
  UBDS cd(L,I) = {1};

  LBDS fg(L) = {0};
  STP  fg(L) = {100};
  UBDS fg(L) = {1000};

  LBDS cg(L,I) = {0};
  STP  cg(L,I) = {ca(R,I)};
  UBDS cg(L,I) = {1};

  LBDS fgd(L,L) = {0};
  STP  fgd(L,L) = <l E L| <l' E L| (l == l'-1)*100 > >;
  UBDS fgd(L,L) = <l E L| <l' E L| (l != l')*100 > >;

  LBDS fgh(L,P) = {0};
  STP  fgh(L,P) = {0,0,0,0,100};
  UBDS fgh(L,P) = {1000};

  LBDS fh(P) = {0};
  STP  fh(P) = {100};
  UBDS fh(P) = {1000};

  LBDS ch(P,I) = {0,0,0,0.01,0};
  STP  ch(P,I) = {ca(R,I)};
  UBDS ch(P,I) = {1};

  LBDS vm(L) = {0};
  STP  vm(L) = {1200};
  UBDS vm(L) = {10000};

  LBDS rm(L,J) = {0};
  STP  rm(L,J) = {0};
  UBDS rm(L,J) = {10000};
}}

CONS {{
  MIN: -ch(1,3)/ch(1,4);

#feed splitter
  splita(r E R): fa(r) =e= <<l E L| fad(r,l)>>;

#pre-reaction mixers
  mixdt(l E L):                fd(l) =e= <<r  E R|          fad(r,l) >>
                                       + <<l' E L|          fgd(l',l)>>;
  mixdc(l E L, i E I): cd(l,i)*fd(l) =e= <<r  E R| ca(r,i) *fad(r,l) >>
                                       + <<l' E L| cg(l',i)*fgd(l',l)>>;

#CSTRs
  cstrt(l E L):                fg(l) =e=         fd(l);
  cstrc(l E L, i E I): cg(l,i)*fg(l) =e= cd(l,i)*fd(l) + vm(l)*<<j E J|nu(i,j)*rm(l,j)>>;

#CSTR reaction rates
  rate1(l E L): rm(l,1) =e= k0(1)*cg(l,1);
  rate2(l E L): rm(l,2) =e= k0(2)*cg(l,1);
  rate3(l E L): rm(l,3) =e= k0(3)*cg(l,4);
  rate4(l E L): rm(l,4) =e= k0(4)*cg(l,4)*cg(l,4);
  rate5(l E L): rm(l,5) =e= k0(5)*cg(l,1)*cg(l,1);
  rate6(l E L): rm(l,6) =e= k0(6)*cg(l,2);
  rate7(l E L): rm(l,7) =e= k0(7)*cg(l,3);

#post-reactor splitters
  splitg(l E L): fg(l) =e= <<l' E L| fgd(l,l')>> + <<p E P|fgh(l,p)>>;

#product mixer
  mixht(p E P):        fh(p)         =e= <<l E L| fgh(l,p)         >>;
  mixhc(p E P, i E I): fh(p)*ch(p,i) =e= <<l E L| fgh(l,p)*cg(l,i) >>;

#Volume constraint
  vol: <<l E L| vm(l) >> =l= 6000;
}}