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.

References

  • 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.

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

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

$Offtext Eject
$Ontext

 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

======================================================================
$Offtext
$Offsymxref Offsymlist

*======================================================================
* Set Options
*======================================================================
Option limcol=0, limrow=0;
Option bratio=1, domlim=1000;
Option optcr=0;
* reduce default factorization frequency of zoom

*=====================================================================
* Declaration of sets
*=====================================================================
Sets
* 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);
Alias (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) lt 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) eq 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 /;

Parameter yhup(hu,j) current binary proposal for hot utility matches
  /lp.r-1          1 /;

Parameter ycup(i,cu) current binary proposal for cold utility matches
  /c-1.cw   1
   c-3.cw   1 /;

Parameter 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;
Parameter yhuk(hu,j,km)  binary storage parameter yhu;
Parameter ycuk(i,cu,km)  binary storage parameter ycu;
Parameter 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
Parameter 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
Parameter 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/;

Parameter vc(i)  variable charge for distillation columns
    /c-1     0.003375
     c-2     0.000893
     c-3     0.004458
     c-4     0.003176/;

Parameter 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
*=====================================================================

Scalars
   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 soluton 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;
Parameter qrk(j,km)      storage of reboiler duties;
Parameter qck(i,km)      storage of condenser duties;
Parameter qcrk(i,j,km)   storage of heat integrated exchanges;
Parameter qhuk(hu,j,km)  storage of hot utility usage;
Parameter qcuk(i,cu,km)  storage of cold utility usage;
Parameter tck(i,km)      storage of condenser temperatures;
Parameter trk(j,km)      storage of reboiler temperatures;
Parameter lmtdk(i,km)    storage of lmtds;

Scalar zoaup             single value storage of upper bound /inf/;
Scalar 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;
Parameter yhuopt(hu,j)   optimal hot utility match;
Parameter ycuopt(i,cu)   optimal cold utility match;
Parameter 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
*=====================================================================

Variables
        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 Variables
        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 Variables
        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.
*=====================================================================

Equations
    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.
*=====================================================================

nlpobj.. zoau =e=
* capital costs
 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) eq ord(j)
                    and ord(jp) eq 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) eq ord(j)) and (ord(jp) eq 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) eq 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 lt 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 equalties
* to the storage parameter lmtdmar
*======================================================================
  lmtdmar(i,kiter) = -sign(lmtdsn.m(i))$(lmtdsn.m(i) ne 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; );
                                                                   );
* end of file