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