管网平差的python程序

在市政给水管网当中，管网平差的目的是在已知节点流量、管段长度的情况下，求得各管段流量和对应的经济管径。本科生学习阶段了解并掌握管网平差原理及方法是必不可少的环节。

在下面的程序当中，将利用哈代克罗斯法求得管段流量。

程序分为三个.py文件：

init_input.py主要是用来输入管网节点编号列表、节点流量列表、管段列表（两端节点列表表示）、管段长度列表。所有的元素应当在列表的对应位置。第一个节点流量为负，表示水厂输入水至管网当中。管段列表当中应当是小节点编号在前，大节点编号在后。

pingcha.py主要是用来对init_input.py文件中的数据进行处理，通过numpy当中的np.linalg.lstsq得出初始流量，准备在back_cal.py当中进行哈代克罗斯法迭代求解管段流量。另外在文件中用算法求出所有环路。算法的原理近乎暴力。首先遍历搜索出所有始、终两点的路径。再判断始、终两点是否直接连接，若是，则成环保留，反之舍去路径。

back_cal.py主要用来哈代克罗斯法迭代，输出结果，并计算程序运行时间。

运用此程序，应先在init_input.py仿照格式填写出管网特征，保存。之后可双击运行back_cal.py，也可在python自带的IDLE下运行。

值得注意的是当有deadlink时，应把支状部分简化到最近的环路节点上。

init_input.py

‘‘‘
dot_set = [0, 1, 2, 3, 4]
dot_set_quantity = [-100, 10, 20, 50, 20]
path_set = [[0, 1],
            [0, 2],
            [1, 3],
            [2, 3],
            [2, 4],
            [3, 4]]
diameter_set = [200, 200, 200, 200, 200, 200]#管径列表
length_set = [300, 300, 300, 300, 300, 300]#管长列表
#生成图
graph = {}
for k in dot_set:
    graph[k] = []
print(graph)
for k, v in path_set:
    graph[k].append(v)
    graph[v].append(k)
print(graph)
‘‘‘


dot_set = [i for i in range(13)]
dot_set_quantity = [-440,
                    50, 40, 50, 30, 20,
                    30, 25, 50, 25, 50,
                    30, 40]
path_set = [[0, 11],
            [0, 12],
            [1,2],
            [1,4],
            [2,3],
            
            [2,8],
            [3,11],
            [4,5],
            [4,7],
            [5,6],
            
            [5,9],
            [6,10],
            [7,8],
            [7,9],
            [8,11],
            
            [8,12],
            [9,10],
            [10,12]
            ]
diameter_set = [250 for i in range(18)]#管径列表
length_set = [400, 200, 400, 200, 200,
              200, 200, 200, 200, 200,
              200, 200, 200, 200, 200,
              400, 200, 200]#管长列表
#生成图
graph = {}
for k in dot_set:
    graph[k] = []
#print(graph)
for k, v in path_set:
    graph[k].append(v)
    graph[v].append(k)
#print(graph)

pingcha.py

import numpy as np
import init_input

class dot(object):
    def __init__(self, number, quantity_of_flow):
        self.n = number
        self.q = quantity_of_flow
class path(object):
    def __init__(self, number, start, finish, quantity_of_flow, diameter, length):
        self.n = number
        self.s = start
        self.f = finish
        self.q = quantity_of_flow
        self.d = diameter
        self.l = length
class net(object):
    def __init__(self, dot_set, path_set):
        self.dote_set = dot_set
        self.path_set = path_set

‘‘‘管网模型参数‘‘‘
dot_set = init_input.dot_set
dot_set_quantity = init_input.dot_set_quantity
path_set = init_input.path_set
diameter_set = init_input.diameter_set
length_set = init_input.length_set
graph = init_input.graph

‘‘‘创建节点类‘‘‘
dot_list = []  # dot类列表
for i in range(len(dot_set)):
    dot(i, dot_set_quantity[i]).n = i
    dot(i, dot_set_quantity[i]).q = dot_set_quantity[i]
    #print(dot(i, dot_set_quantity[i]).n, dot(i, dot_set_quantity[i]).q)
    D = dot(i, dot_set_quantity[i])
    dot_list.append(D)
#print(D)

‘‘‘创建管道类，初始流量为0‘‘‘
path_list = []  # path类列表
for i in range(len(path_set)):
    path(i, path_set[i][0], path_set[i][1], 0, diameter_set[i], length_set[i]).s = path_set[i][0]
    path(i, path_set[i][0], path_set[i][1], 0, diameter_set[i], length_set[i]).f = path_set[i][1]
    #print(i, path_set[i][0], path_set[i][1], 0)
    L = dot(i, path(i, path_set[i][0], path_set[i][1], 0, diameter_set[i], length_set[i]))
    path_list.append(L)
A = np.zeros((len(dot_set), len(dot_set)))  # A是节点连接矩阵
for i in path_set:
    A[i[0]][i[1]] = 1
    A[i[1]][i[0]] = 1
#print(A)

#****************************************************************************
# 假设某些管段流量，再求初始管段流量
v = len(dot_set)
e = len(path_set)
r = 2 + e - v
#print(v, e, r)
C = np.zeros((v, e))
for i in range(len(path_set)):
    m = 1
    for j in path_set[i]:
        C[j][i] = m
        m = -1
