$title Synthesis: Hydrodealkylation of Toluene (HDA,SEQ=123)

$onText
The chemical plant performed the hydro-dealkylation of toluene
into benzene and methane. The flowsheet model was used to make
decisions on choosing between alternative process units at
various stages of the process. The resulting model is nonlinear
and mixed integer.


Kocis, G R, and Grossmann, I E, Computational Experience with DICOPT
Solving Minlp Problems in Process Synthesis. Computers and Chemical
Engineering 13, 3 (1989), 307-315.


              -----------------
              - process units -
              -----------------

   absorber                           abs
   compressors                        comp
   distillation                       dist
   flash                              flsh
   furnace                            furn
   utility-stream exchangers          hec & heh
   stream-stream exchanger            exch
   membrane separator                 memb
   mixers (single inlet stream)       mxr1
   mixers                             mxr
   pump                               pump
   reactor                            rct
   splitters (single outlet stream)   spl1
   splitters                          spl
   expansion valve                    valve


              -----------------------
              - process unit models -
              -----------------------

   absorber                         kremser equation
                                    (douglas, page 431)
                                       neglect heat effects
                                       assume:  pure solvent
                                                33 % tray efficiency
                                                l / m g = 1.4 (benzene)
                                                fixed recoveries for
                                                      h2 , ch4 , dip

   compressors                      isentropic compression of ideal gas
                                    (mccabe and smith, page 201)

   distillation                     rmin from underwood equation
                                    (douglas, page 441)

                                    r = 1.2 * rmin  (heuristic)

                                    nmin from fenske equation
                                    (douglas, page 441)

                                    nt = 2 * nmin (gillilands approx.)
                                    (douglas, page 443)

                                    n = nt / efficiency

   expansion valve                  isentropic expansion of ideal gas
                                    (mccabe and smith, page 200)

   flash                            ideal flash (raoult's law)
                                    (mccabe and smith, page 484)

   furnace                          50 % efficiency

   heat exchangers

   membrane separator               shortcut method (driving force
                                    approximated as arithmetic mean)
                                    (grossmann, design notes)

   mixers (single inlet stream)     linear model for heat & mass balances
                                    (kocis and grossmann, 1988)

   mixers                           nonlinear heat balance
                                    (kocis and grossmann, 1988)

   pump                             pout > pin

   reactor                          plug flow reactor
                                    1.  isothermal
                                    2.  adiabatic - use average temp

   splitters (single outlet stream) linear model for heat & mass balances
                                    (kocis and grossmann, 1988)

   splitters                        split fraction model
                                    (kocis and grossmann, 1988)

Keywords: mixed integer nonlinear programming, chemical engineering, toluene
          hydrodealkylation, process synthesis
$offText

Scalar
   alpha    'compressor coefficient'             / 0.3665 /
   compeff  'compressor efficiency'              / 0.750  /
   gam      'ratio of cp to cv'                  / 1.300  /
   abseff   'absorber tray efficiency'           / 0.333  /
   disteff  'column tray efficiency'             / 0.500  /
   uflow    'upper bound - flow logicals'        / 50     /
   upress   'upper bound - pressure logicals'    / 4.0    /
   utemp    'upper bound - temperature logicals' / 7.0    /

*  want costs in units of  $1000. per year
   costelec 'electricity cost ($340. per kw-yr)' / 0.340 /
   costqc   'cooling cost    ($700  per 1e9 kj)' / 0.700 /
   costqh   'heating cost    ($8000 per 1e9 kj)' / 8.000 /
   costfuel 'fuel cost furnace ($4 per 1e6 btu)' / 4.0   /;

$onText
user must fill in the data for the following sets to define
the number of process streams, components, and process units
$offText

Set
   str    'process streams'     / 1*35, 37*73 /
   compon 'chemical components' / h2 'hydrogen', ch4 'methane', ben 'benzene', tol 'toluene', dip 'diphenyl' /
   abs    'absorber'                   / 1   /
   comp   'compressors'                / 1*4 /
   dist   'distillation'               / 1*3 /
   flsh   'flash separators'           / 1*3 /
   furn   'furnace'                    / 1   /
   hec    'cooling units'              / 1*2 /
   heh    'heating units'              / 1*4 /
   exch   'exchanger'                  / 1   /
   memb   'membrane separators'        / 1*2 /
   mxr1   'mixers (1 input stream)'    / 1*5 /
   mxr    'mixers'                     / 1*5 /
   pump   'pump'                       / 1*2 /
   rct    'reactors'                   / 1*2 /
   spl1   'splitters (1 input stream)' / 1*6 /
   spl    'splitters'                  / 1*3 /
   valve  'expansion valves'           / 1*6 /

$onText
user must fill in the data for the following sets to
define the input-output structure of the superstructure
by defining which streams are input and/or output streams
for the process units
$offText

* =======================
* =                     =
* =   a b s o r b e r   =
* =                     =
* =======================
   ilabs(abs,str)    'abs-stream (inlet liquid) matches'  / 1.67           /
   olabs(abs,str)    'abs-stream (outlet liquid) matches' / 1.68           /
   ivabs(abs,str)    'abs-stream (inlet vapor) matches'   / 1.63           /
   ovabs(abs,str)    'abs-stream (outlet vapor) matches'  / 1.64           /
   asolv(abs,compon) 'abs-solvent component matches'      / 1.tol          /
   anorm(abs,compon) 'abs-comp matches (normal model)'    / 1.ben          /
   asimp(abs,compon) 'abs-heavy component matches'        / 1.(h2,ch4,dip) /

* ===========================
* =                         =
* =   c o m p r e s s o r   =
* =                         =
* ===========================
   icomp(comp,str) 'compressor-stream (inlet) matches'  / 1.5, 2.59, 3.64, 4.56 /
   ocomp(comp,str) 'compressor-stream (outlet) matches' / 1.6, 2.60, 3.65, 4.57 /

* ===============================
* =                             =
* =   d i s t i l l a t i o n   =
* =                             =
* ===============================
   idist(dist,str)    'dist-stream (inlet) matches'      / 1.25, 2.30, 3.33          /
   vdist(dist,str)    'dist-stream (vapor) matches'      / 1.26, 2.31, 3.34          /
   ldist(dist,str)    'dist-stream (liquid) matches'     / 1.27, 2.32, 3.35          /
   dl(dist,compon)    'dist-light components matches'    / 1.h2       , 2.ch4, 3.ben /
   dlkey(dist,compon) 'dist-light key component matches' / 1.ch4      , 2.ben, 3.tol /
   dhkey(dist,compon) 'dist-heavy key component matches' / 1.ben      , 2.tol, 3.dip /
   dh(dist,compon)    'dist-heavy components matches'    / 1.(tol,dip), 2.dip        /
   dkey(dist,compon)  'dist-key component matches';

dkey(dist,compon) = dlkey(dist,compon) + dhkey(dist,compon);

* =================
* =               =
* =   f l a s h   =
* =               =
* =================
Set
   iflsh(flsh,str)   'flsh-stream (inlet) matches'  / 1.17 , 2.46 , 3.39  /
   vflsh(flsh,str)   'flsh-stream (vapor) matches'  / 1.18 , 2.47 , 3.40  /
   lflsh(flsh,str)   'flsh-stream (liquid) matches' / 1.19 , 2.48 , 3.41  /
   fkey(flsh,compon) 'flash-key component matches'  / 1.ch4, 2.ch4, 3.tol /

* =======================
* =                     =
* =    f u r n a c e    =
* =                     =
* =======================
   ifurn(furn,str) 'furn-stream (inlet) matches'  / 1.70 /
   ofurn(furn,str) 'furn-stream (outlet) matches' / 1.9  /

* ===================
* =                 =
* =   c o o l e r   =
* =                 =
* ===================
   ihec(hec,str) 'hec-stream (inlet) matches'  / 1.71, 2.45 /
   ohec(hec,str) 'hec-stream (outlet) matches' / 1.17, 2.46 /

* ===================
* =                 =
* =   h e a t e r   =
* =                 =
* ===================
   iheh(heh,str) 'heh-stream (inlet) matches'  / 1.24, 2.23, 3.37, 4.61 /
   oheh(heh,str) 'heh-stream (outlet) matches' / 1.25, 2.44, 3.38, 4.73 /

* =========================
* =                       =
* =   e x c h a n g e r   =
* =                       =
* =========================
   icexch(exch,str) 'exch-cold stream (inlet)  matches' / 1.8  /
   ocexch(exch,str) 'exch-cold stream (outlet) matches' / 1.70 /
   ihexch(exch,str) 'exch-hot  stream (inlet)  matches' / 1.16 /
   ohexch(exch,str) 'exch-hot  stream (outlet) matches' / 1.71 /

* =======================
* =                     =
* =   m e m b r a n e   =
* =                     =
* =======================
   imemb(memb,str)    'memb-stream (inlet) matches'        / 1.3, 2.54 /
   nmemb(memb,str)    'memb-stream (non-permeate) matches' / 1.4, 2.55 /
   pmemb(memb,str)    'memb-stream (permeate) matches'     / 1.5, 2.56 /
   mnorm(memb,compon) 'normal components'                  / 1.(h2,ch4)     , 2.(h2,ch4)      /
   msimp(memb,compon) 'simplified flux components'         / 1.(ben,tol,dip), 2.(ben,tol,dip) /

* ==================
* =                =
* =   m i x e r    =
* = (single input) =
* =                =
* ==================
   imxr1(mxr1,str) 'mixer-stream (inlet) matches'  / 1.(2,6), 2.(11,13), 3.(27,48), 4.(34,40), 5.(49,50) /
   omxr1(mxr1,str) 'mixer-stream (outlet) matches' / 1.7    , 2.14     , 3.30     , 4.42     , 5.51      /
   mxr1spl1(mxr1,str,str) '1-mxr-inlet 1-spl-outlet matches'  / 1.(2.2,6.3)    , 2.(11.10,13.12)
                                                                3.(27.24,48.23), 4.(34.33,40.37)
                                                                5.(49.23,50.24)                          /

* =================
* =               =
* =   m i x e r   =
* =               =
* =================
   imxr(mxr,str) '-stream (inlet) matches'       / 1.(7,43,66,72), 2.(15,20)
                                                   3.(21,69)     , 4.(51,62)
                                                   5.(57,60,65)                /
   omxr(mxr,str) 'mixer-stream (outlet) matches' / 1.8, 2.16, 3.22, 4.63, 5.72 /

* ===============
* =             =
* =   p u m p   =
* =             =
* ===============
   ipump(pump,str) 'pump-stream (inlet) matches' / 1.42, 2.68 /
   opump(pump,str) 'pump-stream (outlet) matches'/ 1.43, 2.69 /

* =====================
* =                   =
* =   r e a c t o r   =
* =                   =
* =====================
   irct(rct,str)    'reactor-stream (inlet) matches'  / 1.10 , 2.12  /
   orct(rct,str)    'reactor-stream (outlet) matches' / 1.11 , 2.13  /
   rkey(rct,compon) 'reactor-key component matches'   / 1.tol, 2.tol /

* =======================
* =                     =
* =   s p l i t t e r   =
* =   (single output)   =
* =                     =
* =======================
   ispl1(spl1,str) 'splitter-stream (inlet) matches'  / 1.1      , 2.9      , 3.22
                                                        4.32     , 5.52     , 6.58      /
   ospl1(spl1,str) 'splitter-stream (outlet) matches' / 1.(2,3)  , 2.(10,12), 3.(23,24)
                                                        4.(33,37), 5.(53,54), 6.(59,61) /

* =======================
* =                     =
* =   s p l i t t e r   =
* =                     =
* =======================
   ispl(spl,str) 'splitter-stream (inlet) matches'  / 1.19, 2.18, 3.26 /
   ospl(spl,str) 'splitter-stream (outlet) matches' / 1.(20,21), 2.(52,58), 3.(28,29) /

* =========================
* =                       =
* =   e x p   v a l v e   =
* =                       =
* =========================
   ival(valve,str) 'exp.valve-stream (inlet) matches'  / 1.44, 2.38, 3.14
                                                         4.47, 5.29, 6.73 /
   oval(valve,str) 'exp.valve-stream (outlet) matches' / 1.45, 2.39, 3.15
                                                         4.49, 5.50, 6.62 /;

Alias (str,str2), (compon,compon2);

$onText
heat capacities for pure components  ( kj per kg-mol k ) -

              vapor  liquid
   h2            30
   ch4           40
   ben & tol    220     230
   dip          440     460

  (h2, ch4, and tol vapor values from douglas
   tol      liquid value from mccabe and smith
   ben      liquid and vapor values assumed the same as tol
   dip      liquid and vapor values assumed 2 times ben
            i.e. molecular weight of dip = 2 * ben )

approximate bulk stream cp by guessing a stream composition
and calculating a weighted average.

heat of vaporization of toluene taken from reid, prausnitz & sherwood

   30890 kj per kg-mole

(same value assumed for ben if needed)
$offText

Parameter
   cp(str)         'heat capacities ( kj per kgmole-k)'
   heatvap(compon) 'heat of vaporization (kj per kg-mol)'
   cppure(compon)  'pure component heat capacities' / h2 30, ch4 40, ben 225, tol 225, dip 450 /;

heatvap('tol') = 30890.;

Table gcomp(str,compon) 'guess composition values'
          h2   ch4   ben   tol  dip
    7   0.95  0.05
    8   0.50  0.40        0.10
    9   0.50  0.40        0.10
   10   0.50  0.40        0.10
   11   0.45  0.45  0.05  0.05
   12   0.50  0.40        0.10
   13   0.45  0.45  0.05  0.05
   14   0.45  0.45  0.05  0.05
   15   0.45  0.45  0.05  0.05
   16   0.40  0.40  0.10  0.10
   17   0.40  0.40  0.10  0.10
   20   0.03  0.07  0.55  0.35
   21   0.03  0.07  0.55  0.35
   22   0.03  0.07  0.55  0.35
   24   0.03  0.07  0.55  0.35
   25   0.03  0.07  0.55  0.35
   37                     1.00
   38                     1.00
   43               0.05  0.95
   44   0.03  0.07  0.55  0.35
   45   0.03  0.07  0.55  0.35
   46   0.03  0.07  0.55  0.35
   51   0.30  0.70
   57   0.80  0.20
   60   0.50  0.50
   62   0.50  0.50
*  Line for stream 63 estimated by Adrud based on Mixer data
   63   0.47  0.40  0.01  0.12
   65   0.50  0.50
   66                     1.00
   69                     1.00
   70   0.50  0.40        0.10
   71   0.40  0.40  0.10  0.10
   72   0.50  0.50            ;

cp(str) = sum(compon, cppure(compon)*gcomp(str,compon));

display cp;

Parameter
   anta(compon)  'antoine coefficient'
                 / h2 13.6333, ch4 15.2243, ben 15.9008, tol 16.0137, dip 16.6832 /
   antb(compon)  'antoine coefficient'
                 / h2 164.9, ch4 897.84, ben 2788.51, tol 3096.52, dip 4602.23    /
   antc(compon)  'antoine coefficient'
                 / h2 3.19, ch4 -7.16, ben -52.36, tol -53.67, dip -70.42 /
   perm(compon)  'permeability'
                 / h2 55.0e-06, ch4 2.3e-06 /
   cbeta(compon) 'constant values (exp(beta)) in absorber'
                 / h2 1.0003, ch4 1.0008, dip 1.0e+04 /
   aabs(compon)  'absorption factors'
                 / ben 1.4, tol 4.0 /;

$onText
aabs(ben)  fixed at 1.4 (hueristic) and
aabs(tol)  calculated as:  aabs(ben)*vp(tol)/vp(ben) at
           t = 3.0 (temperature of absorber)  --->  4.0

units: cc (stp) / cm**2 - sec - cm hg
       vol      / area - time - pressure

       1 g-mol = 22.4 l = 22.4 * 1000 cc

       1000 g-mol = 1 kg-mole

       1 m**2  = ( 100 cm )**2 = 1.0e+04 cm**2

       76 cm hg = 1 atm = 0.101325 mpa
       --->     750.062 cm hg = 1 mpa

       60 sec = 1 min

convert units to: kg-mol / m**2 - min - mpa

1 psia = 6.894757e+03 pa = 6.894757e-03 mpa
$offText

perm(compon) = perm(compon)*(1./22400.)*1.0e+04*750.062*60./1000.;

Scalar eps1 'small number to avoid div. by zero' / 1.e-04 /;

Parameter
   heatrxn(rct)    'heat of reaction   (kj per kg-mol)'
                   / 1 50100., 2 50100. /
   f1comp(compon)  'feedstock compositions   (h2 feed)'
                   / h2 0.95, ch4 0.05, ben 0.00, tol 0.00, dip 0.00 /
   f66comp(compon) 'feedstock compositions  (tol feed)'
                   / h2 0.00, ch4 0.00, ben 0.00, tol 1.00, dip 0.00 /
   f67comp(compon) 'feedstock compositions  (tol feed)'
                   / h2 0.00, ch4 0.00, ben 0.00, tol 1.00, dip 0.00 /;

$onText
reaction: tol + h2 ---> ben + ch4 (desired)
          2 ben    ---> dip + h2  (undesired)
$offText

* =======================
* =                     =
* =   a b s o r b e r   =
* =                     =
* =======================
Positive Variable nabs(abs) 'number of absorber trays';

Variable gamma(abs,compon), beta(abs,compon);

* ===========================
* =                         =
* =   c o m p r e s s o r   =
* =                         =
* ===========================
Positive Variable
   elec(comp)    'electricity requirement (kw)'
   presrat(comp) 'ratio of outlet to inlet pressure'

* ===============================
* =                             =
* =   d i s t i l l a t i o n   =
* =                             =
* ===============================
   nmin(dist)   'minimum number of trays in column'
   ndist(dist)  'number of trays in column'
   rmin(dist)   'minimum reflux ratio'
   reflux(dist) 'reflux ratio'
   distp(dist)  'column pressure'
   distt(dist)  'column temperature'
   avevlt(dist) 'average volatility'

* =================
* =               =
* =   f l a s h   =
* =               =
* =================
   flsht(flsh)        'flash temperature    (100 k)'
   flshp(flsh)        'flash pressure (mega-pascal)'
   eflsh(flsh,compon) 'vapor phase recovery in flash'

* =======================
* =                     =
* =    f u r n a c e    =
* =                     =
* =======================
   qfuel(furn) 'heating requied (1.e+12 kj per yr)'

* ===================
* =                 =
* =   c o o l e r   =
* =                 =
* ===================
   qc(hec) 'utility requirement (1.e+12 kj per yr)'

* ===================
* =                 =
* =   h e a t e r   =
* =                 =
* ===================
   qh(heh) 'utility requirement (1.e+12 kj per yr)'

* =========================
* =                       =
* =   e x c h a n g e r   =
* =                       =
* =========================
   qexch(exch) 'heat exchanged (1.e+12 kj per yr)'

* =======================
* =                     =
* =   m e m b r a n e   =
* =                     =
* =======================
   a(memb) 'surface area for mass transfer (m**2)'

* =================
* =               =
* =   m i x e r   =
* =   (1 input)   =
* =               =
* =================
   mxr1t(mxr1) 'mixer temperature (100 k)'
   mxr1p(mxr1) 'mixer pressure     (m-pa)'

* =================
* =               =
* =   m i x e r   =
* =               =
* =================
   mxrt(mxr) 'mixer temperature (100 k)'
   mxrp(mxr) 'mixer pressure     (m-pa)'

* ===============
* =             =
* =   p u m p   =
* =             =
* ===============

* =====================
* =                   =
* =   r e a c t o r   =
* =                   =
* =====================
   rctt(rct)          'reactor temperature    (100 k)'
   rctp(rct)          'reactor pressure        (m-pa)'
   rctvol(rct)        'reactor volume   (cubic meter)'
   krct(rct)          'rate constant'
   conv(rct,compon)   'conversion of key component'
   sel(rct)           'selectivity to Benzene'
   consum(rct,compon) 'consumption rate of key'
   q(rct)             'heat removed (1.e+9 kj per yr)'

* =======================
* =                     =
* =   s p l i t t e r   =
* =     (1 output)      =
* =                     =
* =======================
   spl1t(spl1) 'splitter temperature (100 k)'
   spl1p(spl1) 'splitter pressure     (m-pa)'

* =======================
* =                     =
* =   s p l i t t e r   =
* =                     =
* =======================
   splt(spl) 'splitter temperature (100 k)'
   splp(spl) 'splitter pressure     (m-pa)'

* =====================
* =                   =
* =   s t r e a m s   =
* =                   =
* =====================
   vp(str,compon) 'vapor pressure          (mega-pascal)'
   e(str)         'split fraction'
   f(str)         'stream flowrates    (kg-mole per min)'
   fc(str,compon) 'component flowrates (kg-mole per min)'
   p(str)         'stream pressure         (mega_pascal)'
   t(str)         'stream temperature            (100 k)';

Binary Variable y(str) 'binary variable';

Variable
   const  'constant term in obj fcn'
   profit 'overall profit ($1000 per year)';

* =======================
* =                     =
* =   a b s o r b e r   =
* =                     =
* =======================
Equation
   absfact(abs,compon)  'absorbption factor equation'
   gameqn(abs,compon)   'definition of gamma'
   betaeqn(abs,compon)  'definition of beta'
   abssvrec(abs,compon) 'recovery of solvent'
   absrec(abs,compon)   'recovery of non-solvent'
   abssimp(abs,compon)  'recovery of simplified components'
   abscmb(abs,compon)   'overall component mass balance'
   abslogc              'logical constraint'
   abspl(abs)           'pressure relation for liquid'
   abstl(abs)           'temperature relation for liquid'
   abspv(abs)           'pressure relation for vapor'
   absp(abs)            'pressure relation for inlets'
   abst(abs)            'temperature relation at top';

absfact(abs,compon)$anorm(abs,compon)..
   sum(str$ilabs(abs,str), f(str)*p(str))     =e=
   sum(str$ivabs(abs,str), f(str))*aabs(compon)*
   sum(str$ilabs(abs,str), vp.l(str,compon));

gameqn(abs,compon)$asolv(abs,compon)..
   gamma(abs,compon) =e= log((1 - aabs(compon)**(nabs(abs)*abseff + eps1))/(1 - aabs(compon)));                                ;

betaeqn(abs,compon)$(not asimp(abs,compon))..
   beta(abs,compon)  =e= log((1 - aabs(compon)**((nabs(abs)*abseff) + 1))/(1 - aabs(compon)));                                 ;

abssvrec(abs,compon)$asolv(abs,compon)..
   sum(str$ovabs(abs,str), fc(str,compon))*exp(beta(abs,compon)) =e=
   sum(str$ivabs(abs,str), fc(str,compon)) +
   exp(gamma(abs,compon))*sum(str$ilabs(abs,str), fc(str,compon));

absrec(abs,compon)$anorm(abs,compon)..
   sum(str$ovabs(abs,str), fc(str,compon))*exp(beta(abs,compon)) =e=
   sum(str$ivabs(abs,str), fc(str,compon));

abssimp(abs,compon)$asimp(abs,compon)..
   sum(str$ovabs(abs,str), fc(str,compon)) =e=
   sum(str$ivabs(abs,str), fc(str,compon))/cbeta(compon);

abscmb(abs,compon)..
   sum(str$ilabs(abs,str), fc(str,compon))  +
   sum(str$ivabs(abs,str), fc(str,compon)) =e=
   sum(str$olabs(abs,str), fc(str,compon))  +
   sum(str$ovabs(abs,str), fc(str,compon));

abslogc..
   f('63') + f('67') =l= 25.*y('63');

abspl(abs)..
   sum(str$ilabs(abs,str), p(str))  =e=
   sum(str$olabs(abs,str), p(str));

abstl(abs)..
   sum(str$ilabs(abs,str), t(str)) =e= sum(str$olabs(abs,str), t(str));

abspv(abs)..
   sum(str$ivabs(abs,str), p(str)) =e= sum(str$ovabs(abs,str), p(str));

absp(abs)..
   sum(str$ilabs(abs,str), p(str)) =e= sum(str$ivabs(abs,str), p(str));

abst(abs)..
   sum(str$ilabs(abs,str), t(str)) =e= sum(str$ovabs(abs,str), t(str));

* ===========================
* =                         =
* =   c o m p r e s s o r   =
* =                         =
* ===========================
Equation
   compcmb(comp,compon) 'component balance for compressor'
   comphb(comp)         'heat balance for compressor'
   compelec(comp)       'energy balance for compressor'
   ratio(comp)          'pressure ratio (out to in)';

compcmb(comp,compon)..
   sum(str$ocomp(comp,str), fc(str,compon)) =e= sum(str$icomp(comp,str), fc(str,compon));

comphb(comp)..
   sum(str$ocomp(comp,str), t(str)) =e= presrat(comp)*sum(str$icomp(comp,str), t(str));

compelec(comp)..
   elec(comp) =e= alpha*(presrat(comp) - 1)
               *  sum(str$icomp(comp,str), 100.*t(str)*f(str)/60.)
               * (1/compeff)*(gam/(gam - 1.));

ratio(comp)..
   presrat(comp)**(gam/(gam - 1.)) =e= sum(str$ocomp(comp,str), p(str))
                                    /  sum(str$icomp(comp,str), p(str));

* ===============================
* =                             =
* =   d i s t i l l a t i o n   =
* =                             =
* ===============================
Equation
   distcmb(dist,compon)      'component mass balance'
   antdistb(dist,str,compon) 'vapor pressure correlation (bot)'
   antdistt(dist,str,compon) 'vapor pressure correlation (top)'
   relvol(dist)              'average relative volatilty'
   undwood(dist)             'minimum reflux ratio equation'
   actreflux(dist)           'actual reflux ratio'
   fenske(dist)              'minimum number of trays'
   acttray(dist)             'actual number of trays'
   distspec(dist,str,compon) 'recovery specification'
   distheav(dist,compon)     'heavy components'
   distlite(dist,compon)     'light components'
   distpi(dist,str)          'inlet pressure relation'
   distvpl(dist,str)         'bottom vapor pressure relation'
   distvpv(dist,str)         'top vapor pressure relation'
   distpl(dist,str)          'outlet pressure relation(liquid)'
   distpv(dist,str)          'outlet pressure relation(vapor)';

distcmb(dist,compon)..
   sum(str$idist(dist,str), fc(str,compon)) =e=
   sum(str$vdist(dist,str), fc(str,compon))  +
   sum(str$ldist(dist,str), fc(str,compon));

antdistb(dist,str,compon)$(ldist(dist,str)$dkey(dist,compon))..
   log(vp(str,compon)*7500.6168) =e= anta(compon) - antb(compon)/(t(str)*100. + antc(compon));

antdistt(dist,str,compon)$(vdist(dist,str)$dkey(dist,compon))..
   log(vp(str,compon)*7500.6168) =e= anta(compon) - antb(compon)/(t(str)*100. + antc(compon));

relvol(dist)..
   avevlt(dist) =e= sqrt(sum(str$   ldist(dist,str)   ,
                         sum(compon$dlkey(dist,compon), vp(str,compon))  /
                         sum(compon$dhkey(dist,compon), vp(str,compon))) *
                         sum(str$   vdist(dist,str)   ,
                         sum(compon$dlkey(dist,compon), vp(str,compon))  /
                         sum(compon$dhkey(dist,compon), vp(str,compon))));

undwood(dist)..
   sum(str$idist(dist,str), sum(compon$dlkey(dist,compon), fc(str,compon)))*
   rmin(dist)*(avevlt(dist) - 1) =e= sum(str$idist(dist,str), sum(compon$dlkey(dist,compon), f(str)));

actreflux(dist).. reflux(dist) =e= 1.2*rmin(dist);

fenske(dist)..
   nmin(dist)*log(avevlt(dist)) =e=
   log(sum(str$vdist(dist,str),
   sum(compon$dhkey(dist,compon), (f(str) + eps1)/(fc(str,compon) + eps1)))*
   sum(str$ldist(dist,str),
   sum(compon$dlkey(dist,compon), (f(str) + eps1)/(fc(str,compon) + eps1))));

acttray(dist)..
   ndist(dist) =e= nmin(dist)*2./disteff;

distspec(dist,str,compon)$(vdist(dist,str)$dhkey(dist,compon))..
   fc(str,compon) =l= 0.05*sum(str2$idist(dist,str2), fc(str2,compon));

distheav(dist,compon)$dh(dist,compon)..
   sum(str$idist(dist,str), fc(str,compon)) =e= sum(str$ldist(dist,str), fc(str,compon));

distlite(dist,compon)$dl(dist,compon)..
   sum(str$idist(dist,str), fc(str,compon)) =e= sum(str$vdist(dist,str), fc(str,compon));

distpi(dist,str)$idist(dist,str)..
   distp(dist) =l= p(str);

distvpl(dist,str)$ldist(dist,str)..
   distp(dist) =e= sum(compon$dhkey(dist,compon), vp(str,compon));

distvpv(dist,str)$((ord(dist) <> 1)$vdist(dist,str))..
   distp(dist) =e= sum(compon$dlkey(dist,compon), vp(str,compon));

distpl(dist,str)$ldist(dist,str)..
   distp(dist) =e= p(str);

distpv(dist,str) $ vdist(dist,str)..
   distp(dist) =e= p(str);

* =================
* =               =
* =   f l a s h   =
* =               =
* =================
Equation
   flshcmb(flsh,compon)     'component mass balance'
   flshpr(flsh,str)         'flash pressure relation'
   flshpi(flsh,str)         'inlet pressure relation'
   flshpl(flsh,str)         'outlet pressure relation(liquid)'
   flshpv(flsh,str)         'outlet pressure relation(vapor)'
   flshti(flsh,str)         'inlet temp. relation'
   flshtl(flsh,str)         'outlet temp. relation(liquid)'
   flshtv(flsh,str)         'outlet temp. relation(vapor)'
   flshrec(flsh,str,compon) 'vapor recovery relation'
   flsheql(flsh,compon)     'equilibrium relation'
   antflsh(flsh,str,compon) 'vapor pressure correlation';

flshcmb(flsh,compon)..
   sum(str$iflsh(flsh,str), fc(str,compon)) =e=
   sum(str$vflsh(flsh,str), fc(str,compon))  +
   sum(str$lflsh(flsh,str), fc(str,compon));

antflsh(flsh,str,compon)$lflsh(flsh,str)..
   log(vp(str,compon)*7500.6168) =e= anta(compon) - antb(compon)/(t(str)*100. + antc(compon));

flshrec(flsh,str,compon)$lflsh(flsh,str)..
   sum(compon2$fkey(flsh,compon2), eflsh(flsh,compon2))*
   (eflsh(flsh,compon)*sum(compon2$fkey(flsh,compon2), vp(str,compon2)) +
   (1. - eflsh(flsh,compon))*vp(str,compon))                           =e=
   sum(compon2$fkey(flsh,compon2), vp(str,compon2))*eflsh(flsh,compon);

flsheql(flsh,compon)..
   sum(str$vflsh(flsh,str), fc(str,compon)) =e=
   sum(str$iflsh(flsh,str), fc(str,compon))*eflsh(flsh,compon);

flshpr(flsh,str)$lflsh(flsh,str)..
   flshp(flsh)*f(str) =e= sum(compon, vp(str,compon)*fc(str,compon));

flshpi(flsh,str)$iflsh(flsh,str)..
   flshp(flsh) =e= p(str);

flshpl(flsh,str)$lflsh(flsh,str)..
   flshp(flsh) =e= p(str);

flshpv(flsh,str)$vflsh(flsh,str)..
   flshp(flsh) =e= p(str);

flshti(flsh,str)$iflsh(flsh,str)..
   flsht(flsh) =e= t(str);

flshtl(flsh,str)$lflsh(flsh,str)..
   flsht(flsh) =e= t(str);

flshtv(flsh,str)$vflsh(flsh,str)..
   flsht(flsh) =e= t(str);

* =======================
* =                     =
* =    f u r n a c e    =
* =                     =
* =======================
Equation
   furnhb(furn)         'heat balance'
   furncmb(furn,compon) 'component mass balance'
   furnp(furn)          'pressure relation';

furnhb(furn)..
   qfuel(furn) =e= (sum(str$ofurn(furn,str), cp(str)*f(str)*100.*t(str))
                   -sum(str$ifurn(furn,str), cp(str)*f(str)*100.*t(str)))
                   *3600.*8500.*1.0e-12/60.;

furncmb(furn,compon)..
   sum(str$ofurn(furn,str), fc(str,compon)) =e= sum(str$ifurn(furn,str), fc(str,compon));

furnp(furn)..
   sum(str$ofurn(furn,str), p(str)) =e= sum(str$ifurn(furn,str), p(str)) - 0.4826;

* ===================
* =                 =
* =   c o o l e r   =
* =                 =
* ===================
Equation
   heccmb(hec,compon) 'component balance in cooler'
   hechb(hec)         'heat balance for cooler'
   hecp(hec)          'no pressure drop thru cooler';

heccmb(hec,compon)..
   sum(str$ohec(hec,str), fc(str,compon)) =e= sum(str$ihec(hec,str), fc(str,compon));

hechb(hec)..
   qc(hec) =e= (sum(str$ihec(hec,str), cp(str)*f(str)*100.*t(str))
               -sum(str$ohec(hec,str), cp(str)*f(str)*100.*t(str)))
               *3600.*8500.*1.0e-12/60.;
hecp(hec)..
   sum(str$ihec(hec,str), p(str)) =e= sum(str$ohec(hec,str), p(str));

* ===================
* =                 =
* =   h e a t e r   =
* =                 =
* ===================
Equation
   hehcmb(heh,compon) 'component balance in heater'
   hehhb(heh)         'heat balance for heater'
   hehp(heh)          'no pressure drop thru heater';

hehcmb(heh,compon)..
   sum(str$oheh(heh,str), fc(str,compon)) =e= sum(str$iheh(heh,str), fc(str,compon));

hehhb(heh)..
    qh(heh) =e= (sum(str$oheh(heh,str), cp(str)*f(str)*100.*t(str))
                -sum(str$iheh(heh,str), cp(str)*f(str)*100.*t(str)))
                *3600.*8500.*1.0e-12/60.;

hehp(heh).. sum(str$iheh(heh,str), p(str)) =e= sum(str$oheh(heh,str), p(str));

* =========================
* =                       =
* =   e x c h a n g e r   =
* =                       =
* =========================
Equation
   exchcmbc(exch,compon) 'component balance (cold)'
   exchcmbh(exch,compon) 'component balance (hot)'
   exchhbc(exch)         'heat balance for cold stream'
   exchhbh(exch)         'heat balance for hot  stream'
   exchdtm1(exch)        'delta t min condition'
   exchdtm2(exch)        'delta t min condition'
   exchpc(exch)          'pressure relation  (cold)'
   exchph(exch)          'pressure relation  (hot)';

exchcmbc(exch,compon)..
   sum(str$ocexch(exch,str), fc(str,compon)) =e= sum(str$icexch(exch,str), fc(str,compon));

exchcmbh(exch,compon)..
   sum(str$ohexch(exch,str), fc(str,compon)) =e= sum(str$ihexch(exch,str), fc(str,compon));

exchhbc(exch)..
   qexch(exch) =e= (sum(str$ocexch(exch,str), cp(str)*f(str)*100.*t(str))
                   -sum(str$icexch(exch,str), cp(str)*f(str)*100.*t(str)))
                   *3600.*8500.*1.0e-12/60.;

exchhbh(exch)..
   qexch(exch) =e= (sum(str$ihexch(exch,str), cp(str)*f(str)*100.*t(str))
                   -sum(str$ohexch(exch,str), cp(str)*f(str)*100.*t(str)))
                   *3600.*8500.*1.0e-12/60.;

exchdtm1(exch)..
   sum(str$ohexch(exch,str), t(str)) =g= sum(str$icexch(exch,str), t(str)) + 0.25;

exchdtm2(exch)..
   sum(str$ocexch(exch,str), t(str)) =l= sum(str$ihexch(exch,str), t(str)) - 0.25;

exchpc(exch)..
   sum(str$ocexch(exch,str), p(str)) =e= sum(str$icexch(exch,str), p(str));

exchph(exch)..
   sum(str$ohexch(exch,str), p(str)) =e= sum(str$ihexch(exch,str), p(str));

* =======================
* =                     =
* =   m e m b r a n e   =
* =                     =
* =======================
Equation
   memcmb(memb,str,compon) 'component mass balance'
   flux(memb,str,compon)   'mass flux relation'
   simp(memb,str,compon)   'mass flux relation (simplified)'
   memtp(memb,str)         'temp relation for permeate'
   mempp(memb,str)         'pressure relation for permeate'
   memtn(memb,str)         'temp relation for non-permeate'
   mempn(memb,str)         'pressure relation for non-permeate'
   rec(memb,str)           'recovery spec'
   pure(memb,str)          'purity spec';

memcmb(memb,str,compon)$imemb(memb,str)..
   fc(str,compon) =e= sum(str2$pmemb(memb,str2), fc(str2,compon))
                   +  sum(str2$nmemb(memb,str2), fc(str2,compon));

flux(memb,str,compon)$(pmemb(memb,str)$mnorm(memb,compon))..
   fc(str,compon) =e= a(memb)*perm(compon)/2.0
                      *(sum(str2$imemb(memb,str2), p(str2))
                      *(sum(str2$imemb(memb,str2), (fc(str2,compon) + eps1)/(f(str2) + eps1))  +
                        sum(str2$nmemb(memb,str2), (fc(str2,compon) + eps1)/(f(str2) + eps1))) -
                        2.*p(str)*(fc(str,compon) + eps1)/(f(str) + eps1));

simp(memb,str,compon)$(pmemb(memb,str)$msimp(memb,compon))..
   fc(str,compon) =e= 0.0;

memtp(memb,str)$pmemb(memb,str)..
   t(str) =e= sum(str2$imemb(memb,str2), t(str2));

mempp(memb,str)$pmemb(memb,str)..
   p(str) =l= sum(str2$imemb(memb,str2), p(str2));

memtn(memb,str)$nmemb(memb,str)..
   t(str) =e= sum(str2$imemb(memb,str2), t(str2));

mempn(memb,str)$nmemb(memb,str)..
   p(str) =e= sum(str2$imemb(memb,str2), p(str2));

rec(memb,str)$pmemb(memb,str)..
   fc(str,'h2') =g= 0.50*sum(str2$imemb(memb,str2), fc(str2,'h2'));

pure(memb,str)$pmemb(memb,str)..
   fc(str,'h2') =g= 0.50*f(str);

* ==================
* =                =
* =   m i x e r    =
* = (single input) =
* =                =
* ==================
Equation
   mxr1cmb(mxr1,str,compon) 'component balance in mixer'
   mxr1hbl(mxr1,str)        'heat balance in mixer (ls than)'
   mxr1hbg(mxr1,str)        'heat balance in mixer (grt than)'
   mxr1pl(mxr1,str)         'pressure relation (less than)'
   mxr1pg(mxr1,str)         'pressure relation (greater than)';

mxr1cmb(mxr1,str,compon)$omxr1(mxr1,str)..
   fc(str,compon) =e= sum(str2$imxr1(mxr1,str2), fc(str2,compon));

mxr1hbl(mxr1,str)$imxr1(mxr1,str)..
   sum(str2$omxr1(mxr1,str2), t(str2)) =l= t(str) + utemp*sum(str2$mxr1spl1(mxr1,str,str2), (1 - y(str2)));

mxr1hbg(mxr1,str)$imxr1(mxr1,str)..
   sum(str2$omxr1(mxr1,str2), t(str2)) =g= t(str) - utemp*sum(str2$mxr1spl1(mxr1,str,str2), (1 - y(str2)));

mxr1pl(mxr1,str)$imxr1(mxr1,str)..
   sum(str2$omxr1(mxr1,str2), p(str2)) =l= p(str) + upress*sum(str2$mxr1spl1(mxr1,str,str2), (1 - y(str2)));

mxr1pg(mxr1,str)$imxr1(mxr1,str)..
   sum(str2$omxr1(mxr1,str2), p(str2)) =g= p(str) - upress*sum(str2$mxr1spl1(mxr1,str,str2), (1 - y(str2)));

* =================
* =               =
* =   m i x e r   =
* =               =
* =================
Equation
   mxrcmb(mxr,compon) 'component balance in mixer'
   mxrhb(mxr)         'heat balance in mixer'
   mxrhbq(mxr)        'heat balance in quench'
   mxrpi(mxr,str)     'inlet pressure relation'
   mxrpo(mxr,str)     'outlet pressure relation';

mxrcmb(mxr,compon)..
   sum(str$omxr(mxr,str), fc(str,compon)) =e= sum(str$imxr(mxr,str), fc(str,compon));

mxrhb(mxr)$(ord(mxr) <> 2)..
   sum(str$imxr(mxr,str), f(str)*t(str)*cp(str)) =e= sum(str$omxr(mxr,str), f(str)*t(str)*cp(str));

mxrhbq(mxr)$(ord(mxr) = 2)..
   f('16')*t('16') =e= f('15')*t('15') - (fc('20','ben') + fc('20','tol'))*heatvap('tol')/(100.* cp('15'));

mxrpi(mxr,str)$imxr(mxr,str)..
   mxrp(mxr) =e= p(str);

mxrpo(mxr,str)$omxr(mxr,str)..
   mxrp(mxr) =e= p(str);

* ===============
* =             =
* =   p u m p   =
* =             =
* ===============
Equation
   pumpcmb(pump,compon) 'component balance'
   pumphb(pump)         'heat balance'
   pumppr(pump)         'pressure relation';

pumpcmb(pump,compon)..
   sum(str$opump(pump,str), fc(str,compon)) =e= sum(str$ipump(pump,str), fc(str,compon));

pumphb(pump)..
   sum(str$opump(pump,str), t(str)) =e= sum(str$ipump(pump,str), t(str));

pumppr(pump)..
   sum(str$opump(pump,str), p(str)) =g= sum(str$ipump(pump,str), p(str));

* =====================
* =                   =
* =   r e a c t o r   =
* =                   =
* =====================
Equation
   rctspec(rct,str)        'spec. on reactor feed stream'
   rxnrate(rct)            'reaction rate constant'
   rctconv(rct,str,compon) 'conversion of key component'
   rctsel(rct)             'selectivity to benzene'
   rctcns(rct,str,compon)  'consumption rate of key comp.'
   rctmbtol(rct)           'mass balance in reactor (tol)'
   rctmbben(rct)           'mass balance in reactor (ben)'
   rctmbdip(rct)           'mass balance in reactor (dip)'
   rctmbh2(rct)            'mass balance in reactor (h2)'
   rctmbch4(rct)           'mass balance in reactor (ch4)'
   rcthbadb(rct)           'heat balance (adiabatic)'
   rcthbiso(rct)           'heat balance (isothermal)'
   rctisot(rct)            'temp relation (isothermal)'
   rctpi(rct,str)          'inlet pressure relation'
   rctpo(rct,str)          'outlet pressure relation'
   rcttave(rct)            'average temperature relation';

rctspec(rct,str)$irct(rct,str)..
   fc(str,'h2') =g= 5.*(fc(str,'ben') + fc(str,'tol') + fc(str,'dip'));

rxnrate(rct)..
   krct(rct) =e= 6.3e+10*exp(-26167./(rctt(rct)*100.));

rctconv(rct,str,compon)$(rkey(rct,compon)$irct(rct,str))..
   1. - conv(rct,compon) =e= sqr(1./(1. + 0.372*krct(rct)*rctvol(rct)*sqrt(fc(str,compon)
                                   /60 + eps1)*(f(str)/60. + eps1)**(-3./2.)));

rctsel(rct)..
   (1. - sel(rct)) =e= 0.0036*(1. - conv(rct,'tol'))**(-1.544);

rctcns(rct,str,compon)$(rkey(rct,compon)$irct(rct,str))..
   consum(rct,compon) =e= conv(rct,compon)*fc(str,compon);

rctmbtol(rct)..
   sum(str$orct(rct,str), fc(str,'tol')) =e= sum(str$irct(rct,str), fc(str,'tol')) - consum(rct,'tol');

rctmbben(rct)..
   sum(str$orct(rct,str), fc(str,'ben')) =e= sum(str$irct(rct,str), fc(str,'ben')) + consum(rct,'tol')*sel(rct);

rctmbdip(rct)..
   sum(str$orct(rct,str), fc(str,'dip')) =e=
   sum(str$irct(rct,str), fc(str,'dip'))  +
   consum( rct,'tol')*0.5                 +
  (sum(str$irct(rct,str), fc(str,'ben'))  -
   sum(str$orct(rct,str), fc(str,'ben')))*0.5;

rctmbh2(rct)..
   sum(str$orct(rct,str), fc(str,'h2')) =e=
   sum(str$irct(rct,str), fc(str,'h2'))  -
   consum(rct,'tol')                     +
   sum(str$orct(rct,str), fc(str,'dip')) -
   sum(str$irct(rct,str), fc(str,'dip'));

rctmbch4(rct)..
   sum(str$orct(rct,str), fc(str,'ch4')) =e=
   sum(str$irct(rct,str), fc(str,'ch4'))  +
   consum(rct,'tol');

rcthbadb(rct)$(ord(rct) = 1)..
   heatrxn(rct)*consum(rct,'tol')/100.          =e=
   sum(str$orct(rct,str), cp(str)*f(str)*t(str)) -
   sum(str$irct(rct,str), cp(str)*f(str)*t(str));

rcthbiso(rct)$(ord(rct) = 2)..
   heatrxn(rct)*consum(rct,'tol')*60.*8500*1.0e-09 =e= q(rct);

rctisot(rct)$(ord(rct) = 2)..
   sum(str$orct(rct,str), t(str)) =e= sum(str$irct(rct,str), t(str));

rctpi(rct,str)$irct(rct,str)..
   rctp(rct) =e= p(str);

rctpo(rct,str)$orct(rct,str)..
   p(str) =e= rctp(rct) - 0.20684;

rcttave(rct)..
   rctt(rct) =e= (sum(str$irct(rct,str), t(str)) +
                  sum(str$orct(rct,str), t(str)))/2.;

* =======================
* =                     =
* =   s p l i t t e r   =
* =   (single output)   =
* =                     =
* =======================
Equation
   spl1logc(spl1,str)   'splitter logical constraint'
   spl1ysum(spl1)       'single choice'
   spl1cmb(spl1,compon) 'component balance in splitter'
   spl1pi(spl1,str)     'inlet pressure relation'
   spl1po(spl1,str)     'outlet pressure relation'
   spl1ti(spl1,str)     'inlet temperature relation'
   spl1to(spl1,str)     'outlet temperature relation';

spl1logc(spl1,str)$ospl1(spl1,str)..
   f(str) =l= uflow*y(str);

spl1ysum(spl1)..
   sum(str$ospl1(spl1,str), y(str)) =e= 1;

spl1cmb(spl1,compon)..
   sum(str$ospl1(spl1,str), fc(str,compon)) =e= sum(str$ispl1(spl1,str), fc(str,compon));

spl1pi(spl1,str)$ispl1(spl1,str)..
   spl1p(spl1) =e= p(str);

spl1po(spl1,str)$ospl1(spl1,str)..
   spl1p(spl1) =e= p(str);

spl1ti(spl1,str)$ispl1(spl1,str)..
   spl1t(spl1) =e= t(str);

spl1to(spl1,str)$ospl1(spl1,str)..
   spl1t(spl1) =e= t(str);

* =======================
* =                     =
* =   s p l i t t e r   =
* =                     =
* =======================
Equation
   splcmb(spl,str,compon) 'component balance in splitter'
   esum(spl)              'split fraction relation'
   splpi(spl,str)         'inlet pressure relation'
   splpo(spl,str)         'outlet pressure relation'
   splti(spl,str)         'inlet temperature relation'
   splto(spl,str)         'outlet temperature relation';

splcmb(spl,str,compon)$ospl(spl,str)..
   fc(str,compon) =e= sum(str2$ispl(spl,str2), fc(str2,compon)*e(str));

esum(spl)..
   sum(str$ospl(spl,str), e(str)) =e= 1.0;

splpi(spl,str)$ispl(spl,str)..
   splp(spl) =e= p(str);

splpo(spl,str)$ospl(spl,str)..
   splp(spl) =e= p(str);

splti(spl,str)$ispl(spl,str)..
   splt(spl) =e= t(str);

splto(spl,str)$ospl(spl,str)..
   splt(spl) =e= t(str);

* ==========================
* =                        =
* =   e x p    v a l v e   =
* =                        =
* ==========================
Equation
   valcmb(valve,compon) 'component mass balance'
   valt(valve)          'temperature relation'
   valp(valve)          'pressure relation';

valcmb(valve,compon)..
   sum(str$oval(valve,str), fc(str,compon)) =e= sum(str$ival(valve,str), fc(str,compon));

valt(valve)..
   sum(str$oval(valve,str), t(str)/(p(str)**((gam - 1.)/gam))) =e=
   sum(str$ival(valve,str), t(str)/(p(str)**((gam - 1.)/gam)));

valp(valve)..
   sum(str$oval(valve,str), p(str)) =l= sum(str$ival(valve,str), p(str));

* ===================
* =                 =
* =   s t r e a m   =
* =                 =
* ===================
Equation
   h2feed(compon)  'define h2  feed component flows'
   tolfeed(compon) 'define tol feed component flows'
   tolabs(compon)  'define tol absorber compon flows'
   specrec         'spec on h2 recycle'
   specprod        'spec on ben product'
   fbal(str)       'stream flow = sum of comp. flows';

specrec..
   fc('72','h2')  =g= 0.5*f('72');

specprod..
   fc('31','ben') =g= 0.9997*f('31');

fbal(str)..
   f(str) =e= sum(compon, fc(str,compon));

h2feed(compon)..
   fc('1',compon)  =e=  f('1')*f1comp(compon);

tolfeed(compon)..
   fc('66',compon) =e= f('66')*f66comp(compon);

tolabs(compon)..
   fc('67',compon) =e= f('67')*f67comp(compon);

* =========================
* =                       =
* =   o b j e c t i v e   =
* =                       =
* =========================
Equation
   abound 'fake bound on a(1)'
   obj    'cost function definition';

abound..
   a('1') =g= 0.0;

obj..
   profit =e= 510.*(- 2.5 *f('1')  - 14.  *(f('66') + f('67'))
                    + 19.9*f('31') + 11.84* f('35')
                    + 1.08*(fc('4','h2')   + fc('28','h2')  +
                            fc('53','h2')  + fc('55','h2'))
                    + 3.37*(fc('4','ch4')  + fc('28','ch4') +
                            fc('53','ch4') + fc('55','ch4')))
           -  0.815*(elec('1') + elec('2') + elec('3'))
           -  0.887* elec('4')
           -  7.155*(y('3') + y('59') + y('63'))
           -  4.866* y('54')
           -  sum(comp, costelec*elec(comp))
           - (74.3*y('10') + 1.257*rctvol('1') )
           -  1.25*(74.3*y('12') + 1.257*rctvol('2')) - 0.7*q('2')
           - (1.126*y('24') + 0.375*ndist('1')) - (+ 1.55*ndist('2'))
           - (3.9*y('33')   + 1.12*ndist('3'))
           - (43.24*y('3')  + 49.0*f('3'))
           - (43.24*y('54') + 49.0*f('54'))
           - (3.0  *y('63') + 1.2 *nabs('1'))
           - (    + 1171.7 * qfuel('1') + 4000.*qfuel('1'))
           -  sum(hec, 700 * qc(hec))
           -  sum(heh, 8000* qh(heh))
           -  const;


* + + + + + + + + + + + + + + + +
* +                             +
* +       b  o  u  n  d  s      +
* +                             +
* + + + + + + + + + + + + + + + +
profit.up = 1.0e+08;
const.fx  = 22.5;

* ===================
* =   s t r e a m   =
* ===================
y.lo(str) = 0.0;
y.up(str) = 1.0;

e.up(str)    = 1.0;
e.up('20')   = 0.50;
e.lo('21')   = 0.50;
* e.up('52') = 0.50;
* e.lo('58') = 0.50;

f.lo(str)  = 0.;
f.up(str)  = 10.;
f.up('8')  = 50.;
f.up('9')  = 50.;
f.up('10') = 50.;
f.up('11') = 50.;
f.up('12') = 50.;
f.up('13') = 50.;
f.up('14') = 50.;
f.up('15') = 50.;
f.up('16') = 50.;
f.up('17') = 50.;
f.up('18') = 50.;
f.up('52') = 50.;
f.up('54') = 50.;
f.up('56') = 50.;
f.up('57') = 50.;
f.up('58') = 50.;
f.up('59') = 50.;
f.up('60') = 50.;
f.up('70') = 50.;
f.up('71') = 50.;
f.up('72') = 50.;

fc.lo(str,compon) = 0.;
fc.up(str,compon) = 10.;

fc.up('8',compon)  = 30.;
fc.up('9',compon)  = 30.;
fc.up('10',compon) = 30.;
fc.up('11',compon) = 30.;
fc.up('12',compon) = 30.;
fc.up('13',compon) = 30.;
fc.up('14',compon) = 30.;
fc.up('15',compon) = 30.;
fc.up('16',compon) = 30.;
fc.up('17',compon) = 30.;
fc.up('18',compon) = 30.;
fc.up('52',compon) = 30.;
fc.up('54',compon) = 30.;
fc.up('56',compon) = 30.;
fc.up('57',compon) = 30.;
fc.up('58',compon) = 30.;
fc.up('59',compon) = 30.;
fc.up('60',compon) = 30.;
fc.up('70',compon) = 30.;
fc.up('71',compon) = 30.;
fc.up('72',compon) = 30.;

p.lo(str)  = 0.1;
p.up(str)  = 4.0;
t.lo(str)  = 3.00;
t.up(str)  = 10.00;
t.lo('49') = 2.0;
t.lo('50') = 2.0;
t.lo('51') = 2.0;

vp.lo(str,compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo(str)*100. + antc(compon)));
vp.up(str,compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up(str)*100. + antc(compon)));

