Description
This model determines the optimum location of feed plates in a distillation column with 7 ideal stages, a total condenser and a kettle-type reboiler. The feed consists of a mixture of benzene and toluene entering at bubble point.
Large Model of Type : MINLP
Category : GAMS Model library
Main file : feedtray.gms
$title Optimum Feed Plate Location (FEEDTRAY,SEQ=122)
$onText
This model determines the optimum location of feed plates in
a distillation column with 7 ideal stages, a total condenser and
a kettle-type reboiler. The feed consists of a mixture of benzene
and toluene entering at bubble point.
Morari, M, and Grossmann, I E, Eds, Chemical Engineering Optimization
Models with GAMS. Computer Aids for Chemical Engineering Corporation,
1991.
Vishwanathan, J, and Grossmann, I E, A Combined Penalty Function and
Outer Approximation Method for Mixed Integer Nonlinear Programming.
Computers and Chemical Engineering 14, 7 (1990), 769-782.
Section 1
Basic Data
___________
Basic Thermodynamic Data have been taken from
R.C. Reid, J.M. Prausnitz and B.E. Poling :
"The Properties of gases and liquids"
4th Edition, McGraw-Hill (1987)
Units: Pressure : in bars i.e. 0.1 Mpa
Temperature : in kelvins
Heat capacity : kJ/kmol K
Abbreviations: mm Molar mass
Tb Normal boiling point
Tc Critical temperature
Pc Critical Pressure
Keywords: mixed integer nonlinear programming, feed plate location, chemical
engineering
$offText
$sTitle System --- Benzene - Toluene at about 1 bar
Set j 'components' / benzene, toluene /;
Table prcon (j,*) 'basic physical properties'
mm tb tc pc omega liq-den tden
benzene 78.114 353.2 562.2 48.9 0.212 0.885 289
toluene 92.141 383.8 591.8 41.0 0.263 0.867 293;
Table vpcon(j,*) 'constants in the weird equation for vap.pressure'
a b c d tmin
benzene -6.98273 1.33213 -2.62863 -3.33399 288
toluene -7.28607 1.38091 -2.83433 -2.79168 309;
Table cpcon(j,*) 'constants for the isobaric heat capacity equation'
a b c d
benzene -3.392e+1 4.739e-1 -3.017e-4 7.130e-8
toluene -2.435e+1 5.125e-1 -2.765e-4 4.911e-8;
Scalar
treb 'guess temperature for reboiler'
tbot 'guess temperature for bottom-most tray'
ttop 'guess temperature for top-most tray'
tcon 'guess temperature for condenser';
treb = 380;
tbot = 375;
ttop = 360;
tcon = 355;
Scalar rg 'universal gas constant' / 8.314 /;
Parameter h0(j) 'integration constant for enthalpy';
h0(j) = tcon*(cpcon(j,'a') + tcon*(cpcon(j,'b')/2
+ tcon*(cpcon(j,'c')/3 + tcon* cpcon(j,'d')/4)));
$onText
section 2
_________
thermal condition of the feed stream
$offText
Scalar
f 'total no. of moles of feed'
tf 'temperature of the feed (in kelvins)'
pf 'pressure of the feed stream'
vf 'vapour fraction in feed(before expansion)'
shf 'specific enthalpy of feed'
preb 'pressure in the reboiler'
pbot 'pressure in the bottom-most tray'
ptop 'pressure in the top-most tray'
pcon 'pressure in the condenser';
f = 100;
preb = 1.2;
pbot = 1.12;
ptop = 1.08;
pcon = 1.05;
pf = 1.12;
vf = 0.0;
Parameter xf(j) 'molefractions in feed-stream' / benzene 0.70, toluene 0.30 /;
$onText
Assume that feed enters at bubble point determination of bubble temperature
use nlp to solve equation commented out for convenience.
Variable tbub, z1;
Equation bubble, obj1;
bubble.. sum(j, xf(j)*prcon(j,'pc')*exp(prcon(j,'tc')/tbub
* (vpcon(j,'a')*(1 - tbub/prcon(j,'tc'))
+ vpcon(j,'b')*(1 - tbub/prcon(j,'tc'))**1.5
+ vpcon(j,'c')*(1 - tbub/prcon(j,'tc'))**3
+ vpcon(j,'d')*(1 - tbub/prcon(j,'tc'))**6))/pf)
=e= 1;
obj1.. z1 =e= 1;
tbub.lo = prcon('benzene','tb');
tbub.up = prcon('toluene','tb');
tbub.l = 0.5*(tbub.lo +tbub.up);
Model bubt / all /;
solve bubt minimizing z1 using nlp;
$offText
* Note: Solution of the above step indicates tbub = 363.368
tf = 363.368;
shf = sum(j, xf(j)*(tf*(cpcon(j,'a') + tf*(cpcon(j,'b')/2
+ tf*(cpcon(j,'c')/3 + tf*cpcon(j,'d')/4))) - h0(j)
+ rg* prcon(j,'tc')*(vpcon(j,'a')*(1 - tf/prcon(j,'tc'))
+ vpcon(j,'b')*(1 - tf/prcon(j,'tc'))**1.5
+ vpcon(j,'c')*(1 - tf/prcon(j,'tc'))**3
+ vpcon(j,'d')*(1 - tf/prcon(j,'tc'))**6)
+ rg*tf*(vpcon(j,'a')
+ 1.5*vpcon(j,'b')*(1 - tf/prcon(j,'tc'))**0.5
+ 3*vpcon(j,'c')*(1 - tf/prcon(j,'tc'))**2
+ 6*vpcon(j,'d')*(1 - tf/prcon(j,'tc'))**5)));
display shf;
$onText
section 3
_________
modeling equations
description of the column
Note: the stages are numbered bottom upwards (like the
floors of a building). Reboiler is stage (tray) no. 1 and the
condenser is the last tray.
$offText
Set
i 'stages' / 1*9 /
reb(i) 'reboiler'
con(i) 'condenser'
col(i) 'stages in the col'
floc(i) 'possible locations of feed stage' / 2*8 /
abovef(i) 'stages above the feed stage'
belowf(i) 'feed stage and those below it';
reb(i) = yes$(ord(i) = 1);
con(i) = yes$(ord(i) = card(i));
col(i) = yes - (reb(i) + con(i));
Parameter p(i) 'pressure prevailing in tray i';
p(i)$reb(i) = preb;
p(i)$con(i) = pcon;
p(i)$col(i) = pbot - ((pbot - ptop)/(card(i) - 1 - 2))*(ord(i) - 2);
Scalar
hllo 'limit on enthalpies and scaling'
hlhi 'limit on enthalpies and scaling'
hvlo 'limit on enthalpies and scaling'
hvhi 'limit on enthalpies and scaling'
hscale 'limit on enthalpies and scaling';
hllo = tcon*(cpcon('benzene','a') + tcon*(cpcon('benzene','b')/2
+ tcon*( cpcon('benzene','c')/3 + tcon* cpcon('benzene','d')/4)))
- h0('benzene') + rg*prcon('benzene','tc')
* (vpcon('benzene','a')*(1 - tcon/prcon('benzene','tc'))
+ vpcon('benzene','b')*(1 - tcon/prcon('benzene','tc'))**1.5
+ vpcon('benzene','c')*(1 - tcon/prcon('benzene','tc'))**3
+ vpcon('benzene','d')*(1 - tcon/prcon('benzene','tc'))**6)
+ rg*tcon*(vpcon('benzene','a')
+ 1.5*vpcon('benzene','b')*(1 - tcon/prcon('benzene','tc'))**0.5
+ 3*vpcon('benzene','c')*(1 - tcon/prcon('benzene','tc'))**2
+ 6*vpcon('benzene','d')*(1 - tcon/prcon('benzene','tc'))**5);
hlhi = treb*(cpcon('toluene','a') + treb*(cpcon('toluene','b')/2
+ treb*(cpcon('toluene','c')/3 + treb*cpcon('toluene','d')/4)))
- h0('toluene') + rg*prcon('toluene','tc')
* (vpcon('toluene','a')*(1 - treb/prcon('toluene','tc'))
+ vpcon('toluene','b')*(1 - treb/prcon('toluene','tc'))**1.5
+ vpcon('toluene','c')*(1 - treb/prcon('toluene','tc'))**3
+ vpcon('toluene','d')*(1 - treb/prcon('toluene','tc'))**6)
+ rg*treb*(vpcon('toluene','a')
+ 1.5*vpcon('toluene','b')*(1 - treb/prcon('toluene','tc'))**0.5
+ 3*vpcon('toluene','c')*(1 - treb/prcon('toluene','tc'))**2
+ 6*vpcon('toluene','d')*(1 - treb/prcon('toluene','tc'))**5);
hvlo = tcon*(cpcon('benzene','a') + tcon*(cpcon('benzene','b')/2
+ tcon*(cpcon('benzene','c')/3 + tcon*cpcon('benzene','d')/4)))
- h0('benzene');
hvhi = treb*(cpcon('toluene','a') + treb*(cpcon('toluene','b')/2
+ treb*(cpcon('toluene','c')/3 + treb*cpcon('toluene','d')/4)))
- h0('toluene');
hscale = max(abs(hllo), abs(hlhi), abs(hvlo), abs(hvhi));
Positive Variable
x(i,j) 'mole-fraction of j-th component in liquid on i-th tray'
y(i,j) 'mole-fraction of j-th component in vapour on i-th tray'
l(i) 'molar flow rate of liquid leaving tray i'
v(i) 'molar flow rate of vapour leaving tray i'
t(i) 'temperature of tray i'
feed(i) 'feed stream entering tray i'
r 'reflux ratio'
p1 'top product rate'
p2 'bottom product rate';
Variable
hl(i) 'molar sp.enthalpy of liquid in tray i'
hv(i) 'molar sp.enthalpy of vapour in tray i';
Equation
phe(i,j) 'phase equilibrium relation'
errk(i) 'phase equilibrium error function'
cmb(i,j) 'component material balance(1<i<n)'
cmb1(i,j) 'component material balance(i=1)'
cmbn(i,j) 'component material balance(i=n)'
tmb(i) 'total material balance for trays in the column'
tmb1(i) 'total material balance for reboiler'
tmbn(i) 'total material balance for condenser'
defln(i) 'definition of l(n)'
defp2(i) 'definition of p2'
defhl(i) 'definition of hl(i)'
defhv(i) 'definition of hv(i)'
eb(i) 'enthalpy balance'
purcon 'purity constraint'
sumf 'sum of feeds';
phe(i,j).. y(i,j)*p(i) - x(i,j)*prcon(j,'pc')*exp(prcon(j,'tc')/t(i)
* (vpcon(j,'a')*(1 - t(i)/prcon(j,'tc'))
+ vpcon(j,'b')*(1 - t(i)/prcon(j,'tc'))**1.5
+ vpcon(j,'c')*(1 - t(i)/prcon(j,'tc'))**3
+ vpcon(j,'d')*(1 - t(i)/prcon(j,'tc'))**6))
=e= 0;
errk(i).. sum(j,y(i,j)) - sum(j,x(i,j)) =e= 0;
cmb(i,j)$col(i).. l(i)*x(i,j) + v(i)*y(i,j) - l(i+1)*x(i+1,j)
- v(i-1)*y(i-1,j) - (feed(i)*xf(j))$floc(i) =e= 0;
cmb1(i,j)$reb(i).. l(i)*x(i,j) + v(i)*y(i,j) - l(i+1)*x(i+1,j) =e= 0;
cmbn(i,j)$con(i).. (l(i)+ p1)*x(i,j) - v(i-1)*y(i-1,j) =e= 0;
tmb(i)$col(i).. l(i) + v(i) - l(i+1) - v(i-1) - feed(i)$floc(i) =e= 0;
tmb1(i)$reb(i).. l(i) + v(i) - l(i+1) =e= 0;
tmbn(i)$con(i).. l(i) + p1 - v(i-1) =e= 0;
defln(i)$con(i).. l(i) =e= r*p1;
defp2(i)$reb(i).. l(i) - p2 =e= 0;
defhl(i).. hl(i) - sum(j, x(i,j)*(t(i)*(cpcon(j,'a') + t(i)*(cpcon(j,'b')/2
+ t(i)*(cpcon(j,'c')/3 + t(i)*cpcon(j,'d')/4))) - h0(j)
+ rg*prcon(j,'tc')*(vpcon(j,'a')*(1 - t(i)/prcon(j,'tc'))
+ vpcon(j,'b')*(1 - t(i)/prcon(j,'tc'))**1.5
+ vpcon(j,'c')*(1 - t(i)/prcon(j,'tc'))**3
+ vpcon(j,'d')*(1 - t(i)/prcon(j,'tc'))**6)
+ rg*t(i)*(vpcon(j,'a')
+ 1.5*vpcon(j,'b')*(1 - t(i)/prcon(j,'tc'))**0.5
+ 3 *vpcon(j,'c')*(1 - t(i)/prcon(j,'tc'))**2
+ 6 *vpcon(j,'d')*(1 - t(i)/prcon(j,'tc'))**5)))/hscale
=e= 0;
defhv(i).. hv(i) - sum(j, y(i,j)*(t(i)*(cpcon(j,'a') + t(i)*( cpcon(j,'b')/2
+ t(i)*(cpcon(j,'c')/3 + t(i)*cpcon(j,'d')/4))) - h0(j)))/hscale
=e= 0;
eb(i)$col(i).. l(i)*hl(i) + v(i)*hv(i) - l(i+1)*hl(i+1) - v(i-1)*hv(i-1)
- (feed(i)*shf/hscale)$floc(i) =e= 0;
purcon.. x('9','benzene') =g= 0.95;
sumf.. sum(i$floc(i), feed(i)) =e= f;
* Initialization of variables
Set ifeed(i) / 6 /;
belowf(i) = yes$(ord(i) le 6);
abovef(i) = yes - belowf(i);
feed.l(i)$ifeed(i) = f;
r.l = 0.45; r.up = 0.95; r.lo = 0.1;
p1.l = 60 ; p1.up = 80 ; p1.lo = 50;
p2.l = 40 ; p2.up = 50 ; p2.lo = 20;
x.up(i,j) = 1.0;
y.up(i,j) = 1.0;
v.l(i) = ((r.l + 1)*p1.l);
l.l(i)$reb(i) = p2.l;
l.l(i)$(belowf(i)- reb(i)) = p1.l*r.l + (1 - vf)*f;
l.l(i)$abovef(i) = p1.l*r.l;
t.l(i)$reb(i) = treb;
t.l(i)$con(i) = tcon;
t.l(i)$col(i) = tbot + (ttop - tbot)*(ord(i) - 2)/(card(i) - 1 - 2);
t.lo(i) = tcon - 10;
t.up(i) = treb + 10;
x.l(i,'benzene') = 0.90 - 0.5*(ord(i) - card(i))/(1 - card(i));
x.l(i,'toluene') = 1 - x.l(i,'benzene');
y.l(i,'benzene') = 1.0 - 0.8*(ord(i) - card(i))/(1 - card(i));
y.l(i,'toluene') = 1 - y.l(i,'benzene');
hl.l(i) = sum(j, x.l(i,j)*(t.l(i)*(cpcon(j,'a') + t.l(i)*(cpcon(j,'b')/2
+ t.l(i)*( cpcon(j,'c')/3 + t.l(i) *cpcon(j,'d')/4))) - h0(j)
+ rg*prcon(j,'tc')*(vpcon(j,'a')*(1 - t.l(i)/prcon(j,'tc'))
+ vpcon(j,'b')*(1 - t.l(i)/prcon(j,'tc'))**1.5
+ vpcon(j,'c')*(1 - t.l(i)/prcon(j,'tc'))**3
+ vpcon(j,'d')*(1 - t.l(i)/prcon(j,'tc'))**6)
+ rg*t.l(i)*(vpcon(j,'a')
+ 1.5*vpcon(j,'b')*(1 - t.l(i)/prcon(j,'tc'))**0.5
+ 3 *vpcon(j,'c')*(1 - t.l(i)/prcon(j,'tc'))**2
+ 6 *vpcon(j,'d')*(1 - t.l(i)/prcon(j,'tc'))**5)))/hscale;
hv.l(i) = sum(j, y.l(i,j)*(t.l(i)*(cpcon(j,'a') + t.l(i)*(cpcon(j,'b')/2
+ t.l(i)*(cpcon(j,'c')/3 + t.l(i)*cpcon(j,'d')/4))) - h0(j)))/hscale;
hl.lo(i) = (1 - 0.5*sign(hllo))*hllo/ hscale;
hv.lo(i) = (1 - 0.5*sign(hvlo))*hvlo/ hscale;
hl.up(i) = (1 + 0.5*sign(hlhi))*hlhi/ hscale;
hv.up(i) = (1 + 0.5*sign(hvhi))*hvhi/ hscale;
* Additional variables and equations for the feed tray location problem;
Binary Variable yf(i);
Variable zf;
Equation
sumb 'sum of binary variables'
confeed(i) 'constraint on feed'
obj2 'second objective function';
sumb.. sum(i$floc(i),yf(i)) =e= 1;
confeed(i)$floc(i).. feed(i) - f*yf(i) =l= 0;
obj2.. zf =e= p1 - 50*r;
zf.up = 100;
yf.l(i)$floc(i) = 1.0/7.0;
feed.l(i)$floc(i) = f/7.0;
feed.up(i)$floc(i) = f;
Model column 'optimization case' / all /;
solve column using minlp maximizing zf;