gastrans.gms : Gas Transmission Problem - Belgium
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.
Reference:
- 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.
Small Model of Types: MINLP nlp dnlp
$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.
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.
$offtext
$eolcom //
Sets 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
Scalars 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;
Variables f(a,i,j) Arc flow (1e6 SCM)
sig(a,i,j) flow direction (-1 0 +1)
b(a,i,j) flow direction ( 1 = i to j )
s(i) supply - demand (1e6 SCM)
pi(i) squared pressure
sc supply cost
h Smeers obj
Binary Variables b
Equations eq11(i) flow conservation
defsc definition of supply cost
defh definition of Smeers obj
eq12(a,i,j) flow pressure relationship - passive
eq13(a,i,j) flow pressure relationship - active
defsig(a,i,j) definition of flow direction
flo,fup,pilo,piup
parameter flow(a,i,j,*) flow bounds
pirange(a,i,j,*) squared pressure bounds ;
eq11(i).. sum(aij(a,i,j), f(aij)) =e= sum(aij(a,j,i), f(aij)) + s(i);
defsc.. sc =e= sum(i, ndata(i,'c')*s(i));
defh.. h =e= sum(aij(a,i,j), abs(f(aij))*sqr(f(aij))/3/c2(aij));
eq12(aij(ap,i,j)).. sig(aij)*sqr(f(aij)) =e= c2(aij)*(pi(i)-pi(j));
eq13(aij(as,i,j)).. - sqr(f(aij)) =g= c2(aij)*(pi(i)-pi(j));
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);
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;
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;
model one / defh,eq11 /
two / defsc,eq11,eq12,eq13 /
three / defsc,eq11,eq12,eq13,defsig,pilo,piup,flo,fup /;
option limrow=0,limcol=0;
*Initialize f to avoid getting trapped at abs(0)
f.l(aij) = uniform(-1,1);
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;
Parameter frep Flow Report
sdp supply demand and pressure;
frep(aij,'Flow') = f.l(aij);
sdp(i,'Supply') = s.l(i)$(s.l(i) > 0);
sdp(i,'Demand') = -s.l(i)$(s.l(i) < 0);
sdp(i,'Pressure') = sqrt(pi.l(i));
option frep:6:3:1; display frep, sdp;