Description
This GAMS file contains the implementation of the Generalized Benders' Decomposition method for the mathematical formulation based on the superstructure proposed in the paper "Synthesis of General Distillation Sequences - Nonsharp Separations" (Aggarwal and Floudas, 1990). The superstructure features simple columns (single feed and top and bottom products) and no heat integration between columns is considered. The formulation is a mixed-integer nonlinear programming (MINLP) problem. The implementation is based on the algorithmic procedure APROS (Paules and Floudas, 1989).
Small Model of Types : MIP nlp
Category : GAMS Model library
Main file : nonsharp.gms
$title Synthesis of General Distillation Sequences (NONSHARP,SEQ=120)
$onText
This GAMS file contains the implementation of the Generalized Benders'
Decomposition method for the mathematical formulation based on the
superstructure proposed in the paper "Synthesis of General Distillation
Sequences - Nonsharp Separations" (Aggarwal and Floudas, 1990). The
superstructure features simple columns (single feed and top and bottom
products) and no heat integration between columns is considered. The
formulation is a mixed-integer nonlinear programming (MINLP) problem.
The implementation is based on the algorithmic procedure APROS (Paules
and Floudas, 1989).
Morari, M, and Grossmann, I E, Eds, Chemical Engineering Optimization
Models with GAMS. Computer Aids for Chemical Engineering Corporation,
1991.
Aggarwal, A, and Floudas, C A, Synthesis of General Distillation
Sequences: Nonsharp Separation. Computers and Chemical Engineering 14,
6 (1990), 631-653.
Keywords: nonlinear programming, mixed integer linear programming, chemical
engineering, Benders decomposition, distillation sequences, process
synthesis, nonsharp distillation
$offText
$onText
===========================================================================
copyright(C) 1989 A. Aggarwal & C.A. Floudas
Dept. of Chemical Engineering
Princeton University
Princeton, NJ 08544-5263
===========================================================================
-----------------------------
** Three Component Example **
-----------------------------
_____________________________________________________
| Comp. | Feed | Product 1 | Product 2 |
|___________|____________|______________|_____________|
| propane | 100.0 | 30.0 | 70.0 |
| i-butane | 100.0 | 50.0 | 50.0 |
| n-butane | 100.0 | 30.0 | 70.0 |
|___________|____________|______________|_____________|
===========================================================================
$offText
$eolCom //
option limRow = 0, limCol = 0, decimals = 4, solPrint = off;
*==========================================================================
* Declare Sets
*==========================================================================
Set
* Basic Sets - These sets define the problem i.e. number of components
* in the feed stream, columns in the superstructure and no. of product
* streams
cp 'components' / a, b, c /
col 'columns' / col-1*col-2 /
pr 'products' / prod-1*prod-2 /
stm 'column streams' / top, bot /
* Derived Sets - These sets are used to define the mapping of the
* superstructure for the problem at any stage
bypass(pr) 'overall pass to product pr from initial feed'
acol(col) 'existing columns'
prstream(col,stm,pr) 'stream to product pr from column col'
inter(col,col,stm) 'column i fed by column j top or bottom'
key(col,stm,cp) 'key components for column col'
link(col,stm,cp) 'nondistribution of nonkeys'
* Duplicate Sets For Saving Full Mapping - These sets are used for
* saving the full mapping of the original superstructure
savprst(col,stm,pr) 'save full representation of prstream'
savinter(col,col,stm) 'save full representation of feedtop'
savkey(col,stm,cp) 'save full representation of keys'
savlink(col,stm,cp) 'save full representation of link'
* Sets For Algorithm - These are the static and dynamic sets needed
* to implement the iterative algorithm
km 'static iterations' / 1*150 /
kloop 'loop counter' / 1*20 /
k(km) 'dynamic iterations' / 1 /
count(km) 'dynamic counter' / 1 /;
Alias (col,colp);
*==========================================================================
* Define Mapping Of Sets
*==========================================================================
* The set "bypass" defines the set of streams of overall bypasses to the
* final products. If an overall bypass stream for a particular product
* does not exist in the superstructure then that element for this set is
* set to "no"
bypass(pr) = yes;
* The set "inter" defines interconnections between columns based upon the
* superstructure. e.g. "inter('col-1','col-2','top')" refers to a stream
* to column 1 from column 2 top
inter(colp,col,stm) = no;
inter('col-1','col-2','top') = yes;
inter('col-2','col-1','bot') = yes;
* The set "key" defines the key components for each column in the
* superstructure. e.g. "key('col-1','top','a') = yes" defines the light key
* for column 1 to be component a
key(col,stm,cp) = no;
key('col-1','top','a') = yes;
key('col-2','top','b') = yes;
key('col-1','bot','b') = yes;
key('col-2','bot','c') = yes;
* The set "link" is used to restrict distribution of nonkeys in the top
* and bottom streams of a column. Each element of this set represents
* the nonexistence of a particular component in a particular stream
link(col,stm,cp) = no;
link('col-1','top','c') = yes;
link('col-2','bot','a') = yes;
* The set "prstream" defines the set of streams from various columns to
* the final products. If it is known that one or more of such streams
* cannot exist then they can be deleted from the formulation by setting
* the corresponding elements of this set to "no"
prstream(col,stm,pr) = yes;
display bypass, inter, key, link, prstream;
* The full mapping representation of the superstructure is saved in
* the save sets defined earlier
savprst(prstream(col,stm,pr)) = prstream(col,stm,pr);
savinter(inter(colp,col,stm)) = inter(colp,col,stm);
savkey(key(col,stm,cp)) = key(col,stm,cp);
savlink(col,stm,cp) = link(col,stm,cp);
*==========================================================================
* Declare Parameters & Scalars
*==========================================================================
Parameter
* Define the component flowrates for the feedstream }}
feed(cp) 'feed of each component' / a 100, b 100, c 100 /
* Define the various coefficients for the cost expression for each column.
* These coefficients are determined by regression analysis of cost data
* generated by doing column simulations for each column for a range of
* flowrates, feed composition and key component recoveries
cost(col) 'alpha coeff.' / col-1 0.23947, col-2 0.75835 /
a1(col) 'a1 coeff.' / col-1 -0.0139904, col-2 -0.0661588 /
a2(col,stm) 'a2 coeff.' / col-1.top 0.0093514
col-2.top 0.0338147
col-1.bot 0.0077308
col-2.bot 0.0373349 /
a3(col,cp) 'a3 coeff.(corresponds to feed composition)' / col-1.a -0.0005719
col-2.a 0.0016371
col-1.b 0.0042656
col-2.b 0.0288996
col-1.c 0.0
col-2.c 0.0 /
* Define parameters for storage of cost coefficients for each lagrange function
mcost(km,col)
ma1(km,col)
ma2(km,col,stm)
ma3(km,col,cp)
* Define parameters for storage of lagrange multipliers for primal problems
mint(km,col)
* Define parameters for storage of lagrange multipliers for relaxed primal problems
mlint(km,col)
* Define parameters for storage of variable values for primal problems
zk(km)
finitk(km,col)
fink(km,col)
fk(km,col,stm)
fintk(km,colp,col,stm)
fprk(km,col,stm,pr)
fbyk(km,pr)
cfink(km,col,cp)
xkin(km,col,cp)
xk(km,col,stm,cp)
reck(km,col,stm,cp)
* Define parameters to represent current values of complicating variables
* given by the master problem
yp(col)
* Miscellaneous parameters
cut(km)
cutcol(km,col)
nfeas(km) 'feasibility of primal problem'
xinit(cp) 'composition of streams from initial splitter';
Scalar
totfeed 'total flow to the superstructure'
flag 'check to run primal' / 1 /
kend 'stopping criteria for algorithm' / 1 /
iter 'iteration number' / 0 /
zup 'upper bound'
zlo 'lower bound'
np 'number of products';
np = card(pr);
* The table gives the compositions of the desired product streams
Table out(cp,pr) 'product amounts'
prod-1 prod-2
a 30.0 70.0
b 50.0 50.0
c 30.0 70.0;
*==========================================================================
* Declare Variables
*==========================================================================
Variable
c 'total cost of distillation sequence'
alp 'sum of infeasibilities'
mu 'master solution';
Positive Variable
* The positive variables define the material balance for the superstructure
* and the names are self explanatory
finit(col) 'flow from initial splitter to column col'
fin(col) 'total flow to column col'
f(col,stm) 'flow rate of column top & bottom product streams'
fint(colp,col,stm) 'flow rate of interconnecting streams'
fpr(col,stm,pr) 'flow rate of streams to products'
fby(pr) 'flow rate of bypass streams'
cfin(col,cp) 'component flow to column col'
xin(col,cp) 'composition of stream to column col'
x(col,stm,cp) 'composition of column product streams'
rec(col,stm,cp) 'recovery of key components'
* The following are slack variables used for relaxing the constraints of
* the primal problem for the formulation of the relaxed primal problem.
* Each inequality constraint f(x) < 0 is relaxed by the addition of a
* positive slack variables sa to give: f(x) < sa
saint(col) 'sa for logical constraints';
* The binary variables define the existence or nonexistence of a particular
* column in the sequence. A value of 1 means that the column exists,
* a value of 0 denotes nonexistence
Binary Variable y(col) 'existence of column col';
*==========================================================================
* Declare Equations
*==========================================================================
Equation
* PRIMAL PROBLEM EQUATIONS
lpobj 'objective function'
spblinit 'initial splitter balance'
spblcol(col,stm) 'splitter balances at each column top & bottom outlets'
mixbal(col,cp) 'mixer balance equations at inlet of each column'
colbal(col,cp) 'component balance equations around each column'
keybal(col,stm,cp) 'key component balances for each column'
probal(pr,cp) 'component balance equations for each product'
cfloin(col,cp) 'equations defining component flowrates for column inlet'
molsum(col,stm) 'constraint for sum of mole fractions for column products'
molsumin(col) 'constraint for sum of mole fractions for column inlet'
dist(col,stm,cp) 'constraint defining the nondistribution of nonkeys'
intcon(col) 'logical constraints'
* RELAXED PRIMAL PROBLEM EQUATIONS
infeas 'sum of infeasibilities as objective function'
lintcon(col) 'relaxed logical constraints'
* MASTER PROBLEM EQUATIONS
bnd 'lower bound on solution of master problem'
lagrange(km) 'lagrange functions for primal problems'
laerr(km) 'lagrange functions for relaxed primal problems'
intcut(km) 'integer cuts';
*==========================================================================
* Define Equations
*==========================================================================
* PRIMAL PROBLEM EQUATIONS
lpobj.. c =e= sum(acol(col), (cost(col)*yp(col) + (a1(col)
+ sum(key(col,stm,cp), a2(col,stm)*rec(col,stm,cp))
+ sum(cp, a3(col,cp)*xin(col,cp)))*fin(col)));
spblinit.. sum(acol(col), finit(col)) + sum(pr, fby(pr)) =e= totfeed;
spblcol(acol(col),stm)..
sum(inter(colp,col,stm), fint(colp,col,stm))
+ sum(prstream(col,stm,pr), fpr(col,stm,pr)) - f(col,stm)
=e= 0;
mixbal(acol(col),cp)..
finit(col)*xinit(cp) + sum(inter(col,colp,stm), fint(col,colp,stm)*
x(colp,stm,cp)) - cfin(col,cp) =e= 0;
colbal(acol(col),cp)..
cfin(col,cp) - sum(stm,f(col,stm)*x(col,stm,cp)) =e= 0;
keybal(key(col,stm,cp))..
cfin(col,cp)*rec(col,stm,cp) - f(col,stm)*x(col,stm,cp) =e= 0;
probal(pr,cp)$(ord(pr) <> np)..
sum(prstream(col,stm,pr), fpr(col,stm,pr)*x(col,stm,cp)) +
fby(pr)*xinit(cp) - out(cp,pr) =e= 0;
cfloin(acol(col),cp).. fin(col)*xin(col,cp) - cfin(col,cp) =e= 0;
molsum(acol(col),stm).. sum(cp, x(col,stm,cp)) - 1 =e= 0;
molsumin(acol(col)).. sum(cp, xin(col,cp)) - 1 =e= 0;
dist(link(col,stm,cp)).. x(col,stm,cp) =e= 0;
intcon(col).. fin(col) - totfeed*yp(col) =l= 0;
* RELAXED PRIMAL PROBLEM EQUATIONS
infeas.. alp =e= sum(col, saint(col));
lintcon(col).. fin(col) - totfeed*yp(col) =l= saint(col);
* FULL NLP EQUATIONS
* MASTER PROBLEM EQUATIONS
bnd.. mu =g= zlo;
lagrange(k)$(nfeas(k) = 1)..
sum(col, (mcost(k,col)*y(col) + (ma1(k,col)
+ sum(savkey(col,stm,cp), ma2(k,col,stm)*reck(k,col,stm,cp))
+ sum(cp, ma3(k,col,cp)*xkin(k,col,cp)))
* fink(k,col))) // objective function
+ sum(col, mint(k,col)*(fink(k,col) - totfeed*y(col))) =l= mu; // logical const.
laerr(k)$(nfeas(k) = 0)..
sum(col, mlint(k,col)*(fink(k,col) - totfeed*y(col))) =l= 0; // logical const.
intcut(k).. sum(col, cutcol(k,col)*y(col)) =l= cut(k);
* Bounds for all material balance variables
totfeed = sum(cp, feed(cp));
xinit(cp) = feed(cp)/totfeed;
rec.up(col,stm,cp) = 1;
rec.lo(col,stm,cp) = 0.85;
finit.lo(col) = 0;
fin.lo(col) = 0;
f.lo(col,stm) = 0;
fint.lo(savinter(colp,col,stm)) = 0;
fpr.lo(savprst(col,stm,pr)) = 0;
fby.lo(pr) = 0;
finit.up(col) = totfeed;
fin.up(col) = totfeed;
f.up(col,stm) = totfeed;
fint.up(savinter(colp,col,stm)) = totfeed;
fpr.up(savprst(col,stm,pr)) = totfeed;
fby.up('prod-1') = 90;
fby.up('prod-2') = 150;
cfin.lo(col,cp) = 0;
cfin.up(col,cp) = feed(cp);
*==========================================================================
* Declare Models
*==========================================================================
Model
primal / lpobj, spblinit, spblcol, mixbal, colbal, keybal, probal, cfloin,
molsum, molsumin, dist, intcon /
relax / infeas, spblinit, spblcol, mixbal, colbal, keybal, probal, cfloin,
molsum, molsumin, dist, lintcon /
master / bnd, lagrange, laerr, intcut /;
*==========================================================================
* Begin Solving Problem
*==========================================================================
* STARTING POINT - The user has to provide a set of values of the
* complicating variables for the algorithm to start. For finding a
* feasible starting point, the full NLP can be solved (by fixing the
* binary variables). Even though the algorithm is designed to start from
* an infeasible starting point also, it helps the performance to provide
* a feasible starting point
yp('col-1') = 1;
yp('col-2') = 1;
finit.l('col-1') = 60;
finit.l('col-2') = 0;
fin.l('col-1') = 60;
fin.l('col-2') = 40;
f.l('col-1','top') = 20;
f.l('col-2','top') = 20;
f.l('col-1','bot') = 40;
f.l('col-2','bot') = 20;
fint.l('col-1','col-2','top') = 0;
fint.l('col-2','col-1','bot') = 40;
fpr.l('col-1','top','prod-1') = 0;
fpr.l('col-1','top','prod-2') = 20;
fpr.l('col-1','bot','prod-1') = 0;
fpr.l('col-1','bot','prod-2') = 0;
fpr.l('col-2','top','prod-1') = 20;
fpr.l('col-2','top','prod-2') = 0;
fpr.l('col-2','bot','prod-1') = 0;
fpr.l('col-2','bot','prod-2') = 20;
fby.l('prod-1') = 90;
fby.l('prod-2') = 150;
cfin.l('col-1','a') = 20;
cfin.l('col-1','b') = 20;
cfin.l('col-1','c') = 20;
cfin.l('col-2','a') = 0;
cfin.l('col-2','b') = 20;
cfin.l('col-2','c') = 20;
* The upper and lower bounds are initialized to large positive and negative
* numbers respectively
zup = 1000;
zlo = -1000000;
*==========================================================================
* Begin Iterations
*==========================================================================
k(km) = no;
loop(kloop$kend, // START OF LOOP
iter = iter + 1; // update iteration number
count(km) = yes$(ord(km) = iter); // current element
k(km) = k(km) + count(km); // add current element to dynamic set
* Modify sets to current configuration based on current values of the
* binary variables. If a binary variable corresponding to a particular
* column is not active then all streams connected to that column are
* dropped from the superstructure
acol(col) = yes$yp(col);
prstream(savprst(col,stm,pr)) = yes$yp(col);
inter(savinter(colp,col,stm)) = yes$yp(col);
inter(inter(col,colp,stm)) = yes$yp(col);
key(savkey(col,stm,cp)) = yes$yp(col);
link(savlink(col,stm,cp)) = yes$yp(col);
* Store the values of the cost coefficients for the Lagrange function
* based on the current mapping
mcost(count,acol(col)) = cost(col);
ma1(count,acol(col)) = a1(col);
ma2(count,acol(col),stm) = a2(col,stm);
ma3(count,acol(col),cp) = a3(col,cp);
* Solve the relaxed primal problem
solve relax using nlp minimizing alp;
* The next statement checks for the feasiblity of the relaxed primal
* problem and sets the stopping criteria flag "kend" to 0 if the relaxed
* primal problem dees not have a feasible solution
if(((relax.modelStat <> %modelStat.optimal%) and
(relax.modelStat <> %modelStat.locallyOptimal%) and
(relax.modelStat <> %modelStat.feasibleSolution%)), kend = 0);
abort.noerror$(not relax.marginals) "solver did not provide marginals, cannot continue";
* Store the lagrange multipliers as parameters
mlint(count,acol(col)) = -lintcon.m(col);
* Store the values of noncomplicating variables as parameters
finitk(count,acol(col)) = finit.l(col);
fink(count,acol(col)) = fin.l(col);
fk(count,acol(col),stm) = f.l(col,stm);
fintk(count,inter(colp,col,stm)) = fint.l(colp,col,stm);
fprk(count,prstream(col,stm,pr)) = fpr.l(col,stm,pr);
fbyk(count,pr) = fby.l(pr);
cfink(count,acol(col),cp) = cfin.l(col,cp);
xkin(count,acol(col),cp) = xin.l(col,cp);
xk(count,acol(col),stm,cp) = x.l(col,stm,cp);
reck(count,key(col,stm,cp)) = rec.l(col,stm,cp);
nfeas(count) = 0;
flag = 0;
* The value of sum of infeasibilities is checked at this stage
if((alp.l <= .00001),
solve primal using nlp minimizing c ;
* If primal is feasible store the lagrange multipliers as parameters,
* objective function and update upper bound
if(((primal.modelStat = %modelStat.optimal%) or
(primal.modelStat = %modelStat.locallyOptimal%)) and primal.marginals,
mint(count,acol(col)) = -intcon.m(col);
zk(count) = c.l;
nfeas(count) = 1;
if((c.l <= zup), zup = c.l;);
* Store the values of noncomplicating variables as parameters and display them
finitk(count,acol(col)) = finit.l(col);
fink(count,acol(col)) = fin.l(col);
fk(count,acol(col),stm) = f.l(col,stm);
fintk(count,inter(colp,col,stm)) = fint.l(colp,col,stm);
fprk(count,prstream(col,stm,pr)) = fpr.l(col,stm,pr);
fbyk(count,pr) = fby.l(pr);
cfink(count,acol(col),cp) = cfin.l(col,cp);
xkin(count,acol(col),cp) = xin.l(col,cp);
xk(count,acol(col),stm,cp) = x.l(col,stm,cp);
reck(count,key(col,stm,cp)) = rec.l(col,stm,cp);
display finit.l, fin.l, f.l, fint.l, fpr.l, fby.l, cfin.l, xin.l, x.l, rec.l;
);
);
* Store the values of the complicating variables as parameters to formulate
* the integer cuts
cut(count) = sum(col, yp(col)) - 1;
cutcol(count,col) = 2*yp(col) - 1;
* Solve the master problem
solve master using mip minimizing mu;
* The next six statements check for the feasiblity of the master problem
* and set the stopping criteria flag "kend" to 0 if the master problem has
* become infeasible
kend$(master.modelStat = %modelStat.infeasible%) = 0;
kend$(master.modelStat = %modelStat.locallyInfeasible%) = 0;
kend$(master.modelStat = %modelStat.intermediateInfeasible%) = 0;
kend$(master.modelStat = %modelStat.feasibleSolution%) = 0;
kend$(master.modelStat = %modelStat.intermediateNon-Integer%) = 0;
kend$(master.modelStat = %modelStat.integerInfeasible%) = 0;
* If the master problem has a feasible solution, the lower bound isupdated
zlo$kend = mu.l;
* The stopping criteria is checked next. If the upper and lower bounds have
* crossed or are almost equal, "kend" is set to 0
kend$((zup <= zlo) or ((zup - zlo) <= .001)) = 0;
* The values of the complicating variables produced by the master problem
* are stored as parameters (to be used in the next iteration) and displayed
yp(col) = y.l(col);
display yp;
* Display iteration number, sum of infeasibilites for this iteration, and
* the lower and upper bounds at the end of this iteration
display iter, alp.l, zlo, zup;
);