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

OPTIONS {{
  EVTOL = "1e-10";
  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:4|; #reaction components
  SET J = |1:3|; #reaction paths
  SET K = |1:1|; #crossflow reactors (CFRs)
  SET L = |1:1|; #mixed reactors (CSTRs)
  SET R = |1|; #feeds
  SET P = |1|; #products

  SET CI = |0:9|;  #control intervals
  SET CP = |0:10|; #control points (intervals+1)

#Parameters
#reaction rate constants [1/s,1/s,L/gmol-s]
  PARAMETER k0(J) = {1,2,10};
#stoichiometric coefficients
  PARAMETER nu(I,J) = {-1, 0, -2,
                        1,-1,  0,
                        0, 1,  0,
                        0, 0,  1};

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

  XVAR {fab(R,K), #flowrate from feed splitters to CFR side mixers
        fac(R,K), #flowrate from feed splitters to CFR main mixers
        fad(R,L), #flowrate from feed splitters to CSTR mixers
        fb(K),    #flowrate of CFR side inlets
        cb(K,I),  #composition of CFR side inlets
        fc(K),    #flowrate of CFR main inlets
        cc(K,I),  #composition of CFR main inlets
        fd(L),    #flowrate of CSTR inlets
        cd(L,I),  #composition of CSTR inlets
        fe(K),    #flowrate of CFR main outlets
        ce(K,I),  #composition of CFR main outlets
        ff(K),    #flowrate of CFR side outlets
        cf(K,I),  #composition of CFR side outlets
        fg(L),    #flowrate of CSTR outlets
        cg(L,I),  #composition of CSTR outlets
        feb(K,K), #flowrate from CFR main splitters to CFR side mixers
        fec(K,K), #flowrate from CFR main splitters to CFR main mixers
        fed(K,L), #flowrate from CFR main splitters to CSTR mixers
        ffb(K,K), #flowrate from CFR side splitters to CFR side mixers
        ffc(K,K), #flowrate from CFR side splitters to CFR main mixers
        ffd(K,L), #flowrate from CFR side splitters to CSTR
        fgb(L,K), #flowrate from CSTR splitters to CFR side mixers
        fgc(L,K), #flowrate from CSTR splitters to CFR main mixers
        fgd(L,L), #flowrate from CSTR splitters to CSTR mixers
        feh(K,P), #flowrate from CFR main splitters to product mixers
        ffh(K,P), #flowrate from CFR side splitters to product 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

        vt(K),    #PFR reactor volumes
        ks(K,CP), #control parameters for dsdv
        ko(K,CP)  #control parameters for dldv
       };

#Bounds and starting points
  LBDS fab(R,K) = {0};
  STP  fab(R,K) = {0};
  UBDS fab(R,K) = {1000};

  LBDS fac(R,K) = {0};
  STP  fac(R,K) = {50};
  UBDS fac(R,K) = {1000};

  LBDS fad(R,L) = {0};
  STP  fad(R,L) = {50};
  UBDS fad(R,L) = {1000};

  LBDS fb(K) = {0};
  STP  fb(K) = {0};
  UBDS fb(K) = {100};

  LBDS cb(K,I) = {0};
  STP  cb(K,I) = {ca(R,I)};
  UBDS cb(K,I) = {1};

  LBDS fc(K) = {0};
  STP  fc(K) = {50};
  UBDS fc(K) = {1000};

  LBDS cc(K,I) = {0};
  STP  cc(K,I) = {ca(R,I)};
  UBDS cc(K,I) = {1};

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

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

  LBDS fe(K) = {0};
  STP  fe(K) = {50};
  UBDS fe(K) = {1000};

  LBDS ce(K,I) = {0};
  STP  ce(K,I) = {ca(R,I)};
  UBDS ce(K,I) = {1};

  LBDS ff(K) = {0};
  STP  ff(K) = {0};
  UBDS ff(K) = {1000};

  LBDS cf(K,I) = {0};
  STP  cf(K,I) = {ca(R,I)};
  UBDS cf(K,I) = {1};

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

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

  LBDS feb(K,K) = {0};
  STP  feb(K,K) = {0};
  UBDS feb(K,K) = {1000};

  LBDS fec(K,K) = {0};
  STP  fec(K,K) = {0};
  UBDS fec(K,K) = {1000};

  LBDS fed(K,L) = {0};
  STP  fed(K,L) = {0};
  UBDS fed(K,L) = {1000};

  LBDS ffb(K,K) = {0};
  STP  ffb(K,K) = {0};
  UBDS ffb(K,K) = {1000};

  LBDS ffc(K,K) = {0};
  STP  ffc(K,K) = {0};
  UBDS ffc(K,K) = {1000};

  LBDS ffd(K,L) = {0};
  STP  ffd(K,L) = {0};
  UBDS ffd(K,L) = {1000};

  LBDS fgb(L,K) = {0};
  STP  fgb(L,K) = {0};
  UBDS fgb(L,K) = {1000};

  LBDS fgc(L,K) = {0};
  STP  fgc(L,K) = {0};
  UBDS fgc(L,K) = {1000};

  LBDS fgd(L,L) = {0};
  STP  fgd(L,L) = {0};
  UBDS fgd(L,L) = {0};

  LBDS feh(K,P) = {0};
  STP  feh(K,P) = {50};
  UBDS feh(K,P) = {1000};

  LBDS ffh(K,P) = {0};
  STP  ffh(K,P) = {0};
  UBDS ffh(K,P) = {1000};

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

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

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

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

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

  LBDS vt(K) = {0};
  STP  vt(K) = {10};
  UBDS vt(K) = {1000};

  LBDS ks(K,CP) = {0};
  STP  ks(K,CP) = {0};
  UBDS ks(K,CP) = {1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,0};

  LBDS ko(K,CP) = {0};
  STP  ko(K,CP) = {0};
  UBDS ko(K,CP) = {0,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000};

  ZVAR {ft(K),   #CFR flowrate
        ct(K,I), #CFR concentration
        co(K,I), #CFR side exit concentration
        fs(K),   #CFR sidestream flowrate
        fo(K),   #CFR side exit flowrate
        rt(K,J), #CFR reaction rates
        dsdv(K),
        dodv(K)
       };

#States which specify the initial condition
  ISPE {ft(K), ct(K,I), co(K,I)};
#States whose initial conditions are to be optimized
  ICS  {ft(K), ct(K,I), co(K,I)};
#Parameters corresponding to the initial conditions
  ICP  {fc(K), cc(K,I), cc(K,I)};

#10 intervals
  PARA time(CP) = {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};
  TIME            {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};

}}

