gastrans.gms : Gas Transmission Problem - Belgium

Description

The problem of distributing gas through a network of pipelines is formulated as
a cost minimization subject to nonlinear flow-pressure relations, material
balances, and pressure bounds. The Belgian gas network is used as an example.

First, we model a straight-forward NLP formulation that can be solved fine
by todays NLP solvers.
Afterwards, the 3-stage approach by Wolf & Smeers is implemented (line 160ff).

de Wolf, D, and Smeers, Y, The Gas Transmission Problem Solved by
and Extension of the Simplex Algorithm. Management Science 46, 11
(2000), 1454-1465.

Keywords: nonlinear programming, discontinuous derivatives, engineering, distribution
          problem, gas transmission, network problem


Small Model of Type : NLP


Category : GAMS Model library


Main file : gastrans.gms

$title Gas Transmission Problem - Belgium (GASTRANS,SEQ=217)

$onText
The problem of distributing gas through a network of pipelines is formulated as
a cost minimization subject to nonlinear flow-pressure relations, material
balances, and pressure bounds. The Belgian gas network is used as an example.

First, we model a straight-forward NLP formulation that can be solved fine
by todays NLP solvers.
Afterwards, the 3-stage approach by Wolf & Smeers is implemented (line 160ff).

de Wolf, D, and Smeers, Y, The Gas Transmission Problem Solved by
and Extension of the Simplex Algorithm. Management Science 46, 11
(2000), 1454-1465.

Keywords: nonlinear programming, discontinuous derivatives, engineering, distribution
          problem, gas transmission, network problem
$offText

$eolCom //

Set
   i 'towns'
     / Anderlues, Antwerpen, Arlon,   Berneau,   Blaregnies
       Brugge,    Dudzele,   Gent,    Liege,     Loenhout
       Mons,      Namur,     Petange, Peronnes,  Sinsin
       Voeren,    Wanze,     Warnand, Zeebrugge, Zomergem   /
   a          'arcs' / 1*24 /
   as(a)      'active arcs'
   ap(a)      'passive arcs'
   aij(a,i,i) 'arc description';

Alias (i,j);

Set nc 'node data column headers'
       / slo 'supply lower bound (mill M3 per day)'
         sup 'supply upper bound (mill M3 per day)'
         plo 'pressure lower bound           (bar)'
         pup 'pressure upper bound           (bar)'
         c   'cost                    ($ per MBTU)' /;

Table Ndata(i,nc) 'node data'
               slo       sup   plo    pup   c
   Anderlues     0     1.2       0   66.2   0
   Antwerpen  -inf    -4.034    30   80.0   0
   Arlon      -inf    -0.222     0   66.2   0
   Berneau       0     0         0   66.2   0
   Blaregnies -inf   -15.616    50   66.2   0
   Brugge     -inf    -3.918    30   80.0   0
   Dudzele       0     8.4       0   77.0   2.28
   Gent       -inf    -5.256    30   80     0
   Liege      -inf    -6.385    30   66.2   0
   Loenhout      0     4.8       0   77.0   2.28
   Mons       -inf    -6.848     0   66.2   0
   Namur      -inf    -2.120     0   66.2   0
   Petange    -inf    -1.919    25   66.2   0
   Peronnes      0     0.96      0   66.2   1.68
   Sinsin        0     0         0   63.0   0
   Voeren     20.344  22.012    50   66.2   1.68
   Wanze         0     0         0   66.2   0
   Warnand       0     0         0   66.2   0
   Zeebrugge   8.870  11.594     0   77.0   2.28
   Zomergem      0     0         0   80.0   0   ;

Set ac 'arc data column headers' / D   'diameter (mm)'
                                   L   'length   (km)'
                                   act 'type indicator (1 = active)' /;

Table AData(a,i,j,*) 'arc data'
                                  D     L  act
    1.Zeebrugge. Dudzele      890.0   4.0
    2.Zeebrugge. Dudzele      890.0   4.0
    3.Dudzele  . Brugge       890.0   6.0
    4.Dudzele  . Brugge       890.0   6.0
    5.Brugge   . Zomergem     890.0  26.0
    6.Loenhout . Antwerpen    590.1  43.0
    7.Antwerpen. Gent         590.1  29.0
    8.Gent     . Zomergem     590.1  19.0
    9.Zomergem . Peronnes     890.0  55.0
   10.Voeren   . Berneau      890.0   5.0    1
   11.Voeren   . Berneau      395.0   5.0    1
   12.Berneau  . Liege        890.0  20.0
   13.Berneau  . Liege        395.0  20.0
   14.Liege    . Warnand      890.0  25.0
   15.Liege    . Warnand      395.0  25.0
   16.Warnand  . Namur        890.0  42.0
   17.Namur    . Anderlues    890.0  40.0
   18.Anderlues. Peronnes     890.0   5.0
   19.Peronnes . Mons         890.0  10.0
   20.Mons     . Blaregnies   890.0  25.0
   21.Warnand  . Wanze        395.5  10.5
   22.Wanze    . Sinsin       315.5  26.0    1
   23.Sinsin   . Arlon        315.5  98.0
   24.Arlon    . Petange      315.5   6.0     ;

Scalar
   T   'gas temperature                (K)' / 281.15  /
   e   'absolute rugosity             (mm)' /   0.05  /
   den 'density of gas relative to air (-)' /   0.616 /
   z   'compressibility factor         (-)' /   0.8   /;

Parameter
   lam(a,i,j)    'lambda constant (page 1464)'
   c2(a,i,j)     'pipe constant   (page 1463)'
   arep(a,i,j,*) 'arc report';