* =======================
* =   a b s o r b e r   =
* =======================
nabs.lo(abs) = 0.;
nabs.up(abs) = 40.;

gamma.lo(abs,'tol') = log((1 - aabs('tol')**( nabs.lo('1')*abseff + eps1))/(1 - aabs('tol')));
gamma.up(abs,'tol') = min(15, log((1 - aabs('tol')**( nabs.up('1')*abseff + eps1))/(1 - aabs('tol'))));
beta.lo(abs,compon) = log((1 - aabs(compon)**( nabs.lo('1')*abseff + eps1 + 1.))/(1 - aabs(compon)));
beta.up(abs,compon) = min(15, log((1 - aabs(compon)**( nabs.up('1')*abseff + eps1 + 1.))/(1 - aabs(compon))));

t.fx('67') = 3.0;

vp.fx('67',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.l('67')*100. + antc(compon)));

* ===========================
* =   c o m p r e s s o r   =
* ===========================
presrat.lo(comp) = 1.;
presrat.up(comp) = 8./3.;

elec.lo(comp)= 0.;
elec.up(comp)= 100.;

* ===============================
* =   d i s t i l l a t i o n   =
* ===============================
nmin.lo('1') = 0.;
nmin.up('1') = 4.;
nmin.lo('2') = 8.;
nmin.up('2') = 14.;
nmin.lo('3') = 0.;
nmin.up('3') = 4.;

