#Nonisothermal Reactor Network Synthesis
#Constant Density Assumption
#CSTRs and CFRs
#linear control of sidestream flowrates
#quadratic control of CFR temperature with continuous 1st derivatives
#
# 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:2|; #reaction paths
  SET K = |1:1|; #crossflow reactors (CFRs)
  SET L = |1:1|; #mixed reactors (CSTRs)
  SET R = |1:2|; #feeds
  SET P = |1|; #products

  SET CI = |0:9|;  #control intervals
  SET CP = |0:10|; #control points (intervals+1)
  SET CO = |0:2|;  #collocation points
  SET OR = |0:1|;  #order or the collocation equations

#Parameters
#Gas constant
  PARAMETER rg = 1.987; #cal/gmol/K
#reaction rate constants []
  PARAMETER k0(J) = {5.4e7,3.6e5};
#activation energy []
  PARAMETER e(J) = {19138,9569};
#heat of reaction (-H/rho/Cp) []
  PARAMETER h(J)  = {10,20};
#stoichiometric coefficients
  PARAMETER nu(I,J) = {-1,-1,
                       -1,-1,
                        1, 0,
                        0, 1};

#feed flowrate []
  PARAMETER fa(R) = {50,50};
#feed composition []
  PARAMETER ca(R,I) = {1,0,0,0,
                       0,1,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
        tm(L),    #CSTR reactor temperatures
        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
        kt(K,CP,OR) #control parameters for dtdv
       };

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

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

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

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

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

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

  LBDS cc(K,I) = {1e-6,1e-6,0,0};
  STP  cc(K,I) = {0.5,0.5,0,0};
  UBDS cc(K,I) = {1};

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

  LBDS cd(L,I) = {0};
  STP  cd(L,I) = {0.5,0.5,0,0};
  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) = {0.5,0.5,0,0};
  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) = {0.5,0.5,0,0};
  UBDS cf(K,I) = {1};

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

  LBDS cg(L,I) = {1e-6,1e-6,0,0};
  STP  cg(L,I) = {0.5,0.5,0,0};
  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};
  UBDS fec(K,K) = {0};

  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) = {0.5,0.5,0,0};
  UBDS ch(P,I) = {1};

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

  LBDS tm(L) = {450};
  STP  tm(L) = {800};
  UBDS tm(L) = {800};

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

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

  LBDS ks(K,CP) = {0};
  STP  ks(K,CP) = {50,45,40,35,30,25,20,15,10,5,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};

  LBDS kt(K,CP,OR) = {450};
  STP  kt(K,CP,OR) = {800};
  UBDS kt(K,CP,OR) = {800};

  ZVAR {ft(K),   #CFR flowrate
        ct(K,I), #CFR concentration
        tt(K),   #CFR temperature
        co(K,I), #CFR side exit concentration
        fs(K),   #CFR sidestream flowrate
        fo(K),   #CFR side exit flowrate
        rt(K,J), #CFR reaction rates
        dtdv(K),
        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)};
#Parameters corresponding to the initial conditions
  ICP  {fc(K), cc(K,I)};
#Initial conditions of remaining z-varaibles
  STP co(K,I) = {0};

#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};

#collocation points
  PARA pt(CI,CO) = {0.00,0.05,0.10,
                    0.10,0.15,0.20,
                    0.20,0.25,0.30,
                    0.30,0.35,0.40,
                    0.40,0.45,0.50,
                    0.50,0.55,0.60,
                    0.60,0.65,0.70,
                    0.70,0.75,0.80,
                    0.80,0.85,0.90,
                    0.90,0.95,1.00};
}}

