minlphi.gms : Heat Integrated Distillation Sequences

Description

This problem describes a formulation and algorithmic procedure
for obtaining heat-integrated distillation sequences for the separation
of a given multi component feed stream into its pure components products.


Small Model of Types : NLP mip


Category : GAMS Model library


Main file : minlphi.gms

$title Heat Integrated Distillation Sequences (MINLPHI,SEQ=118)

$onText
This problem describes a formulation and algorithmic procedure
for obtaining heat-integrated distillation sequences for the separation
of a given multi component feed stream into its pure components products.


Morari, M, and Grossmann, I E, Eds, Chemical Engineering Optimization
Models with GAMS. Computer Aids for Chemical Engineering Corporation,
1991.

Floudas, C A, and Paules IV, G E, A Mixed-Integer Nonlinear Programming
Formulation for the Synthesis of Heat Integrated Distillation Sequence.
Computers and Chemical Engineering 12, 6 (1988), 531-546.

======================================================================

   A MATHEMATICAL PROGRAMMING FORMULATION FOR PROCESS SYNTHESIS

===================================================================

   copyright    G.E. PAULES IV & C.A. FLOUDAS

            *** Dept. of Chemical Engineering ***
                 *** Princeton University ***
                         May 23, 1987


   Algorithm:  The Outer Approximation with Equality Relaxation
          Full Solution with Starting Point from FIXDT

======================================================================

        This formulation provides the Optimal Heat Integrated
    Distillation Sequence with Pressure as a continuous variable
                for a three component separation.
   The Outer Approximation with Equality Relaxation algorithm is
        used in the automatic solution procedure using GAMS

             Components:     a == Hexane
                             b == Benzene
                             c == Heptane

Total feed to superstructure == 396 kgmol/hr

Multicomponent feed composition:
                        a = 0.80
                        b = 0.10
                        c = 0.10

======================================================================


A Superstructure of the form ...



                       _______               _______
                      _|_    |              _|_    |
                     /   \  ( )            /   \  ( )
                     |   |___|__ A         |   |___|___ B
                     |   |                 |   |
           |---------| 1 |                 | 3 |
           |         |   |       ----------|   |
           |         |   |       |         |   |
           |         |   |_______|         |   |
           |         \___/  |  BC          \___/_______ C
    F      |           |   ( )               |     |
  -------->|           |____|                |----( )
  (ABC)    |
           |           _______               _______
           |          _|_    |              _|_    |
           |         /   \  ( )            /   \  ( )
           |         |   |___| AB          |   |___|___ A
           |         |   |   |_____________|   |
           |---------| 2 |                 | 4 |
                     |   |                 |   |
                     |   |                 |   |
                     |   |______ C         |   |_______ B
                     \___/  |              \___/   |
                       |   ( )               |    ( )
                       |____|                |_____|



 is used with binary variables representing:
    a_  the existence of columns in the sequence.
    b_  the selection of heat exchangers for heat integration.
    c_  the selection of hot and cold utilities.


Associated Reference:
_______________________
 "A Mixed-Integer Nonlinear Programming formulation for the
  synthesis of Heat-Integrated Distillation Sequences"

  C.A. Floudas and G.E. Paules IV,  1988.
  Computers and Chemical Engineering vol. 12 no. 6 pp. 531-546

Keywords: mixed integer linear programming, nonlinear programming, chemical engineering,
          distillation sequences, heat integrated distillation
======================================================================
$offText
$offSymXRef offSymList

*======================================================================
* Set Options
*======================================================================
option limCol = 0, limRow = 0, bRatio = 1, domLim = 1000, optCr = 0;
* reduce default factorization frequency of zoom

*=====================================================================
* Declaration of sets
*=====================================================================
Set
* the set of all columns and their condensers in the superstructure
   i  'condensers-columns' / c-1*c-4 /
* the set of all reboilers in the superstructure
   j  'reboilers' / r-1*r-4 /
* the set of all hot utilities available
   hu 'hot utilities' / lp,ex /
* the set of all cold utilities available
   cu 'cold utilities' / cw /
* an index for linear fit coefficients
   n  'index' / a,b /
* the set of all intermediate products in superstructure
   m  'intermediates' / ab,bc /
* this set maps columns to produced intermediate products
   pm(i,m)  'products' / c-1.bc, c-2.ab /