ndist.lo(dist) = nmin.lo(dist)*2./disteff;
ndist.up(dist) = nmin.up(dist)*2./disteff;

rmin.lo('1') = 0.02;
rmin.up('1') = 0.10;
rmin.lo('2') = 0.50;
rmin.up('2') = 2.00;
rmin.lo('3') = 0.02;
rmin.up('3') = 0.10;

reflux.lo(dist) = rmin.lo(dist)*1.2;
reflux.up(dist) = rmin.up(dist)*1.2;
distp.lo(dist)  = 0.1;
distp.up(dist)  = 4.0;

* fix pressure at 1.02 mpa (150 psia)
distp.fx('1') = 1.020;
distp.up('2') = 0.400;
distp.up('3') = 0.250;

* fix top temperature at 320 k
t.fx('26') = 3.20;

t.lo('27') = (antb('ben')/(anta('ben') - log(distp.lo('1')*7500.6168)) - antc('ben'))/100.;
t.up('27') = (antb('ben')/(anta('ben') - log(distp.up('1')*7500.6168)) - antc('ben'))/100.;
t.lo('31') = (antb('ben')/(anta('ben') - log(distp.lo('2')*7500.6168)) - antc('ben'))/100.;
t.up('31') = (antb('ben')/(anta('ben') - log(distp.up('2')*7500.6168)) - antc('ben'))/100.;
t.lo('32') = (antb('tol')/(anta('tol') - log(distp.lo('2')*7500.6168)) - antc('tol'))/100.;
t.up('32') = (antb('tol')/(anta('tol') - log(distp.up('2')*7500.6168)) - antc('tol'))/100.;
t.lo('34') = (antb('tol')/(anta('tol') - log(distp.lo('3')*7500.6168)) - antc('tol'))/100.;
t.up('34') = (antb('tol')/(anta('tol') - log(distp.up('3')*7500.6168)) - antc('tol'))/100.;
t.lo('35') = (antb('dip')/(anta('dip') - log(distp.lo('3')*7500.6168)) - antc('dip'))/100.;
t.up('35') = (antb('dip')/(anta('dip') - log(distp.up('3')*7500.6168)) - antc('dip'))/100.;

