MODEL:
! Network equilibrium.                              ( NetEqlx):
 Find equilibrium flows over arcs in a network, given
   demands at each node,
   bounds on pressures at each node,
   resistance characteristics of pipes connecting the nodes,
 Direction of flow over each pipe not known beforehand;
! Keywords: Distribution, Electricity Distribution, Equilibrium,
      Gas Distribution, Network Flow, SignPower function;
! Reference: 
   Hansen, C., K. Madsen, and H. Nielsen(1991), "Optimization of Pipe Networks"
  Math. Prog. vol. 52, no.1, pp. 45-58;
! Best performance is with the settings:
   LINGO  | Options... | Global Solver |  on
   LINGO  | Options... | Global Solver | Multistart Attempts |  off;
SETS:
  NODE: DL, DU, PL, PU, P, DELIVER;  ! P = Pressure at this node;
  ARC( NODE, NODE): R, FLO, ! FLO =  Flow on this arc;
        FLOF, FLOB;         ! Forward or backward;
 ENDSETS
DATA: NODE = A, B, C, D, E, F, G, H; ! Lower & upper limits on demand at each node; DL = 1 2 4 6 8 7 -9999 -9999; DU = 1 2 4 6 8 7 9999 9999; ! Lower & upper limits on pressure at each node; PL = 0 0 0 0 0 0 240 240; PU = 9999 9999 9999 9999 9999 9999 240 240; ! The arcs available and their resistance parameter; ARC = B A, C A, C B, D C, E D, F D, G D, F E, H E, G F, H F; R = 1, 25, 1, 3, 18, 45, 1, 12, 1, 30, 1; PPAM = 1; ! Compressibility parameter; !For incompressible fluids and electricity: PPAM = 1, for gases: PPAM = 2; FPAM = 1.852; !Resistance due to flow parameter; ! electrical networks: FPAM = 1; ! other fluids: 1.8 <= FPAM <= 2; ! For simple optimization networks: FPAM=0, for arcs with flow>=0; ENDDATA @FOR( NODE( K): ! For each node K; ! Bound the pressure; @BND( PL(K), P(K), PU(K)); ! Flow in = amount delivered + flow out; @SUM( ARC( I, K): FLO( I, K)) = DELIVER( K) + @SUM( ARC( K, J): FLO( K, J)); ! Bound on amount delivered at each node; @FREE( DELIVER(k)); @BND( DL(K), DELIVER(K), DU(K)); ); @FOR( ARC( I, J): ! Flow can go either way; @FREE( FLO(I,J)); ! Relate pressures at 2 ends to flow over arc; P(I)^ PPAM - P(J)^ PPAM = ! Alternative 1; ! R(I,J) * @SIGN(FLO(I,J))* @ABS( FLO(I,J))^ FPAM; ! Alternative 2, gets rid of @SIGN( ); ! R(I,J) * FLO(I,J)* @ABS( FLO(I,J))^(FPAM-1); ! Alternative 3, get rid of @ABS( ) and @SIGN( ) but add complementarity constraint; ! R(I,J) * ( FLOF(I,J)^FPAM - FLOB(I,J)^FPAM); ! FLO(I,J) = FLOF(I,J) - FLOB(I,J); ! FLOF(I,J)*FLOB(I,J) = 0; ! Alternative 4, use SIGNPOWER function; R(I,J) * @SIGNPOWER(FLO(I,J), FPAM); ); END