* this set maps columns to intermediate product feeds
   fm(i,m)  'feeds' / c-3.bc, c-4.ab /
* these sets are for dynamic control of solution algorithm
   km          'static iterations' / k-1*k-100 /
   k(km)       'dynamic iterations'
   kiter(km)   'dynamic counter'
   kdynmax(km) 'dynamic loop control';

* alias sets for condensers and reboilers
Alias (ip,i), (jp,j);
* alias driving loop index - cant appear in equations
Alias (kloop,km);

*=====================================================================
* Definition of "z" parameters for conditional control of model
* used to map permissible matches between condensers and reboilers
* and the position of columns in the superstructure
*=====================================================================

* defines the set of leading columns in the superstructure
Parameter zlead(i) 'leading columns in superstructure' / c-1 1, c-2 1 /;

* defines allowable matches of heat integration for superstructure
* only permits heat integration between columns in the same sequence
Table zcrhx(i,j) 'condenser to reboiler allowable matches'
         r-1  r-2  r-3  r-4
   c-1               1
   c-2                    1
   c-3     1
   c-4          1          ;

* Parameter used in pure integer constraint to permit only one
* direction of heat integration between two columns
* this would yield an infeasible solution but the constraint
* is included explicitly to reduce milp solution time
Parameter zlim(i,j) 'direction of heat integration';
zlim(i,j) = 1$(zcrhx(i,j) and (ord(i) < ord(j)));

* relates appropriate reboiler to the condenser of same column
* (preferably should use an alias rather than a different set)
Parameter zcr(i,j) 'reboiler-condenser pairs';
zcr(i,j) = 1$(ord(i) = ord(j));

*=====================================================================
* Binary variables are divided into 4 classes and variable/parameter
* names starting with "y"
*     ycol - column selection
*     yhx  - heat integration exchanger matches
*     yhu  - hot utility matches
*     ycup - cold utiltiy matches

* These parameters store first guess combination of binary variables
* used to initialize minlp algorithm and parameterize the minlp
* primal problem throughout the rest of the iterations
*=====================================================================
Parameter
   yhxp(i,j)  'current proposal for heat integration matches'    / c-1.r-3  1 /
   yhup(hu,j) 'current binary proposal for hot utility matches'  / lp.r-1   1 /
   ycup(i,cu) 'current binary proposal for cold utility matches' / c-1.cw   1
                                                                   c-3.cw   1 /
   ycolp(i)   'current storage for columns in superstructure'    / c-1      1
                                                                   c-3      1 /;

*=====================================================================
* These parameters store the values of the binary proposals
* for all the iterations performed for use in integer cuts
* and recovering optimal solution
*=====================================================================
Parameter
   yhxk(i,j,km)  'binary storage parameter yhx'
   yhuk(hu,j,km) 'binary storage parameter yhu'
   ycuk(i,cu,km) 'binary storage parameter ycu'
   ycolk(i,km)   'binary storage parameter ycol';

*=====================================================================
* Declaration of parameters for rest of model
*=====================================================================
* mass balances for each sharp separator
Parameter
   spltfrc(i,m) 'split fraction of distillation columns' / c-1.bc 0.20
                                                           c-2.ab 0.90 /
* minimum condenser temperatures obtained from simulation data
   tcmin(i)     'minimum condenser temperatures'         / c-1  341.92
                                                           c-2  343.01
                                                           c-3  353.54
                                                           c-4  341.92 /
* either hottest hot utility-dtmin or for individual separations
* 2*dtmin below critical temperature of bottoms product
   trmax(j)     'maximum reboiler temperatures';

trmax(j) = 1000;

*====================================================================
* scaled cost coefficients for distillation column fits
* nonlinear fixed-charge cost model
* cost = fc*y + vc*flow*temp
* scaling factor = 1000
*====================================================================
Parameter
   fc(i)   'fixed charge for distillation columns'    / c-1  151.125
                                                        c-2  180.003
                                                        c-3  4.2286
                                                        c-4  213.42   /
   vc(i)   'variable charge for distillation columns' / c-1  0.003375
                                                        c-2  0.000893
                                                        c-3  0.004458
                                                        c-4  0.003176 /
   thu(hu) 'hot utility temperatures'                 / lp   421.0
                                                        ex   373.0    /;