vp.lo(str,compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo(str)*100. + antc(compon)));
vp.up(str,compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up(str)*100. + antc(compon)));

avevlt.lo(dist) = sum(compon$dlkey(dist,compon),
                  sum(str   $ldist(dist,str)   ,vp.up(str,compon)))/
                  sum(compon$dhkey(dist,compon),
                  sum(str   $ldist(dist,str)   ,vp.up(str,compon)));

avevlt.up(dist) = sum(compon$dlkey(dist,compon),
                  sum(str   $vdist(dist,str)   , vp.lo(str,compon)))/
                  sum(compon$dhkey(dist,compon),
                  sum(str   $vdist(dist,str)   , vp.lo(str,compon)));

* ==================
* =   f l a s h    =
* ==================
eflsh.up(flsh,compon) = 1.;
eflsh.lo(flsh,compon) = 0.;

eflsh.lo('1','h2')  = 0.975;
eflsh.lo('1','ch4') = 0.975;
eflsh.up('1','ben') = 0.075;
eflsh.up('1','tol') = 0.050;
eflsh.up('1','dip') = 0.001;

eflsh.lo('2','h2')  = 0.975;
eflsh.lo('2','ch4') = 0.975;
eflsh.up('2','ben') = 0.20;
eflsh.up('2','tol') = 0.10;
eflsh.up('2','dip') = 0.01;

