pump.gms : Pump Network Synthesis


The aim is to identify the least costly configuration of centrifugal
pumps that achieves a pre specified pressure rise based on a given
total flowrate.

Small Model of Type : MINLP

Category : GAMS Model library

Main file : pump.gms

$title Pump Network Synthesis (PUMP,SEQ=205)

The aim is to identify the least costly configuration of centrifugal
pumps that achieves a pre specified pressure rise based on a given
total flowrate.

Floudas, C A, Pardalos, P M, Adjiman, C S, Esposito, W R,
Gumus, Z H, Harding, S T, Klepeis, J L, Meyer, C A, and
Schweiger, C A, Handbook of Test Problems in Local and Global
Optimization. Kluwer Academic Publishers, 1999.

Westerlund, T, Petterson, F, and Grossmann, I E, Optimization
of Pump Configuration Problems as a MINLP Problem. Computers
and Chemical Engineering 18, 9 (1994), 845-858.

The entire collection of models can found at

Keywords: mixed integer nonlinear programming, pump network optimization, pump
          configuration, engineering, process optimization

Set i      'set of levels' / 1*3 /;

   wmax    'maximum rotation speed'              / 2950 /
   Vtot    'total volumetric flowrate'           /  350 /
   dPtot   'total pressure rise'                 /  400 /
   nsmax   'maximum number of pumps in series'   /    3 /
   npmax   'maximum number of pumps in parallel' /    3 /;

   Pmax(i) 'maximum power output'
   C(i)    'fixed cost of pump'
   Cd(i)   'operating cost coefficient';

Table ldata(i,*) 'data of the levels'
      Pmax        C    Cd     m1      m2        m3     m4   m5       m6
   1    80  6329.03  1800  19.90  0.1610  0.000561  0.696  629  0.01160
   2    25  2489.31  1800   1.21  0.0644  0.000564  2.950  215  0.11500
   3    45  3270.27  1800   6.52  0.1020  0.000232  0.530  361  0.00946;

Pmax(i) = ldata(i,'Pmax');
C(i)    = ldata(i,'C');
Cd(i)   = ldata(i,'Cd');

   P(i)    'power output of pumps on level i'
   w(i)    'rotation speed for pumps on level i'
   dp(i)   'pressure rise on level i'
   vdot(i) 'flow through pumps on level i'
   x(i)    'fraction of total flow on level i'
   np(i)   'number of parallel lines on level i'
   ns(i)   'number of pumps in series on level i'
   z(i)    'existence of level i'
   objval  'objective function variable';

Positive Variable P, w, dp, vdot, x;
Integer  Variable np, ns;
Binary   Variable z;

   f        'objective function'
   g(i)     'power output calculation for level i'
   gd(i)    'pressure rise calculation for level i'
   sumx     'constraint on volume fractions'
   gvdot(i) 'volume flowrate calculation for pumps on level i'
   gdp(i)   'constraints on pressure rise'
   lw(i)    'logical constraints on w'
   lP(i)    'logical constraints on P'
   ldp(i)   'logical constraints dp'
   lvdot(i) 'logical constraints on vdot'
   lx(i)    'logical constraints on x'
   lnp(i)   'logical constraints on np'
   lns(i)   'logical constraints on ns';

f..          objval  =e= sum(i, (C(i) + Cd(i)*P(i))*np(i)*ns(i)*z(i));

g(i)..       P(i)    =e= ldata(i,'m1')*power(w(i)/wmax,3)
                      +  ldata(i,'m2')*power(w(i)/wmax,2)*vdot(i)
                      -  ldata(i,'m3')*w(i)/wmax*power(vdot(i),2);

gd(i)..      dp(i)   =e= ldata(i,'m4')*w(i)/wmax*vdot(i)
                      +  ldata(i,'m5')*power(w(i)/wmax,2)
                      -  ldata(i,'m6')*power(vdot(i),2);

sumx..       sum(i,x(i)) =e= 1;

gvdot(i)..   x(i) =e= vdot(i)/Vtot*np(i);

gdp(i)..     z(i) =e= dp(i)/dPtot*ns(i);

lw(i)..      w(i)/wmax    =l= z(i);

lP(i)..      P(i)/Pmax(i) =l= z(i);

ldp(i)..     dp(i)/dPtot  =l= z(i);

lvdot(i)..   vdot(i)/Vtot =l= z(i);

lx(i)..      x(i)  =l= z(i);

lnp(i)..     np(i) =l= npmax*z(i);

lns(i)..     ns(i) =l= nsmax*z(i);

P.up(i)    = Pmax(i);
w.up(i)    = wmax;
dp.up(i)   = dPtot;
vdot.up(i) = Vtot;
x.up(i)    = 1;
np.up(i)   = npmax;
ns.up(i)   = nsmax;

Set h 'variable name headers' / P, dp, vdot, w, x, np, ns, z /;

Table gs(i,h) 'global solution'
               P   dp  vdot        w           x  np  ns  z
   1    28.27034  400   160 2855.102  0.91428570   2   1  1
   2     2.63440  200    30 2950.000  0.08571429   1   2  1;

* Initialize starting point
* Turn on all equipment and let model turn some down
* Otherwise NLP solver doesn't find a feasible point
P.l(i)    = P.up(i);
dp.l(i)   = dp.up(i);
vdot.l(i) = vdot.up(i);
w.l(i)    = w.up(i);
x.l(i)    = 0.33;
z.l(i)    = 1;
np.l(i)   = np.up(i);
ns.l(i)   = ns.up(i);

option optCr = 0.0;

Model pump / all /;

solve pump using minlp minimizing objval;

execError = 0;

* Did we find the global solution?
Parameter rep 'solution report';
rep('P',i,'local')    = P.l(i);
rep('dp',i,'local')   = dp.l(i);
rep('vdot',i,'local') = vdot.l(i);
rep('w',i,'local')    = w.l(i);
rep('x',i,'local')    = x.l(i);
rep('z',i,'local')    = z.l(i);
rep('np',i,'local')   = np.l(i);
rep('ns',i,'local')   = ns.l(i);
rep(h,i,'global')     = gs(i,h);
rep(h,i,'diff')       = rep(h,i,'global') - rep(h,i,'local');

option  rep:8:2:1;
display rep;