* hot utility cost coeff - gives cost in thousands of dollars per year
* ucost = q(10e+6 kj/hr)*costhu(hu)
Parameter costhu(hu) 'hot utility cost coefficients' / lp 24.908, ex 9.139 /;

Table kf(i,n) 'coeff. for heat duty temperature fits'
            a       b
   c-1   32.4  0.0225
   c-2   25.0  0.0130
   c-3   3.76  0.0043
   c-4   35.1  0.0156;

Table af(i,n) 'coeff. for column temperature fits'
             a      b
   c-1   9.541  1.028
   c-2   12.24  1.050
   c-3   8.756  1.029
   c-4   9.181  1.005;

*=====================================================================
* define scalar quantities for rest of model
*=====================================================================
Scalar
   totflow 'total flow to superstructure'               / 396    /
   u       'large number for logical constraints'       / 1500   /
   uint    'upper bound for integer logical'            / 20     /
   fchx    'fixed charge for heat exchangers scaled'    / 3.392  /
   vchx    'variable charge for heat exchangers scaled' / 0.0893 /
   htc     'overall heat transfer coefficient'          / 0.0028 /
   dtmin   'minimum temperature approach'               / 10.0   /
   tcin    'inlet temperature of cold water'            / 305.0  /
   tcout   'outlet temperature of cold water'           / 325.0  /
   costcw  'cooling water cost coefficient'             / 4.65   /
   beta    'income tax correction factor'               / 0.52   /
   alpha   'one over payout time factor in years'       / 0.40   /;

$onText
*=====================================================================
  The parameters declared here are assigned throughout the
  algorithmic procedures.
  They perform the following tasks in the algorithm
    1) transfer of solution data between master and subproblem
    2) storage of solution data
    3) control of upper and lower bounds in milp master
    4) storage of optimal solution
*=====================================================================
$offText

* Storage of variable levels for each iteration
* Identifier derived from name of variable with letter "k" appended
Parameter
   fk(i,km)      'storage of flowrates'
   qrk(j,km)     'storage of reboiler duties'
   qck(i,km)     'storage of condenser duties'
   qcrk(i,j,km)  'storage of heat integrated exchanges'
   qhuk(hu,j,km) 'storage of hot utility usage'
   qcuk(i,cu,km) 'storage of cold utility usage'
   tck(i,km)     'storage of condenser temperatures'
   trk(j,km)     'storage of reboiler temperatures'
   lmtdk(i,km)   'storage of lmtds';

Scalar
   zoaup         'single value storage of upper bound'  /  inf /
   zoalo         'single value storage of lower bounds' / -inf /;

* storage of optimal binary variable combination
* continuous variable levels are not stored separately as they
* can be obtained from the xxxk storage parameters above
Parameter
   yhxopt(i,j)   'optimal heat integration'
   yhuopt(hu,j)  'optimal hot utility match'
   ycuopt(i,cu)  'optimal cold utility match'
   ycolopt(i)    'optimal superstructure';

Scalar kopt      'iteration at which optimal solution was found';

* storage of sign() of Lagrange multiplier from nonlinear equalities
Parameter lmtdmar(i,km) 'direction matrix for nonlinear equalities';

*=====================================================================
* declaration of variables
*=====================================================================
Variable
   zoau         'objective function value of nlp subproblem'
   zoal         'objective function value of milp masters'
   vqcr(km)     'heat integration contribution to milpcon'
   vqhu(km)     'hot utility exchange contribution to milpcon'
   vqcu(km)     'cold utility exchange contribution to milpcon';

Positive Variable
   f(i)         'flowrates to columns'
   qr(j)        'reboiler duties for column with reboiler j'
   qc(i)        'condenser duties for column i'
   qcr(i,j)     'heat integration heat transfer'
   qhu(hu,j)    'hot utility heat transfer'
   qcu(i,cu)    'cold utility heat transfer'
   tc(i)        'condenser temperature for column with cond. i'
   tr(j)        'reboiler temperature for column with reb. j'
   lmtd(i)      'lmtd for cooling water exchanges'
   sl1(i)       'artificial slack variable for lmtd equalities'
   sl2(i)       'artificial slack variable for lmtd equalities'
   s1(i)        'artificial slack variable for reb-con equalities'
   s2(i)        'artificial slack variable for reb-con equalities'
   s3(i)        'artificial slack variable for duty equalities'
   s4(i)        'artificial slack variable for duty equalities';

