* ==>hydroc20.gms $TITLE A Distillation Problem (hydrocarbon-20 data) $offsymxref offsymlist offuellist offuelxref $inlinecom /* */ /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Model of an n-stage distillation column; the stages i = 0 .. N-1 correspond to a reboiler, N - 2 plates, and a condenser, resp. Reference: More's "A collection of nonlinear model problems", Contributed by Roger Fletcher. Data: hydrocarbon-6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ sets pwr / 0 * 2 /, /* powers of enthalpy eqns*/ P / 1 * 18 /, /* number of plates (N-2) */ C / 1 * 3 /; /* numb of components in mixture */ scalars mbscale / 0.01 /, hbscale / 0.000001 /, K / 9 /, /* feed between plates K/K+1 */ feedt / 100 /, /* temperature of feed */ BB / 40 /, /* amount taken from reboiler */ DD / 60 /, /* amount taken from condenser */ q / 2500000 /, /* heat rate applied to reboiler */ pinvB, pinvT, piB / 1 /, /* pressure in reboiler */ piT / 1 /, /* pressure in condenser */ tBinit / 100 /, tTinit / 100 /, vBinit / 300 /; pinvB = 1 / piB; pinvT = 1 / piT; parameters a(C) / 1 9.647, 2 9.953, 3 9.446 /, b(C) / 1 -2998, 2 -3448.1, 3 -3347.25 /, cc(C) / 1 230.66, 2 235.88, 3 215.31 /, alpha(C,pwr), beta(C,pwr), feedl(C) / 1 30, 2 30, 3 40 /, feedv(C) / 1 0, 2 0, 3 0 /, pinv(P), pi(P) / /* pressure in plates */ 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 1 10 1 11 1 12 1 13 1 14 1 15 1 16 1 17 1 18 1 /, tinit(P) / 1 100 2 100 3 100 4 100 5 100 6 100 7 100 8 100 9 100 10 100 11 100 12 100 13 100 14 100 15 100 16 100 17 100 18 100 /, xBinit(C) / 1 .0 2 .3 3 .1 /, xinit(P,C), xTinit(C) / 1 .5 2 .8 3 .0 /, vinit(P) / 1 300 2 300 3 300 4 300 5 300 6 300 7 300 8 300 9 300 10 300 11 300 12 300 13 300 14 300 15 300 16 300 17 300 18 300 /; pinv(P) = 1 / pi(P); table alpha (C,pwr) 0 1 2 1 0 37.6 0 2 0 48.2 0 3 0 45.4 0 ; table beta (C,pwr) 0 1 2 1 8425 24.2 0 2 9395 35.6 0 3 10466 31.9 0 ; table xinit(P,C) 1 2 3 1 .0 .3 .9 2 .01 .3 .9 3 .2 .4 .8 4 .05 .4 .8 5 .07 .45 .8 6 .09 .5 .7 7 .1 .5 .7 8 .15 .5 .6 9 .2 .5 .6 10 .25 .6 .5 11 .3 .6 .5 12 .35 .6 .5 13 .4 .6 .4 14 .4 .7 .4 15 .42 .7 .4 16 .45 .75 .3 17 .45 .75 .2 18 .5 .8 .1 ; variables tB, /* temperature of reboiler */ t(P), /* temperature of stage P */ tT, /* temperature of condenser */ xB(C), /* proportion of component C leaving reboiler */ x(P,C), /* proportion of component C in l(P) */ xT(C), /* proportion of component C leaving condenser */ vB, /* amount of vapor leaving reboiler */ v(P); /* amount of vapor going from stage P to P+1 */ t.lo(P) = 1; tT.lo = 1; tB.lo = 1; v.lo(P) = 0.01; vB.lo = 0.01; equations psumB, /* proportions sum to 1 */ psum(P), /* "" */ psumT, /* "" */ mbalB(C), /* material balance constraints */ mbal(P,C), /* "" */ mbalT(C), /* "" */ hbalB, /* heat balance constraints */ hbal(P); /* "" */ psumB .. sum(C, pinvB * exp(a(C) + b(C)/(cc(C)+tB) ) * xB(C) ) =e= 1; psum(P) .. sum(C, pinv(P) * exp(a(C) + b(C)/(cc(C) + t(P)) ) * x(P,C) ) =e= 1; psumT .. sum(C, pinvT * exp(a(C) + b(C)/(cc(C)+tT) ) * xT(C) ) =e= 1; mbalB(C) .. (vB + BB) * x("1",C) * mbscale =e= ( vB * pinvB * exp(a(C) + b(C)/(cc(C)+tB) ) * xB(C) + BB * xB(C) ) * mbscale; mbal(P,C) .. ( v(P-1) * pinv(P-1) * exp(a(C) + b(C)/(cc(C)+t(P-1)) ) * x(P-1,C) + (vB * pinvB * exp(a(C) + b(C)/(cc(C)+tB) ) * xB(C))$(ord(P) eq 1) + ( (v(P) + BB)$(ord(P) lt K) + (V(P) - DD)$(ord(P) ge K) ) * (x(P+1,C) + xT(C)$(ord(P) eq card(P))) + feedl(C)$(ord(P) eq K) + feedv(C)$(ord(P) eq K+1) ) * mbscale =e= ( v(P) * pinv(P) * exp(a(C) + b(C)/(cc(C)+t(P)) ) * x(P,C) + ( (v(P-1) + BB)$(ord(P) le K) +(v(P-1) - DD)$(ord(P) gt K) + vB$(ord(P) eq 1) ) * x(P,C) ) * mbscale ; mbalT(C) .. xT(C) =e= sum(P$(ord(P) eq card(P)), pinv(P) * exp(a(C) + b(C)/(cc(C)+t(P)) ) * x(P,C) ); hbalB .. ( (vB + BB) * sum(C, x("1", C) * (alpha(C,"0") + alpha(C,"1")*t("1") + alpha(C,"2")*t("1")*t("1"))) + q ) * hbscale =e= ( vB * pinvB * sum(C, exp(a(C) + b(C)/(cc(C)+tB)) * xB(C) * (beta(C,"0") + beta(C,"1")*tB + beta(C,"2")*tB*tB) ) + BB * sum(C, x("1", C) * (alpha(C,"0") + alpha(C,"1")*tB + alpha(C,"2")*tB*tB)) ) * hbscale ; hbal(P) .. ( v(P-1) * pinv(P-1) * sum(C, exp(a(C) + b(C)/(cc(C)+t(P-1))) * x(P-1,C) * (beta(C,"0") + beta(C,"1")*t(P-1) + beta(C,"2")*t(P-1)*t(P-1) ) ) + (vB * pinvB * sum(C, exp(a(C) + b(C)/(cc(C)+tB)) * xB(C) * (beta(C,"0") + beta(C,"1")*tB + beta(C,"2")*tB*tB) ) )$(ord(P) eq 1) + ( (v(P) + BB)$(ord(P) lt K) + (v(P) - DD)$(ord(P) ge K) ) * sum(C, x(P+1,C) * (alpha(C,"0") + alpha(C,"1")*t(P+1) + alpha(C,"2")*t(P+1)*t(P+1)) + (xT(C) * (alpha(C,"0") + alpha(C,"1")*tT + alpha(C,"2")*tT*tT))$(ord(P) eq card(P)) ) + (sum(C, feedl(C) * (alpha(C,"0") + alpha(C,"1")*feedt + alpha(C,"2")*feedt*feedt)))$(ord(P) eq K) + (sum(C, feedv(C) * (beta(C,"0") + beta(C,"1")*feedt + beta(C,"2")*feedt*feedt)))$(ord(P) eq K+1) ) * hbscale =e= ( v(P) * pinv(P) * sum(C, exp(a(C) + b(C)/(cc(C)+t(P))) * x(P,C) * (beta(C,"0") + beta(C,"1")*t(P) + beta(C,"2")*t(P)*t(P) ) ) + ( (v(P-1) + BB)$(ord(P) le K) +(v(P-1) - DD)$(ord(P) gt K) + vB$(ord(P) eq 1) ) * sum(C, x(P,C) * (alpha(C,"0") + alpha(C,"1")*t(P) +alpha(C,"2")*t(P)*t(P)) ) ) * hbscale ; model hydroc20 / psumB.tB, psum.t, psumT.tT, mbalB.xB, mbal.x, mbalT.xT, hbalB.vB, hbal.v /; option limrow = 0; option limcol = 0; option iterlim = 10000; option reslim = 240; t.l(P) = tinit(P); tT.l = tTinit; tB.l = tBinit; x.l(P,C) = xinit(P,C); xT.l(C) = xTinit(C); xB.l(C) = xBinit(C); v.l(P) = vinit(P); vB.l = vBinit; solve hydroc20 using mcp;