$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;