Binary Variable
   yhx(i,j)     'heat integration matches condenser i reboiler j'
   yhu(hu,j)    'hot utility matches hot utility hu reboiler j'
   ycu(i,cu)    'cold utility matches condenser i cold util cu'
   ycol(i)      'columns in superstructure';

*=====================================================================
* declaration of equations
* for solution of the nlp subproblems:
* early versions of GAMS did not permit binary variables to appear
* in the constraints of a nonlinear programming problem even if
* they appeared in linear constraints and were fixed at a bound
* therefore -
* constraints that contain the binary variables are duplicated:
* one form contains the declared binary variable and the  other
* substitutes a parameter that is assigned the current level of
* the binary variable.  constraints that are duplicated and are to
* appear in the nlp subproblem model have the letter "n" prepended
* to the equation name.
*=====================================================================
Equation
   nlpobj         'nlp subproblems objective'
   milpcon(km)    'nonlinear contribution to milp objective'
   evqcr(km)      'heat integration contribution to milpcon'
   evqhu(km)      'hot utility exchange contribution to milpcon'
   evqcu(km)      'cold utility exchange contribution to milpcon'
   lmtdsn(i)      'nonlinear form of lmtd definition'
   lmtdsm(i,km)   'linearization of lmtdsn(i) in milp masters'
   ntempset(i)    'sets temperatures of inactive columns to 0 (nlp)'
   tempset(i)     'sets temperatures of inactive columns to 0 (milp)'
   nartrex1(i)    'relaxes artificial slack variables (nlp)'
   artrex1(i)     'relaxes artificial slack variables (milp)'
   nartrex2(i)    'relaxes artificial slack variables (nlp)'
   artrex2(i)     'relaxes artificial slack variables (milp)'
   material(m)    'material balances for each intermediate product'
   feed           'feed to superstructure'
   nmatlog(i)     'material balance logical constraints (nlp)'
   matlog(i)      'material balance logical constraints'
   duty(i)        'heat duty definition of condenser i'
   rebcon(i,j)    'equates condenser and reboiler duties'
   conheat(i)     'condenser heat balances'
   rebheat(j)     'reboiler heat balances'
   dtminlp(j)     'minimum temp approach for low pressure steam'
   ndtminc(i)     'minimum temp allowable for each condenser (nlp)'
   dtminc(i)      'minimum temp allowable for each condenser'
   trtcdef(i,j)   'relates reboiler and condenser temps of columns'
   ndtmincr(i,j)  'minimum temp approach for heat integration (nlp)'
   ndtminex(j)    'minimum temp approach for exhaust steam (nlp)'
   nhxclog(i,j)   'logical constraint for heat balances (nlp)'
   nhxhulog(hu,j) 'logical constraint for heat balances (nlp)'
   nhxculog(i,cu) 'logical constraint for heat balances (nlp)'
   nqcqrlog(i)    'logical constraint for con-reb duties (nlp)'
   dtmincr(i,j)   'minimum temp approach for heat integration'
   dtminex(j)     'minimum temp approach for exhaust steam'
   hxclog(i,j)    'logical constraint for heat balances'
   hxhulog(hu,j)  'logical constraint for heat balances'
   hxculog(i,cu)  'logical constraint for heat balances'
   qcqrlog(i)     'logical constraint for con-reb duties'
   boundup        'upper bound on milp objective'
   boundlo        'lower bound on milp objective'

* these are the pure binary constraints of the minlp
   sequen(m)      'restricts superstructure to a single sequence'
   lead           'sequence control'
   limutil(j)     'limits columns to have a single hot utility'
   hidirect(i,j)  'requires a single direction of heat integration'
   heat(i)        'logical integer constraint'
   cuts(km)       'integer cuts for kth iteration';