MODEL {{
  MIN: -100*ch(1,3) + <<l E L|vm(l)>> + <<k E K| vt(k)>>;

#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)*exp[-e(1)/rg/tm(l)]*cg(l,1)*cg(l,2)^0.3;
  rate2(l E L): rm(l,2) =e= k0(2)*exp[-e(2)/rg/tm(l)]*cg(l,1)^0.5*cg(l,2)^1.8;

#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);
  control3(k E K): tt(k) =e= interv[i E CI|
    kt(k,i  ,0)*(t-pt(i,1))*(t-pt(i,2))/((pt(i,0)-pt(i,1))*(pt(i,0)-pt(i,2)))
  + kt(k,i  ,1)*(t-pt(i,0))*(t-pt(i,2))/((pt(i,1)-pt(i,0))*(pt(i,1)-pt(i,2)))
  + kt(k,i+1,0)*(t-pt(i,0))*(t-pt(i,1))/((pt(i,2)-pt(i,0))*(pt(i,2)-pt(i,1))) ];

  dcontrol3(k E K): dtdv(k) =e= interv[i E CI|
    kt(k,i  ,0)*((t-pt(i,1))+(t-pt(i,2)))/(pt(i,0)-pt(i,1))/(pt(i,0)-pt(i,2))
  + kt(k,i  ,1)*((t-pt(i,0))+(t-pt(i,2)))/(pt(i,1)-pt(i,0))/(pt(i,1)-pt(i,2))
  + kt(k,i+1,0)*((t-pt(i,0))+(t-pt(i,1)))/(pt(i,2)-pt(i,0))/(pt(i,2)-pt(i,1)) ];

#CFR reaction rates
  rate1(k E K): rt(k,1) =e= k0(1)*exp[-e(1)/rg/tt(k)]*ct(k,1)*ct(k,2)^0.3;
  rate2(k E K): rt(k,2) =e= k0(2)*exp[-e(2)/rg/tt(k)]*ct(k,1)^0.5*ct(k,2)^1.8;

#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  )*(t-pt(i,2))/(pt(i,0)-pt(i,2))
                             + ks(k,i+1)*(t-pt(i,0))/(pt(i,2)-pt(i,0)) ];
  dcontrol1(k E K): dsdv(k) =e= interv[i E CI|
                                 ks(k,i  )/(pt(i,0)-pt(i,2))
                               + ks(k,i+1)/(pt(i,2)-pt(i,0)) ];

#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  )*(t-pt(i,2))/(pt(i,0)-pt(i,2))
                             + ko(k,i+1)*(t-pt(i,0))/(pt(i,2)-pt(i,0)) ];
  dcontrol2(k E K): dodv(k) =e= interv[i E CI|
                                 ko(k,i  )/(pt(i,0)-pt(i,2))
                               + ko(k,i+1)/(pt(i,2)-pt(i,0)) ];

#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)>>;

#continuity of derivatives for temperature control
  cont1(k E K, i E CI & i>0):
                  kt(k,i-1,0)*(pt(i-1,2)-pt(i-1,1))/
                                 (pt(i-1,0)-pt(i-1,1))/(pt(i-1,0)-pt(i-1,2))
                + kt(k,i-1,1)*(pt(i-1,2)-pt(i-1,0))/
                                 (pt(i-1,1)-pt(i-1,0))/(pt(i-1,1)-pt(i-1,2))
                + kt(k,i  ,0)*(2*pt(i-1,2)-pt(i-1,1)-pt(i-1,0))/
                                 (pt(i-1,2)-pt(i-1,0))/(pt(i-1,2)-pt(i-1,1))
                - kt(k,i  ,0)*(2*pt(i,0)-pt(i,2)-pt(i,1))/
                                 (pt(i,0)-pt(i,1))/(pt(i,0)-pt(i,2))
                - kt(k,i  ,1)*(pt(i,0)-pt(i,2))/
                                 (pt(i,1)-pt(i,0))/(pt(i,1)-pt(i,2))
                - kt(k,i+1,0)*(pt(i,0)-pt(i,1))/
                                 (pt(i,2)-pt(i,0))/(pt(i,2)-pt(i,1))
          =e= 0;

#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);
}}