eflsh.lo('3','h2')  = 0.99;
eflsh.lo('3','ch4') = 0.99;
eflsh.lo('3','ben') = 0.90;
eflsh.lo('3','tol') = 0.90;
eflsh.up('3','dip') = 0.40;

flshp.lo('1') = 2.5;
flshp.up('1') = 3.25;
flshp.lo('2') = 0.1;
flshp.up('2') = 3.25;
flshp.lo('3') = 0.1;
flshp.up('3') = 3.25;
flsht.lo('1') = 3.00;
flsht.up('1') = 4.50;
flsht.lo('2') = 3.00;
flsht.up('2') = 4.50;
flsht.lo('3') = 3.00;
flsht.up('3') = 6.50;

t.lo('19') = flsht.lo('1');
t.up('19') = flsht.up('1');

vp.lo('19',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo('19')*100. + antc(compon)));
vp.up('19',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up('19')*100. + antc(compon)));

t.lo('48') = flsht.lo('2');
t.up('48') = flsht.up('2');

vp.lo('48',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo('48')*100. + antc(compon)));
vp.up('48',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up('48')*100. + antc(compon)));

*t.lo('41') = flsht.lo('3');
t.lo('41')  = t.lo('32');
t.up('41')  = flsht.up('3');

vp.lo('41',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo('41')*100. + antc(compon)));
vp.up('41',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up('41')*100. + antc(compon)));