*=====================================================================
* equations for nlp subproblems
* note that some equations are duplicated in structure but
* given different names in the nlp and milp. these equations
* involve both continuous and binary variables. In older
* versions of GAMS, it was not permissible to pose nonlinear
* models with discrete variables present, even when their values
* were held fixed (rmidnlp). This required two forms of the equation
* two be declared: one with the discrete variables present (milp)
* and one with binary variables replaced by parameters that have
* been assigned the current levels of their associated binary
* variables (nlp). These equations start with the letter "n"
* in the nlp subproblems.
*=====================================================================
*                         capital costs
nlpobj.. zoau =e= alpha*( sum(i,fc(i)*ycolp(i) + vc(i)*(tc(i) - tcmin(i))*f(i))
                        + sum((i,j)$zcrhx(i,j),fchx*yhxp(i,j)
                        + (vchx/htc)*(qcr(i,j)/(tc(i) - tr(j) + 1 - ycolp(i))))
                        + sum((i,cu),fchx*ycup(i,cu)
                        + (vchx/htc)*(qcu(i,cu)/(lmtd(i) + 1 - ycolp(i))))
                        + sum((hu,j),fchx*yhup(hu,j)
                        + (vchx/htc)*(qhu(hu,j)/(thu(hu) - tr(j)))))
*                         operating costs
                        + beta*((costcw*sum((i,cu),qcu(i,cu)))+ sum((hu,j),costhu(hu)*qhu(hu,j)));

lmtdsn(i)..   lmtd(i) - (2/3)*sqrt((tc(i) - tcin)*(tc(i) - tcout))
            - (1/6)*((tc(i) - tcin) + (tc(i) - tcout)) - (sl1(i) - sl2(i)) =e= 0;

nartrex1(i).. s1(i) + s2(i) + sl1(i) - u*(1 - ycolp(i)) =l= 0;

nartrex2(i).. s3(i) + s4(i) + sl2(i) - u*(1 - ycolp(i)) =l= 0;

ntempset(i).. tc(i) + lmtd(i) + sum(j$zcr(i,j),tr(j)) - u*ycolp(i) =l= 0;

material(m).. sum(i$pm(i,m),spltfrc(i,m)*f(i)) - sum(i$fm(i,m),f(i)) =e= 0;

feed..        sum(i$zlead(i),f(i)) =e= totflow;

duty(i).. qc(i) - (kf(i,"a") + kf(i,"b")*(tc(i) - tcmin(i))) - (s3(i) - s4(i)) =e= 0;

rebcon(i,j)$zcr(i,j).. qr(j) - qc(i) =e= 0;

conheat(i).. qc(i) =e= sum(j$zcrhx(i,j),qcr(i,j)) + sum(cu,qcu(i,cu));

rebheat(j).. qr(j) =e= sum(i$zcrhx(i,j),qcr(i,j)) + sum(hu,qhu(hu,j));

trtcdef(i,j)$zcr(i,j)..
   tr(j) - (af(i,"a") + af(i,"b")*(tc(i) - tcmin(i))) - (s1(i) - s2(i)) =e= 0;

nmatlog(i).. f(i) - u*ycolp(i) =l= 0;

ndtminc(i).. (tcmin(i) - tc(i) - u*(1 - ycolp(i))) =l= 0;

dtminlp(j).. dtmin - (thu("lp") - tr(j)) =l= 0;

ndtmincr(i,j)$zcrhx(i,j).. tr(j) - tc(i)  - u*(1 - yhxp(i,j)) + dtmin =l= 0;

ndtminex(j).. dtmin - (thu("ex") - tr(j)) - u*(1 - yhup("ex",j)) =l= 0;

nhxclog(i,j)$zcrhx(i,j).. qcr(i,j) =l= u*yhxp(i,j);

nhxhulog(hu,j).. qhu(hu,j) =l= u*yhup(hu,j);

nhxculog(i,cu).. qcu(i,cu) =l= u*ycup(i,cu);

nqcqrlog(i)..    qc(i) + sum(j$zcr(i,j),qr(j)) - u*ycolp(i) =l= 0;

Model nlpsub  '- collection of equations for nlp subproblems'
              / nlpobj , lmtdsn  , nartrex1, nartrex2, ntempset, material, feed
                nmatlog, duty    , rebcon  , conheat , rebheat , ndtminc , dtminlp
                trtcdef, ndtmincr, ndtminex, nhxclog , nhxhulog, nhxculog, nqcqrlog /;

nlpsub.solPrint=%solPrint.report%;