MODEL {{
  MIN: -ch(1,2);

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

#pre-reaction mixers
  mixbt(k E K):                fb(k) =e= <<r  E R|          fab(r,k) >>
                                       + <<k' E K|          feb(k',k) +          ffb(k',k)>>
                                       + <<l  E L| fgb(l,k) >>;
  mixbc(k E K, i E I): cb(k,i)*fb(k) =e= <<r  E R| ca(r,i) *fab(r,k) >>
                                       + <<k' E K| ce(k',i)*feb(k',k) + cf(k',i)*ffb(k',k)>>
                                       + <<l  E L| cg(l,i) *fgb(l,k) >>;

  mixct(k E K):                fc(k) =e= <<r  E R|          fac(r,k) >>
                                       + <<k' E K|          fec(k',k) +          ffc(k',k)>>
                                       + <<l  E L|          fgc(l,k) >>;
  mixcc(k E K, i E I): cc(k,i)*fc(k) =e= <<r  E R| ca(r,i) *fac(r,k) >>
                                       + <<k' E K| ce(k',i)*fec(k',k) + cf(k',i)*ffc(k',k)>>
                                       + <<l  E L| cg(l,i) *fgc(l,k) >>;

  mixdt(l E L):                fd(l) =e= <<r  E R|          fad(r,l) >>
                                       + <<k  E K|          fed(k,l) +         ffd(k,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) >>
                                       + <<k  E K| ce(k,i) *fed(k,l) + cf(k,i)*ffd(k,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,2);
  rate3(l E L): rm(l,3) =e= k0(3)*cg(l,1)*cg(l,1);

#CFRs
  cfrt(k E K): ft'(k) =e= -dsdv(k) - dodv(k);
  cftc(k E K, i E I): ct'(k,i)*(ft(k)+1e-6) =e= vt(k)*<<j E J|nu(i,j)*rt(k,j)>>
                                              - (cb(k,i)-ct(k,i))*dsdv(k);
  point1(k E K)[10]:                fe(k) =e=         ft(k);
  point2(k E K, i E I)[10]: ce(k,i)*fe(k) =e= ct(k,i)*ft(k);

#CFR reaction rates
  rate1(k E K): rt(k,1) =e= k0(1)*ct(k,1);
  rate2(k E K): rt(k,2) =e= k0(2)*ct(k,2);
  rate3(k E K): rt(k,3) =e= k0(3)*ct(k,1)*ct(k,1);

#CFR side feed
  mat1(k E K):  fb(k) =e= ks(k,0);
  control1(k E K):  fs(k) =e= interv[i E CI| ks(k,i) + (ks(k,i+1)
                            - ks(k,i))*(t-time(i))/(time(i+1)-time(i))];
  dcontrol1(k E K): dsdv(k) =e= interv[i E CI|(ks(k,i+1)- ks(k,i))/
                                                   (time(i+1)-time(i))];

#CFR side exit
  exitc(k E K, i E I): co'(k,i)*(fo(k)+1e-6) =e= (ct(k,i)-co(k,i))*dodv(k);
  mat2(k E K):  ff(k) =e= ko(k,10);
  point3(k E K, i E I)[10]: cf(k,i)*ff(k) =e= co(k,i)*fo(k);

  control2(k E K):  fo(k) =e= interv[i E CI| ko(k,i) + (ko(k,i+1)
                            - ko(k,i))*(t-time(i))/(time(i+1)-time(i))];
  dcontrol2(k E K): dodv(k) =e= interv[i E CI|(ko(k,i+1)- ko(k,i))/
                                                   (time(i+1)-time(i))];

#post-reactor splitters
  splite(k E K):  fe(k) =e= <<k' E K| feb(k,k') + fec(k,k')>>
                          + <<l  E L| fed(k,l) >>
                          + <<p  E P| feh(k,p) >>;

  splitf(k E K): ff(k) =e= <<k' E K| ffb(k,k') + ffc(k,k')>>
                         + <<l  E L|  ffd(k,l) >>
                         + <<p  E P|  ffh(k,p) >>;

  splitg(l E L):  fg(l) =e= <<k  E K| fgb(l,k) + fgc(l,k) >>
                          + <<l' E L| fgd(l,l')>>
                          + <<p  E P| fgh(l,p) >>;

#product mixer
  mixht(p E P):                fh(p) =e= <<k E K|         feh(k,p) +         ffh(k,p)>>
                                       + <<l E L|         fgh(l,p)>>;

  mixhc(p E P, i E I): ch(p,i)*fh(p) =e= <<k E K| ce(k,i)*feh(k,p) + cf(k,i)*ffh(k,p)>>
                                       + <<l E L| cg(l,i)*fgh(l,p)>>;

#logical constraints
  log1(k E K, i E CI): ks(k,i+1) =l= ks(k,i);
  log2(k E K, i E CI): ko(k,i+1) =g= ko(k,i);
}}