* ======================
* =   f u r n a c e    =
* ======================
qfuel.up(furn) = 10.;

* ====================
* =   c o o l e r    =
* ====================
qc.up(hec) = 10.;

* ====================
* =   h e a t e r    =
* ====================
qh.up(heh) = 10.;

* ==========================
* =   e x c h a n g e r    =
* ==========================
qexch.up(exch) = 10.;

* =======================
* =   m e m b r a n e   =
* =======================
a.lo(memb) = 100;
a.up(memb) = 10000;

* ==================
* =   m i x e r    =
* =  (1 output)    =
* ==================
mxr1t.lo(mxr1) = 3.0;
mxr1t.up(mxr1) = 10.0;
mxr1p.lo(mxr1) = 0.1;
mxr1p.up(mxr1) = 4.0;

* ==================
* =   m i x e r    =
* ==================
mxrt.lo(mxr) = 3.0;
mxrt.up(mxr) = 10.0;
mxrp.lo(mxr) = 0.1;
mxrp.up(mxr) = 4.0;

* =====================
* =   r e a c t o r   =
* =====================
rctt.lo(rct)   = 8.9427;
rctt.up(rct)   = 9.7760;
rctp.fx(rct)   = 3.4474;
rctvol.up(rct) = 200.;
krct.lo(rct)   = 6.3e+10*exp(-26167./(rctt.lo(rct)*100.));
krct.up(rct)   = 6.3e+10*exp(-26167./(rctt.up(rct)*100.));

conv.up(rct,'tol') = 0.973;
sel.up(rct) = 1.0 - 0.0036;
q.up('2')   = 100.;

* =======================
* =   s p l i t t e r   =
* =     (1 outlet)      =
* =======================
spl1t.lo(spl1) = 3.00;
spl1t.up(spl1) = 10.00;
spl1p.lo(spl1) = 0.1;
spl1p.up(spl1) = 4.0;

* =======================
* =   s p l i t t e r   =
* =======================
splt.lo(spl) = 3.00;
splt.up(spl) = 10.00;
splp.lo(spl) = 0.1;
splp.up(spl) = 4.0;

* =========================
* =   o  t  h  e  r  s    =
* =========================
f.fx('31') = 2.08;
t.lo('9')  = 8.943;
t.up('16') = 8.943;
t.fx('1')  = 3.0;
p.fx('1')  = 3.930;
p.fx('66') = 3.930;
t.fx('66') = 3.00;

* + + + + + + + + + + + + + + + + + +
* +                                 +
* +       i  n  i  t  i  a  l       +
* +                                 +
* +          p  o  i  n  t          +
* +                                 +
* + + + + + + + + + + + + + + + + + +

* ---------------------------------------
* _  b i n a r y     v a r i a b l e s  _
* _______________________________________
y.l(str)   = 0.0;
y.l('2')   = 1.0;
*y.l('10') = 1.0;
y.l('12')  = 1.0;
y.l('24')  = 1.0;
*y.l('33') = 1.0;
y.l('37')  = 1.0;
*y.l('53') = 1.0;
y.l('54')  = 1.0;
y.l('59')  = 1.0;
*y.fx(str) = y.l(str);

* -----------------------------------------------
* _  c o n t i n u o u s    v a r i a b l e s   _
* _______________________________________________
t.l(str) = 3.0;
p.l(str) = 3.0;
f.l('1') = 4.5 - 2.0*y.l('54');
fc.l('1',compon) = f.l('1')*f1comp(compon);

* ---------------
* - b y p a s s -
* ---------------
fc.l('2',compon) = fc.l('1',compon)*y.l('2');
t.l('2') = t.l('1');
p.l('2') = p.l('1');

* -----------------------
* - m e m b r a n e   1 -
* -----------------------
fc.l('3',compon) = fc.l('1',compon)*y.l('3');
fc.l('4',compon) = fc.l('3',compon);
fc.l('5',compon) = fc.l('4',compon);

t.l('3') = t.l('1');
t.l('4') = t.l('1');
t.l('5') = t.l('1');
p.l('3') = p.l('3');
p.l('4') = p.l('3');
p.l('5') = p.l('3')*0.25;

* ---------------------------
* - c o m p r e s s o r   1 -
* ---------------------------
fc.l('6',compon) =  fc.l('5',compon);
p.l('6')         =  p.l('1');
presrat.l('1')   = (p.l('6')/p.l('5'))**0.23;
t.l('6')         =  t.l('5')*presrat.l('1');

* ---------------------
* - 1 - m i x e r   1 -
* ---------------------
fc.l('7',compon) = fc.l('1',compon) + fc.l('6',compon);
t.l('7') = t.l('2')*y.l('2') + t.l('6')*y.l('3');
p.l('7') = p.l('1');

* -----------------------
* - h 2   r e c y c l e -
* -----------------------
fc.l('57','h2')  = 2.5*y.l('54');
fc.l('57','ch4') = 0.5*y.l('54');
fc.l('57','ben') = 0.0*y.l('54');
fc.l('57','tol') = 0.0*y.l('54');
fc.l('57','dip') = 0.0*y.l('54');
p.l('57')        = p.l('1');

fc.l('60','h2')  = 14.0*y.l('59') - 2.0*y.l('54');
fc.l('60','ch4') = 14.0*y.l('59') - 2.0*y.l('54');
fc.l('60','ben') = 0.0 *y.l('59') + 0.0*y.l('54');
fc.l('60','tol') = 0.0 *y.l('59') + 0.0*y.l('54');
fc.l('60','dip') = 0.0 *y.l('59') + 0.0*y.l('54');
p.l('60')        = p.l('1');

fc.l('65','h2')  = 1.0*y.l('63');
fc.l('65','ch4') = 2.5*y.l('63');
fc.l('65','ben') = 0.0*y.l('63');
fc.l('65','tol') = 0.5*y.l('63');
fc.l('65','dip') = 0.0*y.l('63');
p.l('66')        = p.l('1');

fc.l('72',compon) = fc.l('60',compon) + fc.l('57',compon) + fc.l('65',compon);
p.l('72')         = p.l('1');
fc.l('59',compon) = fc.l('60',compon);
fc.l('56',compon) = fc.l('57',compon);
fc.l('64',compon) = fc.l('65',compon);

* -------------------------
* - t o l   r e c y c l e -
* -------------------------
fc.l('43',compon) = 0.0;
fc.l('43','ben')  = 0.1;
fc.l('43','tol')  = 2.0;
p.l('43')         = p.l('1');
fc.l('42',compon) = fc.l('43',compon);

* -------------------
* - t o l   f e e d -
* -------------------
f.l('66') = 2.1;
fc.l('66',compon) = f.l('66')*f66comp(compon);
p.l('66') = p.l('1');

* -----------------
* - m i x e r   1 -
* -----------------
fc.l('8',compon) = fc.l('7',compon) + fc.l('43',compon) + fc.l('66',compon) + fc.l('72',compon);

t.l('8')  =  3.5;
p.l('8')  =  p.l('1');

* -----------------------
* - e x c h a n g e r   1
* -----------------------
fc.l('70',compon) = fc.l('8',compon);
t.l('70') = 8.693;
p.l('70') = p.l('1');

* ------------------
* - f u r n a c e  -
* ------------------
fc.l('9',compon) = fc.l('8',compon);
t.l('9') = 8.943;
p.l('9') = p.l('1') - 0.4826;

* -----------------------------
* - 1 - s p l i t t e r    2  -
* -----------------------------
fc.l('10',compon) = fc.l('9',compon)*y.l('10');
fc.l('12',compon) = fc.l('9',compon)*y.l('12');