*======================================================================
* Define equations for milp master problems
* Note: the nonlinear parts of the objective function related
*       to heat exchanger area have been broken out into separate
*       constraints to perform their linearizations, only a
*       contribution term appears in the linearized objective
*       function milpcon.
*======================================================================
milpcon(k)..
   zoal =g= alpha*(sum(i,fc(i)*ycol(i))
         +  fchx*(sum((i,j)$zcrhx(i,j),yhx(i,j))
         +  sum((hu,j),yhu(hu,j)) + sum((i,cu),ycu(i,cu)))
         +  sum(i, (vc(i)*((tck(i,k) - tcmin(i))*(f(i) - fk(i,k))
                 +  fk(i,k)*(tc(i) - tcmin(i)))))
         +  (vchx/htc)*(vqcr(k) + vqhu(k) + vqcu(k)))
         +  beta*((costcw*sum((i,cu),qcu(i,cu)))
         +  sum((hu,j),costhu(hu)*qhu(hu,j)));

*==========================================================================
* these are the linearized contributions to the objective related
* to heat exchange.  the appearance of the binary variable storage
* parameters in the denominator of some of the expressions is done
* to prevent division by zero during model generation for linearizations
* done at points where the temperatures were set to zero for unused
* columns.  the numerator is zero then also and no error is introduced.
*==========================================================================
evqcr(k)..
   vqcr(k) =e= sum((i,j)$zcrhx(i,j),((qcrk(i,j,k)/(tck(i,k) - trk(j,k) + 1 - ycolk(i,k)))
                    + ((1/(tck(i,k) - trk(j,k) + 1 - ycolk(i,k)))
                    * (qcr(i,j) - qcrk(i,j,k)))*ycolk(i,k)
                    + ((qcrk(i,j,k)/(sqr(tck(i,k) - trk(j,k)) + 1 - ycolk(i,k)))
                    * ((tr(j) - trk(j,k)) - (tc(i) - tck(i,k))))));

evqhu(k)..
   vqhu(k) =e= sum((hu,j),((qhuk(hu,j,k)/(thu(hu) - trk(j,k)))
                    + ((1/(thu(hu) - trk(j,k)))
                    * (qhu(hu,j) - qhuk(hu,j,k)))*sum(i$zcr(i,j), ycolk(i,k))
                    + ((qhuk(hu,j,k)/sqr(thu(hu) - trk(j,k)))*(tr(j) - trk(j,k)))));

evqcu(k)..
   vqcu(k) =e= sum((i,cu),((qcuk(i,cu,k)/(lmtdk(i,k) + 1 - ycolk(i,k)))
                    + ((1/(lmtdk(i,k) + 1 - ycolk(i,k)))
                    * (qcu(i,cu) - qcuk(i,cu,k)))*ycolk(i,k)
                    - ((qcuk(i,cu,k)/(sqr(lmtdk(i,k)) + 1 - ycolk(i,k)))
                    * (lmtd(i) - lmtdk(i,k)))));

lmtdsm(i,k)..
     lmtdmar(i,k)*(lmtd(i) - (2/3)*sqrt((tck(i,k) - tcin)*(tck(i,k) - tcout))
   - (1/6)*((tck(i,k) - tcin) + (tck(i,k) - tcout))
   - ((1/3)*(((2*tck(i,k) - (tcin + tcout))
   / sqrt(sqr(tck(i,k)) - (tcin + tcout)*tck(i,k) + (tcin*tcout))) + 1))
   * (tc(i) - tck(i,k)) - (sl1(i) - sl2(i))) =l= 0;

artrex1(i).. s1(i) + s2(i) + sl1(i) - u*(1 - ycol(i)) =l= 0;

artrex2(i).. s3(i) + s4(i) + sl2(i) - u*(1 - ycol(i)) =l= 0;

tempset(i).. tc(i) + lmtd(i) + sum(j$zcr(i,j), tr(j)) - u*ycol(i) =l= 0;

matlog(i).. f(i) - u*ycol(i) =l= 0;

dtminc(i).. (tcmin(i) - tc(i) - u*(1 - ycol(i))) =l= 0;

dtmincr(i,j)$zcrhx(i,j).. tr(j) - tc(i) - u*(1 - yhx(i,j)) + dtmin =l= 0;

dtminex(j).. dtmin - (thu("ex") - tr(j)) - u*(1 - yhu("ex",j)) =l= 0;

hxclog(i,j)$zcrhx(i,j).. qcr(i,j) =l= u*yhx(i,j);

hxhulog(hu,j).. qhu(hu,j) =l= u*yhu(hu,j);

hxculog(i,cu).. qcu(i,cu) =l= u*ycu(i,cu);