#print(C)
#解方程，得出一组初始流量
init_q, useless1, u2, u3 = np.linalg.lstsq(C, dot_set_quantity, rcond=0)
#print(init_q)

# 搜索所有环算法***************************************************************
# 找到所有从start到end的路径
def findAllPath(graph, start, end, path=[]):
    path = path + [start]
    if start == end:
        return [path]

    paths = []  # 存储所有路径
    for node in graph[start]:
        if node not in path:
            newpaths = findAllPath(graph, node, end, path)
            for newpath in newpaths:
                paths.append(newpath)
    return paths

allpath_list = []
for i in range(len(dot_set)):
    for j in range(i + 1, len(dot_set)):
        if A[i][j] == 1:
            allpath = findAllPath(graph, i, j)
            allpath_list.append(allpath)

i = True
while i == True:
    if len(allpath_list) != 1:
        allpath_list[0].extend(allpath_list[1])
        del allpath_list[1]
    else:
        i = False
allpath_list = allpath_list[0]

for i in range(len(allpath_list) - 1):
    if len(allpath_list[i]) == 2:
        allpath_list[i] = []
        continue
    for j in range(i + 1, len(allpath_list)):
        if sorted(allpath_list[i]) == sorted(allpath_list[j]):
            allpath_list[j] = []

if len(allpath_list[-1]) == 2:
    allpath_list[-1] = []
m = 0
for i in allpath_list:
    if i == []:
        m += 1
for i in range(m):
    allpath_list.remove([])
#print(allpath_list)

back_cal.py

import time
time_1 = time.process_time()

import numpy as np
import pingcha
A = pingcha.C
#print(A)
q = pingcha.init_q
print(‘初始流量:\n‘,q)
Q = pingcha.dot_set_quantity
#print(Q)
#allpath_list = [[0, 2, 3, 1], [0, 2, 4, 3, 1], [2, 4, 3]]
allpath_list = pingcha.allpath_list
#allpath_list点转换成线编号,存入pathlist_p，用管段标号表示环路
pathlist_p = []
for i in allpath_list:
    l = []
    for j in range(len(i)):
        a = [i[j-1], i[j]]
        for n in range(len(pingcha.path_set)):
            if a[1]==pingcha.path_set[n][0] and a[0]==pingcha.path_set[n][1]:
                l.append(n)
                break
            elif a[0]==pingcha.path_set[n][0] and a[1]==pingcha.path_set[n][1]:
                l.append(-n)
            else:
                None
    pathlist_p.append(l)
#print(pathlist_p)#挑选出的闭合环路

l = [[0 for i in range(len(pingcha.path_set))] for j in range(len(pathlist_p))]
for i in range(len(pathlist_p)):
    for j in pathlist_p[i]:
        if j >= 0:
            l[i][j] = 1
        else:
            l[i][-j] = -1
L = np.array(l)
#print(L)
s = [64*pingcha.length_set[i]/(3.1415926**2*100**2*pingcha.diameter_set[i]) for i in range(len(q))]
h = [s[i]*(q[i]**2) for i in range(len(q))]
#print(‘h:‘, h)

x = 0
t = 1#闭合环路的水头误差，当所有环路小于0.01m，即完成平差
while t > 0.01:
    x+=1
    closure_error_list = []#各环的闭合差组成的列表
    for a in pathlist_p:
        closure_error = 0#a环的闭合差
        sum_sq = 0#环路sq之和
        for b in a:
            sum_sq += s[abs(b)]*abs(q[abs(b)])
            if b >= 0:#可能会有bug，0号管没法判定方向
                closure_error += h[abs(b)]
            else:
                closure_error -= h[abs(b)]
        closure_error_list.append(closure_error)
        rivision_q = closure_error/(2*sum_sq)#校正流量
        for b in a:
            if (b>=0 and q[abs(b)]>0) or 62                (b<0 and q[abs(b)]<0):
                q[abs(b)] -= rivision_q
            elif (b<0 and q[abs(b)]>0) or 65                  (b>=0 and q[abs(b)]<0):
                q[abs(b)] += rivision_q

            #根据经济流速选管径
            t1 = 0
            while True:
                t1 += 1
                if t1 == 4:
                    break
                v = abs((q[abs(b)]*1000))/(3.1415926*(pingcha.diameter_set[abs(b)]/2)**2)#流速
                if v<0.6:
                    if pingcha.diameter_set[abs(b)] <= 100:
                        break
                    else:
                        pingcha.diameter_set[abs(b)] -= 50
                elif v>0.9:
                    if pingcha.diameter_set[abs(b)] >= 400:
                        break
                    else:
                        pingcha.diameter_set[abs(b)] += 50
                else:
                    #print(pingcha.diameter_set[abs(b)])
                    #print(v)
                    break
                #print(pingcha.diameter_set[abs(b)])
                #print(v)
                    
        #print(rivision_q)
        h = [s[i]*(q[i]**2) for i in range(len(q))]
    t = max([abs(i) for i in closure_error_list])
    print(‘第‘, x,‘次平差‘)
    #print(‘h:‘, h)
    print(‘管段流量:\n‘, q)
    print(‘管径:‘, [i for i in pingcha.diameter_set])
    #print(‘closure_error_list:‘, closure_error_list)
    print(‘环路最大水头闭合差:‘, t)

time_2 = time.process_time()
print(‘耗时:‘,time_2-time_1)