t.l('10') = t.l('9');
t.l('12') = t.l('9');
p.l('10') = p.l('9');
p.l('12') = p.l('9');

* ---------------------
* -   r e a c t o r   -
* ---------------------
f.l(str)    = sum(compon, fc.l(str,compon));
rctt.l(rct) = 9.2;
krct.l(rct) = 6.3e+10*exp(-26167./(rctt.l(rct)*100.));

rctvol.l('1') = 100.*y.l('10');
rctvol.l('2') = 100.*y.l('12');

conv.l('1','tol') = 1. - sqr(1./(1. + 0.372*krct.l('1')*rctvol.l('1')*sqrt(fc.l('10','tol')/60. + eps1)*(f.l('10')/60. + eps1)**(-3./2.)));
conv.l('2','tol') = 1. - sqr(1./(1. + 0.372*krct.l('2')*rctvol.l('2')*sqrt(fc.l('12','tol')/60. + eps1)*(f.l('12')/60. + eps1)**(-3./2.)));

sel.l(rct) = 1. - 0.0036*(1. - conv.l(rct,'tol'))**(-1.544);

consum.l('1','tol') = conv.l('1','tol')*fc.l('10','tol');
consum.l('2','tol') = conv.l('2','tol')*fc.l('12','tol');

fc.l('11','tol') = fc.l('10','tol') - consum.l('1','tol');
fc.l('11','ben') = fc.l('10','ben') + consum.l('1','tol')* sel.l('1');
fc.l('11','dip') = fc.l('10','dip') + consum.l('1','tol')*(1 - sel.l('1'))/2.;
fc.l('11','h2')  = fc.l('10','h2')  - consum.l('1','tol')*(1 + sel.l('1'))/2.;
fc.l('11','ch4') = fc.l('10','ch4') + consum.l('1','tol');
fc.l('13','tol') = fc.l('12','tol') - consum.l('2','tol');
fc.l('13','ben') = fc.l('12','ben') + consum.l('2','tol')* sel.l('2');
fc.l('13','dip') = fc.l('12','dip') + consum.l('2','tol')*(1 - sel.l('2'))/2.;
fc.l('13','h2')  = fc.l('12','h2')  - consum.l('2','tol')*(1 + sel.l('2'))/2.;
fc.l('13','ch4') = fc.l('12','ch4') + consum.l('2','tol');

f.l(str)  =  sum(compon, fc.l(str,compon));
t.l('11') = (heatrxn('1')*consum.l('1','tol')/(100.*cp('10')) + f.l('10')*t.l('10'))/(f.l('11') + eps1);
t.l('13') = (heatrxn('2')*consum.l('2','tol')/(100.*cp('12')) + f.l('12')*t.l('12'))/(f.l('13') + eps1);
p.l('11') =  p.l('9') - 0.2068;
p.l('13') =  p.l('9') - 0.2068;

* -------------------------
* -   1 - m i x e r   2   -
* -------------------------
fc.l('14',compon) = fc.l('11',compon) + fc.l('13',compon);
t.l('14') = t.l('11')*y.l('10')  +  t.l('13')*y.l('12');
p.l('14') = p.l('11');

* ----------------
* -   e x p   3  -
* ----------------
fc.l('15',compon) = fc.l('14',compon);
t.l('15') = t.l('14');
p.l('15') = p.l('14');

* --------------------
* -   m i x e r   2  -
* --------------------
fc.l('16','h2')  = fc.l('15','h2');
fc.l('16','ch4') = fc.l('15','ch4');
fc.l('16','ben') = fc.l('15','ben')*1.2;
fc.l('16','tol') = fc.l('15','tol')*1.2;
fc.l('16','dip') = fc.l('15','dip')*1.2;

t.l('16') = 8.943;
p.l('16') = p.l('15');

* -----------------------
* - e x c h a n g e r   1
* -----------------------
fc.l('71',compon) = fc.l('16',compon);
t.l('71') = 5.0;
p.l('71') = p.l('16');

* ------------------
* - c o o l e r  1 -
* ------------------
fc.l('17',compon) = fc.l('71',compon);
t.l('17') = 3.0;
p.l('17') = p.l('71');

* ---------------------
* -   f l a s h   1   -
* ---------------------
eflsh.l('1','h2')  = 0.995;
eflsh.l('1','ch4') = 0.99;
eflsh.l('1','ben') = 0.04;
eflsh.l('1','tol') = 0.01;
eflsh.l('1','dip') = 0.0001;
fc.l('18',compon)  = fc.l('17',compon)*eflsh.l('1',compon);
fc.l('19',compon)  = fc.l('17',compon) - fc.l('18',compon);

t.l('17') = 3.0;
t.l('18') = 3.0;
t.l('19') = 3.0;
p.l('18') = p.l('17');
p.l('19') = p.l('17');

vp.l('19',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.l('19')*100. + antc(compon)));

* -------------------------------------
* -   q u e n c h    s p l i t t e r  -
* -------------------------------------
fc.l('20',compon) = fc.l('19',compon)*0.20;
fc.l('21',compon) = fc.l('19',compon)*0.80;

t.l('20') = t.l('19');
t.l('21') = t.l('19');
p.l('20') = p.l('19');
p.l('21') = p.l('19');

* --------------------
* -   m i x e r   3  -
* --------------------
fc.l('69',compon) = 0.0*y.l('63');
fc.l('22',compon) = fc.l('21',compon) + fc.l('69',compon);

t.l('22') = t.l('21');
p.l('22') = p.l('21');

* ---------------------------
* - 1 - s p l i t t e r   3 -
* ---------------------------
fc.l('23',compon) = fc.l('22',compon)*y.l('23');
fc.l('24',compon) = fc.l('22',compon)*y.l('24');

t.l('23') = t.l('22');
t.l('24') = t.l('22');
p.l('23') = p.l('22');
p.l('24') = p.l('22');

* ------------------
* - h e a t e r  1 -
* ------------------
fc.l('25',compon) = fc.l('24',compon);
t.l('25') = t.l('24');
p.l('25') = p.l('24');

* ------------------------
* -   c o l u m n - 1    -
* ------------------------
nmin.l('1')   = 2.5*y.l('24');
ndist.l('1')  = nmin.l('1')*4;
rmin.l('1')   = 0.06;
reflux.l('1') = 1.2*rmin.l('1');
distp.l('1')  = 1.02;
avevlt.l('1') = 250.;

fc.l('26','h2')   = fc.l('25','h2') *1.00;
fc.l('26','ch4')  = fc.l('25','ch4')*0.99;
fc.l('26','ben')  = fc.l('25','ben')*0.01;
fc.l('26','tol')  = fc.l('25','tol')*0.00;
fc.l('26','dip')  = fc.l('25','dip')*0.00;
fc.l('27',compon) = fc.l('25',compon) - fc.l('26',compon);

p.l('26') = distp.l('1');
p.l('27') = distp.l('1');
t.l('26') = t.lo('26');
t.l('27') = t.lo('27') + 0.1 ;

vp.l('26','ch4') = (1./7500.6168)*exp(anta('ch4') - antb('ch4')/(t.l('26')*100. + antc('ch4')));
vp.l('26','ben') = (1./7500.6168)*exp(anta('ben') - antb('ben')/(t.l('26')*100. + antc('ben')));
vp.l('27','ch4') = (1./7500.6168)*exp(anta('ch4') - antb('ch4')/(t.l('27')*100. + antc('ch4')));
vp.l('27','ben') = (1./7500.6168)*exp(anta('ben') - antb('ben')/(t.l('27')*100. + antc('ben')));

* ------------------
* - h e a t e r  1 -
* ------------------
fc.l('44',compon) = fc.l('23',compon);
p.l('44') = p.l('23');
t.l('44') = 3.0*(3.24/0.1)**0.23;

* ----------------
* -   e x p   3  -
* ----------------
fc.l('45',compon) = fc.l('44',compon);
p.l('45') = 0.1;
t.l('45') = 3.0;

* ------------------
* - c o o l e r  2 -
* ------------------
fc.l('46',compon) = fc.l('45',compon);
p.l('46') = p.l('45');
t.l('46') = t.l('45');

* ---------------------
* -   f l a s h   2   -
* ---------------------
t.l('47') = t.l('46');
t.l('48') = t.l('46');
p.l('47') = p.l('46');
p.l('48') = p.l('46');

vp.l('48',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.l('48')*100. + antc(compon)));

eflsh.l('2','h2')  = 0.99;
eflsh.l('2','ch4') = 0.99;
eflsh.l('2','ben') = 0.05;
eflsh.l('2','tol') = 0.02;
eflsh.l('2','dip') = 0.0001;

fc.l('47',compon) = fc.l('46',compon)*eflsh.l('2',compon);
fc.l('48',compon) = fc.l('46',compon) - fc.l('47',compon);

* ----------------
* -   e x p   4  -
* ----------------
fc.l('49',compon) = fc.l('47',compon);
p.l('49') = p.l('47');
t.l('49') = t.l('47');

* -------------------------
* -   1 - m i x e r   3   -
* -------------------------
fc.l('30',compon) = fc.l('27',compon) + fc.l('48',compon);
t.l('30') = t.l('27')*y.l('24') + t.l('48')*y.l('23');
p.l('30') = p.l('27')*y.l('24') + p.l('48')*y.l('23');

* ------------------------
* -  s p l i t t e r   3 -
* ------------------------
e.l('28') = 0.9;
e.l('29') = 1.0 - e.l('28');

fc.l('28',compon) = fc.l('26',compon)*e.l('28');
fc.l('29',compon) = fc.l('26',compon)*e.l('29');

p.l('28') = p.l('26');
p.l('29') = p.l('29');
t.l('28') = t.l('26');
t.l('29') = t.l('29');

* ----------------
* -   e x p   5  -
* ----------------
fc.l('50',compon) = fc.l('29',compon);
p.l('50') = p.l('49');
t.l('50') = t.l('29')*(p.l('50')/p.l('29'))**0.23;

* ----------------------
* -    m i x e r   5   -
* ----------------------
fc.l('51',compon) = fc.l('49',compon) + fc.l('50',compon);
p.l('51') = p.l('49');
t.l('51') = t.l('49')*y.l('23') + t.l('50')*y.l('24');

* ------------------------
* -   c o l u m n - 2    -
* ------------------------
nmin.l('2')   = 12.0;
ndist.l('2')  = nmin.l('2')*4;
rmin.l('2')   = 1.3;
reflux.l('2') = 1.2*rmin.l('2');
distp.l('2')  = 0.1;
avevlt.l('2') = 2.5;

p.l('31') = distp.l('2');
p.l('32') = distp.l('2');
t.l('31') = t.lo('31') + 0.1;
t.l('32') = t.lo('32') + 0.1;

fc.l('31','h2')   = fc.l('30','h2') *1.00;
fc.l('31','ch4')  = fc.l('30','ch4')*1.00;
fc.l('31','ben')  = fc.l('30','ben')*0.99;
fc.l('31','tol')  = fc.l('30','tol')*0.01;
fc.l('31','dip')  = fc.l('30','dip')*0.00;
fc.l('32',compon) = fc.l('30',compon) - fc.l('31',compon);