qcqrlog(i).. qc(i) + sum(j$zcr(i,j), qr(j)) - u*ycol(i) =l= 0;

* pure binary constraints
* material balances determine sequence
sequen(m).. sum(i$pm(i,m), ycol(i)) - sum(i$fm(i,m), ycol(i)) =e= 0;

* select 1 sequence
lead.. sum(i$zlead(i),ycol(i)) =e= 1;

* limit choice of hot utility to 1
limutil(j).. sum(hu,yhu(hu,j)) =l= 1;

* only one of the mutual heat integration binaries can be 1
hidirect(i,j)$zlim(i,j).. yhx(i,j) + sum((ip,jp)$(ord(ip) = ord(j) and
                                                  ord(jp) = ord(i)), yhx(ip,jp)) =l= 1;

* if a column doesn't exist then all binary variables associated
* with it must also be set to zero
heat(i)..    sum(j$zcrhx(i,j),yhx(i,j)
          +  sum((ip,jp)$((ord(ip) = ord(j)) and (ord(jp) = ord(i))),yhx(ip,jp)))
          +  sum((hu,j),yhu(hu,j)$zcr(i,j))
          +  sum(cu,ycu(i,cu)) - uint*ycol(i)
         =l= 0;

* integer cuts
cuts(k)..    sum(i,sign(ycolk(i,k) - 0.5)*ycol(i))
          +  sum((i,j)$zcrhx(i,j),sign(yhxk(i,j,k) - 0.5)*yhx(i,j))
          +  sum((hu,j),sign(yhuk(hu,j,k) - 0.5)*yhu(hu,j))
          +  sum((i,cu),sign(ycuk(i,cu,k) - 0.5)*ycu(i,cu))
         =l= sum(i,ycolk(i,k)) + sum((i,j)$zcrhx(i,j),yhxk(i,j,k))
          +  sum((hu,j),yhuk(hu,j,k)) + sum((i,cu),ycuk(i,cu,k)) - 1;

*======================================================================
* declare the milp master problem
*======================================================================
Model master 'milp master problem'
             / milpcon, evqcr  , evqhu   , evqcu   , lmtdsm , artrex1
               artrex2, tempset, material, feed    , matlog , duty
               rebcon , conheat, rebheat , dtminc  , dtminlp, trtcdef
               dtmincr, dtminex, hxclog  , hxhulog , hxculog, qcqrlog
               sequen , lead   , limutil , hidirect, heat   , cuts    /;

master.solPrint=%solPrint.summary%;

*=====================================================================
* all declarations made, start algorithmic procedures

* initialize the optimal storage parameters to 1st guess
*=====================================================================
yhxopt(i,j)  = yhxp(i,j);
yhuopt(hu,j) = yhup(hu,j);
ycuopt(i,cu) = ycup(i,cu);
ycolopt(i)   = ycolp(i);
kopt         = 1;

*======================================================================
* assign the initial configuration to the binary proposal parameter
*======================================================================
kiter("k-1") = yes;

yhxk(i,j,kiter)  = yhxp(i,j);
yhuk(hu,j,kiter) = yhup(hu,j);
ycuk(i,cu,kiter) = ycup(i,cu);
ycolk(i,kiter)   = ycolp(i);
yhx.l(i,j)       = yhxp(i,j);
yhu.l(hu,j)      = yhup(hu,j);
ycu.l(i,cu)      = ycup(i,cu);
ycol.l(i)        = ycolp(i);

* set an arbitrary initial lower bound
zoal.l = -10e+6;

*======================================================================
* give the continuous variables a starting point for 1st nlp
*======================================================================
tr.l("r-1")   = 410;
tc.l("c-1")   = 390;
tc.l("c-3")   = 360;
tr.l("r-3")   = 380;
tc.l("c-2")   = 0;
tr.l("r-2")   = 0;
tc.l("c-4")   = 0;
tr.l("r-4")   = 0;
f.l("c-1")    = totflow;
lmtd.l("c-1") = 75;
lmtd.l("c-3") = 25;
lmtd.l("c-2") = 0;
lmtd.l("c-4") = 0;
qr.l("r-2")   = 0;
qc.l("c-2")   = 0;
qr.l("r-4")   = 0;
qc.l("c-4")   = 0;

