Gessler method of solution of the pipe network problem

gessler.py

from numpy import matrix, multiply, power, abs, transpose, sign, sum
from numpy.linalg import inv

L_pipe=matrix([3,2,3,2,3,2,3]) #[km]
D_pipe=matrix([30,20,20,20,20,20,20]) #[cm]
CHW_pipe=matrix([100,100,100,100,100,100,100]) # Friction coefficient of Hazen Williams
_=multiply(power(D_pipe,4.87),power(CHW_pipe,1.852))
r_pipe =1.526e7*L_pipe/_ 

# for each pipe the initial guess
Q0_pipe =matrix([110,30,110,110,80,60,140]) #[m^3/hr]

# for each node, the inputs/outputs
Qd = matrix([-220,0,0,0,20,200]).transpose()  #[m^3/hr]

H1_assumed = 100.0 #[m]
Q_pipe=Q0_pipe
H=0.0
OLD_H=1.0
N_iter = 100
# start the loop
while sum(abs(OLD_H-H))> 0.01:
    # for each pipe, get the 'm'
    M = -1.0/multiply(r_pipe,power(abs(Q_pipe),0.852)) 
     
    # construct the LHS
    LHS=matrix([[-(M[0,0]+M[0,3]),M[0,0],0,M[0,3],0,0],
    [M[0,0],-(M[0,0]+M[0,1]+M[0,4]),M[0,1],0,M[0,4],0],
    [0,M[0,1],-(M[0,1]+M[0,2]+M[0,6]),M[0,2],0,M[0,6]],
    [M[0,3],0,M[0,2],-(M[0,2]+M[0,3]),0,0],
    [0,M[0,4],0,0,-(M[0,4]+M[0,5]),M[0,5]],
    [0,0,M[0,6],0,M[0,5],-(M[0,5]+M[0,6])]])
    
    LHS[0,0]=10**10
    # right hand side vector = outputs
    RHS = -Qd
    
    RHS[0,0]=H1_assumed*10**5
    OLD_H=H
    
    # the solution is to invert the matrix and solve the equation
    H=inv(LHS)*RHS
    
    hloss = matrix([H[0,0]-H[1,0],  H[1,0]-H[2,0],  H[3,0]-H[2,0] , H[0,0]-H[3,0] , H[1,0]-H[4,0],  H[4,0]-H[5,0] , H[2,0]-H[5,0]])
    _ = power(abs(hloss)/r_pipe,1.0/1.852)
    Q_pipe=multiply(sign(abs(hloss)),_)
    
    
print "Heads:\n", H
print ""
print "Discharges: \n", Q_pipe.transpose()