vp.l('31','ben') = ( 1./7500.6168)*exp(anta('ben') - antb('ben')/(t.l('31')*100. + antc('ben')));
vp.l('31','tol') = ( 1./7500.6168)*exp(anta('tol') - antb('tol')/(t.l('31')*100. + antc('tol')));
vp.l('32','ben') = ( 1./7500.6168)*exp(anta('ben') - antb('ben')/(t.l('32')*100. + antc('ben')));
vp.l('32','tol') = ( 1./7500.6168)*exp(anta('tol') - antb('tol')/(t.l('32')*100. + antc('tol')));

* ----------------------------
* - 1 - s p l i t t e r   4  -
* ----------------------------
fc.l('33',compon) = fc.l('32',compon)*y.l('33');
fc.l('37',compon) = fc.l('32',compon)*y.l('37');

t.l('33') = t.l('32');
t.l('37') = t.l('32');
p.l('33') = p.l('32');
p.l('37') = p.l('32');

* ------------------------
* -   c o l u m n - 3    -
* ------------------------
nmin.l('3')   = 2.5*y.l('33');
ndist.l('3')  = nmin.l('3')*4;
rmin.l('3')   = 0.03;
reflux.l('3') = 1.2*rmin.l('3');
distp.l('3')  = 0.1;
avevlt.l('3') = 40.;

fc.l('34','h2')   = fc.l('33','h2') *1.00;
fc.l('34','ch4')  = fc.l('33','ch4')*1.00;
fc.l('34','ben')  = fc.l('33','ben')*1.00;
fc.l('34','tol')  = fc.l('33','tol')*0.99;
fc.l('34','dip')  = fc.l('33','dip')*0.00;
fc.l('35',compon) = fc.l('33',compon) - fc.l('34',compon);

p.l('34') = distp.l('3');
p.l('35') = distp.l('3');
t.l('34') = t.lo('34') + 0.1;
t.l('35') = t.lo('35') + 0.1;

vp.l('34','tol') = (1./7500.6168)*exp(anta('tol') - antb('tol')/(t.l('34')*100. + antc('tol')));
vp.l('34','dip') = (1./7500.6168)*exp(anta('dip') - antb('dip')/(t.l('34')*100. + antc('dip')));
vp.l('35','tol') = (1./7500.6168)*exp(anta('tol') - antb('tol')/(t.l('35')*100. + antc('tol')));
vp.l('35','dip') = (1./7500.6168)*exp(anta('dip') - antb('dip')/(t.l('35')*100. + antc('dip')));

* --------------------
* -  h e a t e r  3  -
* --------------------
fc.l('38',compon) = fc.l('37',compon);
p.l('38') = p.l('37');
t.l('38') = t.l('37');

* ----------------
* -   e x p   2  -
* ----------------
fc.l('39',compon) = fc.l('38',compon);
p.l('39') = 0.1;
t.l('39') = t.l('38');

* ---------------------
* -   f l a s h   3   -
* ---------------------
eflsh.l('2','h2')  = 0.999;
eflsh.l('2','ch4') = 0.999;
eflsh.l('2','ben') = 0.999;
eflsh.l('2','tol') = 0.95;
eflsh.l('2','dip') = 0.25;
fc.l('40',compon)  = fc.l('39',compon) * eflsh.l('2',compon);
fc.l('41',compon)  = fc.l('39',compon) - fc.l('40',compon);

t.l('41') = 4.0;
vp.l('41',compon) = (1/7500.6168)*exp(anta(compon) - antb(compon)/(t.l('41')*100. + antc(compon)));

p.l('40') = p.l('39');
p.l('41') = p.l('39');
t.l('40') = t.l('41');

* ------------------------
* -  s p l i t t e r   2 -
* ------------------------
e.l('52') = 0.2;
e.l('58') = 1.0 - e.l('52');

fc.l('52',compon) = fc.l('18',compon)*e.l('52');
fc.l('58',compon) = fc.l('18',compon)*e.l('58');

p.l('52') = p.l('18');
p.l('58') = p.l('18');
t.l('52') = t.l('18');
t.l('58') = t.l('18');

* ----------------------------
* - 1 - s p l i t t e r   5  -
* ----------------------------
fc.l('53',compon) = fc.l('52',compon)*y.l('53');
fc.l('54',compon) = fc.l('52',compon)*y.l('54');

t.l('53') = t.l('52');
t.l('54') = t.l('52');
p.l('53') = p.l('52');
p.l('54') = p.l('52');

* ------------------------
* - m e m b r a n e   2  -
* ------------------------
fc.l('55','h2')   = fc.l('54','h2') *0.2;
fc.l('55','ch4')  = fc.l('54','ch4')*0.9;
fc.l('55','ben')  = fc.l('54','ben')*1.0;
fc.l('55','tol')  = fc.l('54','tol')*1.0;
fc.l('55','dip')  = fc.l('54','dip')*1.0;
fc.l('56',compon) = fc.l('54',compon) - fc.l('55',compon);

t.l('55') = t.l('54');
t.l('56') = t.l('54');
p.l('55') = p.l('54');
p.l('56') = p.l('54')*0.25;

* ----------------------------
* - 1 - s p l i t t e r   6  -
* ----------------------------
fc.l('61',compon) = fc.l('58',compon)*y.l('61');
t.l('61') = t.l('58');
p.l('61') = p.l('58');

* --------------------
* -  h e a t e r  4  -
* --------------------
fc.l('73',compon)  =  fc.l('61',compon);
p.l('73') = p.l('61');
t.l('73') = t.l('61')*(3.24/0.1)**0.23;

* ----------------
* -   e x p   6  -
* ----------------
fc.l('62',compon) = fc.l('73',compon);
p.l('62') = p.l('51');
t.l('62') = 3.0;

* -----------------
* - m i x e r   4 -
* -----------------
fc.l('63',compon) = fc.l('51',compon) - fc.l('62',compon);
t.l('63') = t.l('62');
p.l('63') = p.l('62');

* -----------------------
* -   a b s o r b e r   -
* -----------------------
f.l('63')          = sum( compon, fc.l('63',compon));
f.l('67')          = f.l('63')*0.20;
fc.l('67',compon)  = f.l('67')*f67comp(compon);
nabs.l(abs)        = 30.*y.l('63');
gamma.l(abs,'tol') = log((1 - aabs('tol')**(nabs.l(abs)*abseff + eps1))/(1 - aabs('tol')));
beta.l(abs,compon) = log((1 - aabs(compon)**(nabs.l(abs)*abseff +  1.0))/(1 - aabs(compon)));

fc.l('64',compon)$asimp('1',compon) =  fc.l('63',compon)/cbeta(compon);
fc.l('64',compon)$anorm('1',compon) =  fc.l('63',compon)/exp( beta.l('1',compon));
fc.l('64',compon)$asolv('1',compon) = (fc.l('63',compon) + fc.l('67',compon)*exp(gamma.l('1',compon)))
                                      /exp(beta.l('1',compon));

fc.l('65','ben')  = fc.l('64','ben');
fc.l('68',compon) = fc.l('63',compon) + fc.l('67',compon) - fc.l('64',compon);
f.l(str)          = sum(compon, fc.l(str,compon));

display f.l, fc.l, t.l, p.l, rctvol.l, conv.l, consum.l, krct.l, sel.l;
display nmin.l, ndist.l, vp.l, gamma.l, beta.l, eflsh.l;

option limCol = 0, limRow = 100, domLim = 1000;

Model douglas / all /;

solve douglas maximizing profit using minlp;

$onText
option limRow = 0;

* y.fx(str)  = 0.0;
* y.fx('2')  = 1.0;
* y.fx('10') = 1.0;
* y.fx('24') = 1.0;
* y.fx('37') = 1.0;
* y.fx('53') = 1.0;
* y.fx('59') = 1.0;

* switch to reactor #2
y.l('10') = 0.0;
y.l('12') = 1.0;

* -----------------------------
* - 1 - s p l i t t e r    2  -
* -----------------------------
fc.l('10',compon) = fc.l('9',compon)*y.l('10');
fc.l('12',compon) = fc.l('9',compon)*y.l('12');

* ---------------------
* -   r e a c t o r   -
* ---------------------
f.l(str)      = sum(compon, fc.l(str,compon));
krct.l(rct)   = 6.3e+10*exp(-26167./(rctt.l(rct)*100.));
rctvol.l('1') = 100.*y.l('10');
rctvol.l('2') = 100.*y.l('12');

conv.l('1','tol') = 1. - sqr(1./(1. + 0.372*krct.l('1')*
                             rctvol.l('1')*sqrt(fc.l('10','tol')/60. + eps1)*
                             (f.l('10')/60. + eps1)**(-3./2.)));

conv.l('2','tol') = 1. - sqr(1./(1.+ 0.372*krct.l('2')*
                             rctvol.l('2')*sqrt(fc.l('12','tol')/60. + eps1)*
                             (f.l('12')/60. + eps1)**(-3./2.)));

sel.l(rct) = 1. - 0.0036*(1. - conv.l(rct,'tol'))**(-1.544);

consum.l('1','tol') = conv.l('1','tol')*fc.l('10','tol');
consum.l('2','tol') = conv.l('2','tol')*fc.l('12','tol');

fc.l('11','tol') = fc.l('10','tol') - consum.l('1','tol');
fc.l('11','ben') = fc.l('10','ben') + consum.l('1','tol')*sel.l('1');
fc.l('11','dip') = fc.l('10','dip') + consum.l('1','tol')*(1 - sel.l('1'))/2.;
fc.l('11','h2')  = fc.l('10','h2')  - consum.l('1','tol')*(1 + sel.l('1'))/2.;
fc.l('11','ch4') = fc.l('10','ch4') + consum.l('1','tol');
fc.l('13','tol') = fc.l('12','tol') - consum.l('2','tol');
fc.l('13','ben') = fc.l('12','ben') + consum.l('2','tol')*sel.l('2');
fc.l('13','dip') = fc.l('12','dip') + consum.l('2','tol')*(1 - sel.l('2'))/2.;
fc.l('13','h2')  = fc.l('12','h2')  - consum.l('2','tol')*(1 + sel.l('2'))/2.;
fc.l('13','ch4') = fc.l('12','ch4') + consum.l('2','tol');

f.l(str)  =  sum(compon, fc.l(str,compon));
t.l('11') = (heatrxn('1')*consum.l('1','tol')/(100.*cp('10')) + f.l('10')*t.l('10'))/(f.l('11') + eps1);
t.l('13') = (heatrxn('2')*consum.l('2','tol')/(100.*cp('12')) + f.l('12')*t.l('12'))/(f.l('13') + eps1);

solve douglas maximizing profit using nlp;

* add membrane #2
y.l('53') = 0.0;
y.l('54') = 1.0;

* ----------------------------
* - 1 - s p l i t t e r   5  -
* ----------------------------
fc.l('53',compon) = fc.l('52',compon)*y.l('53');
fc.l('54',compon) = fc.l('52',compon)*y.l('54');

* ------------------------
* - m e m b r a n e   2  -
* ------------------------
fc.l('55','h2')   = fc.l('54','h2') *0.2;
fc.l('55','ch4')  = fc.l('54','ch4')*0.9;
fc.l('55','ben')  = fc.l('54','ben')*1.0;
fc.l('55','tol')  = fc.l('54','tol')*1.0;
fc.l('55','dip')  = fc.l('54','dip')*1.0;
fc.l('56',compon) = fc.l('54',compon) - fc.l('55',compon);

t.l('55') = t.l('54');
t.l('56') = t.l('54');
p.l('55') = p.l('54');
p.l('56') = p.l('54')*0.25;

solve douglas maximizing profit using nlp;
$offText