*======================================================================
* add bounds on tc. A sqrt in equation lmtdsn is defined for tc > tcout
* and for tc < tcin. The relevant interval is determined for each
* element of tc based on the initial values given above.
*======================================================================
tc.lo("c-1") = tcout+1;
tc.up("c-2") = tcin-1;
tc.lo("c-3") = tcout+1;
tc.up("c-4") = tcin-1;

*======================================================================
* bound the reboiler temperatures by their maximum allowable
*======================================================================
tr.up(j) = trmax(j);

*======================================================================
* initialize the dynamic sets for algorithm control
*======================================================================
k(km)       =  no;
kiter(km)   =  no;
kdynmax(km) = yes;

*======================================================================
* major driving loop of algorithm
*======================================================================
loop(kloop$kdynmax(kloop),

*  update the dynamic iteration sets
* -set kiter to contain only the current iteration element
* -add to k the current iteration element
   kiter(km) = yes$(ord(km) = ord(kloop));
   k(kiter)  = yes;

*  store the current binary combination
   yhxk(i,j,kiter)  = yhx.l(i,j);
   yhuk(hu,j,kiter) = yhu.l(hu,j);
   ycuk(i,cu,kiter) = ycu.l(i,cu);
   ycolk(i,kiter)   = ycol.l(i);

*  set the current combination parameters that appear in the nlp constraints
   yhxp(i,j)  = yhx.l(i,j);
   yhup(hu,j) = yhu.l(hu,j);
   ycup(i,cu) = ycu.l(i,cu);
   ycolp(i)   = ycol.l(i);
   zoal.lo    = zoal.l;

*======================================================================
* the current levels of the lmtds are moved away from zero
* to prevent evaluation errors in the next nlp subproblem
*======================================================================
   lmtd.l(i) = lmtd.l(i) + 1;

   nlpsub.resLim = 15;
*  solve the nlp subproblem
   solve nlpsub using nlp minimizing zoau;

*  resolve with Conopt to get marginals for lmtdsn, if not provided by used NLP solver
   if(nlpsub.marginals = 0,
      option nlp = conopt;
      solve nlpsub using nlp minimizing zoau;
      option nlp = %system.nlp%;
   );

*======================================================================
* update the optimal solution storage parameters if new nlp
* objective function value is less than the incumbent
*======================================================================
   if((zoau.l < zoaup),
      yhxopt(i,j)  = yhx.l(i,j);
      yhuopt(hu,j) = yhu.l(hu,j);
      ycuopt(i,cu) = ycu.l(i,cu);
      ycolopt(i)   = ycol.l(i);
      kopt         = ord(kloop);
   );

*======================================================================
* assign the solution levels of the variables that appear in the
* nonlinear equations to their corresponding storage parameters
*======================================================================
   fk(i,kiter)      = f.l(i);
   qrk(j,kiter)     = qr.l(j);
   qck(i,kiter)     = qc.l(i);
   qcrk(i,j,kiter)  = qcr.l(i,j);
   qhuk(hu,j,kiter) = qhu.l(hu,j);
   qcuk(i,cu,kiter) = qcu.l(i,cu);
   tck(i,kiter)     = tc.l(i);
   trk(j,kiter)     = tr.l(j);
   lmtdk(i,kiter)   = lmtd.l(i);

*======================================================================
* assign the sign of marginal values of the nonlinear equalities
* to the storage parameter lmtdmar
*======================================================================
   lmtdmar(i,kiter) = -sign(lmtdsn.m(i))$(lmtdsn.m(i) <> eps);

*======================================================================
* store the smallest nlp objective value for upper bound on master
*======================================================================
   zoaup   = min(zoaup,zoau.l);
   zoal.up = zoaup;
*  protect against numerical errors introduced by the solver
   zoal.lo = min(zoal.lo,zoal.up);

*  now solve the milp master problem
   solve master using mip minimizing zoal;

   display "new binary combination", ycol.l, yhx.l, yhu.l, ycu.l;

*======================================================================
* check stopping criterion:
* master problem integer infeasible
*======================================================================
   if((master.modelStat=%modelStat.infeasible% or
       master.modelStat=%modelStat.integerInfeasible% or
       master.modelStat=%modelStat.infeasibleNoSolution%),
      kdynmax(km) = no;
      display "stopping criterion met", zoaup, yhxopt, yhuopt, ycuopt, ycolopt, kopt;
   );
);