aij(a,i,j) = adata(a,i,j,'L');

as(a) = sum(aij(a,i,j), adata(aij,'act'));
ap(a) = not as(a);

lam(aij(a,i,j)) = 1/sqr(2*log10(3.7*adata(aij,'D')/e));
c2(aij(a,i,j))  = 96.074830e-15*power(adata(aij,'D'),5)/lam(aij)/z/t/adata(aij,'L')/den;
arep(aij,'lam') = lam(aij);
arep(aij,'c2')  = c2(aij);

option  arep:6:3:1;
display arep, as, ap;

Variable
   f(a,i,j) 'arc flow        (1e6 SCM)'
   s(i)     'supply - demand (1e6 SCM)'
   pi(i)    'squared pressure'
   sc       'supply cost';

Equation
   flowbalance(i)   'flow conservation'
   weymouthp(a,i,j) 'flow pressure relationship - passive'
   weymoutha(a,i,j) 'flow pressure relationship - active'
   defsc            'definition of supply cost';

flowbalance(i).. sum(aij(a,i,j), f(aij)) =e= sum(aij(a,j,i), f(aij)) + s(i);

weymouthp(aij(ap,i,j)).. signpower(f(aij),2) =e= c2(aij)*(pi(i) - pi(j));

weymoutha(aij(as,i,j))..       - sqr(f(aij)) =g= c2(aij)*(pi(i) - pi(j));

defsc.. sc =e= sum(i, ndata(i,'c')*s(i));

* supply, demand, pressure, and flow bounds
s.lo(i)  = ndata(i,'slo');
s.up(i)  = ndata(i,'sup');
pi.lo(i) = sqr(ndata(i,'plo'));
pi.up(i) = sqr(ndata(i,'pup'));
f.lo(aij(as,i,j)) = 0;

* initialize flow to avoid getting trapped at signpower(0,2)
f.l(aij) = uniform(-1,1);

Model gastrans / flowbalance, weymouthp, weymoutha, defsc /;

solve gastrans using nlp min sc;

display f.l;

* to run the Wolf & Smeers approach, remove the following $exit
$exit

Parameter
   frep 'flow report'
   sdp  'supply demand and pressure';

frep('NLP',aij,'Flow')  =  f.l(aij);
sdp('NLP',i,'Supply')   =  s.l(i)$(s.l(i) > 0);
sdp('NLP',i,'Demand')   = -s.l(i)$(s.l(i) < 0);
sdp('NLP',i,'Pressure') =  sqrt(pi.l(i));
sdp('NLP','','Obj')     =  sc.l;

Parameter
   flow(a,i,j,*)    'flow bounds'
   pirange(a,i,j,*) 'squared pressure bounds';

flow(aij(ap,i,j),'min') = -sqrt(c2(aij)*(pi.up(j) - pi.lo(i)));
flow(aij(ap,i,j),'max') =  sqrt(c2(aij)*(pi.up(i) - pi.lo(j)));

pirange(aij(ap,i,j),'min') = pi.lo(i) - pi.up(j);
pirange(aij(ap,i,j),'max') = pi.up(i) - pi.lo(j);

option  flow:3:3:1, pirange:3:3:1;
display flow, pirange;

Equation
   defh              'definition of Smeers obj'
   defsig(a,i,j)     'definition of flow direction'
   weymouthp2(a,i,j) 'flow pressure relationship - passive'
   flo
   fup
   pilo
   piup;

Variable
   sig(a,i,j) 'flow direction (-1 or +1 for passive elements)'
   b(a,i,j)   'flow direction ( 1 = i to j)'
   h          'Smeers obj';

Binary Variable b;

weymouthp2(aij(ap,i,j)).. sig(aij)*sqr(f(aij)) =e= c2(aij)*(pi(i) - pi(j));

defh.. h =e= sum(aij(a,i,j), abs(f(aij))*sqr(f(aij))/3/c2(aij));

defsig(aij(ap,i,j)).. sig(aij) =e= 2*b(aij) - 1;

flo(aij(ap,i,j))..  f(aij) =g= flow(aij,'min')*(1 - b(aij));

fup(aij(ap,i,j))..  f(aij) =l= flow(aij,'max')* b(aij);

pilo(aij(ap,i,j)).. pi(i) - pi(j) =g= pirange(aij,'min')*(1 - b(aij));

piup(aij(ap,i,j)).. pi(i) - pi(j) =l= pirange(aij,'max')*b(aij);

* drop the previous solution
pi.l(i)  = 0;
s.l(i)   = 0;
f.l(aij) = uniform(-1,1);

Model
   one   / defh,  flowbalance /
   two   / defsc, flowbalance, weymouthp2, weymoutha /
   three / defsc, flowbalance, weymouthp2, weymoutha, defsig, pilo, piup, flo, fup /;

option limRow = 0, limCol = 0;

solve one using dnlp min h; // provides good initial point

solve two using  nlp min sc;

* assignmenst below fix known solution
* b.fx(*aij)         = 1;
* b.fx(aij('7',i,j)) = 0;
* b.fx(aij('8',i,j)) = 0;

solve three using minlp min sc;

frep('Smeers',aij,'Flow')  =  f.l(aij);
sdp('Smeers',i,'Supply')   =  s.l(i)$(s.l(i) > 0);
sdp('Smeers',i,'Demand')   = -s.l(i)$(s.l(i) < 0);
sdp('Smeers',i,'Pressure') =  sqrt(pi.l(i));
sdp('Smeers','','Obj')     =  sc.l;

option  frep:6:4:1;
display frep, sdp;