管网平差的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
1 ''' 2 dot_set = [0, 1, 2, 3, 4] 3 dot_set_quantity = [-100, 10, 20, 50, 20] 4 path_set = [[0, 1], 5 [0, 2], 6 [1, 3], 7 [2, 3], 8 [2, 4], 9 [3, 4]] 10 diameter_set = [200, 200, 200, 200, 200, 200]#管径列表 11 length_set = [300, 300, 300, 300, 300, 300]#管长列表 12 #生成图 13 graph = {} 14 for k in dot_set: 15 graph[k] = [] 16 print(graph) 17 for k, v in path_set: 18 graph[k].append(v) 19 graph[v].append(k) 20 print(graph) 21 ''' 22 23 24 dot_set = [i for i in range(13)] 25 dot_set_quantity = [-440, 26 50, 40, 50, 30, 20, 27 30, 25, 50, 25, 50, 28 30, 40] 29 path_set = [[0, 11], 30 [0, 12], 31 [1,2], 32 [1,4], 33 [2,3], 34 35 [2,8], 36 [3,11], 37 [4,5], 38 [4,7], 39 [5,6], 40 41 [5,9], 42 [6,10], 43 [7,8], 44 [7,9], 45 [8,11], 46 47 [8,12], 48 [9,10], 49 [10,12] 50 ] 51 diameter_set = [250 for i in range(18)]#管径列表 52 length_set = [400, 200, 400, 200, 200, 53 200, 200, 200, 200, 200, 54 200, 200, 200, 200, 200, 55 400, 200, 200]#管长列表 56 #生成图 57 graph = {} 58 for k in dot_set: 59 graph[k] = [] 60 #print(graph) 61 for k, v in path_set: 62 graph[k].append(v) 63 graph[v].append(k) 64 #print(graph)
pingcha.py
1 import numpy as np 2 import init_input 3 4 class dot(object): 5 def __init__(self, number, quantity_of_flow): 6 self.n = number 7 self.q = quantity_of_flow 8 class path(object): 9 def __init__(self, number, start, finish, quantity_of_flow, diameter, length): 10 self.n = number 11 self.s = start 12 self.f = finish 13 self.q = quantity_of_flow 14 self.d = diameter 15 self.l = length 16 class net(object): 17 def __init__(self, dot_set, path_set): 18 self.dote_set = dot_set 19 self.path_set = path_set 20 21 '''管网模型参数''' 22 dot_set = init_input.dot_set 23 dot_set_quantity = init_input.dot_set_quantity 24 path_set = init_input.path_set 25 diameter_set = init_input.diameter_set 26 length_set = init_input.length_set 27 graph = init_input.graph 28 29 '''创建节点类''' 30 dot_list = [] # dot类列表 31 for i in range(len(dot_set)): 32 dot(i, dot_set_quantity[i]).n = i 33 dot(i, dot_set_quantity[i]).q = dot_set_quantity[i] 34 #print(dot(i, dot_set_quantity[i]).n, dot(i, dot_set_quantity[i]).q) 35 D = dot(i, dot_set_quantity[i]) 36 dot_list.append(D) 37 #print(D) 38 39 '''创建管道类，初始流量为0''' 40 path_list = [] # path类列表 41 for i in range(len(path_set)): 42 path(i, path_set[i][0], path_set[i][1], 0, diameter_set[i], length_set[i]).s = path_set[i][0] 43 path(i, path_set[i][0], path_set[i][1], 0, diameter_set[i], length_set[i]).f = path_set[i][1] 44 #print(i, path_set[i][0], path_set[i][1], 0) 45 L = dot(i, path(i, path_set[i][0], path_set[i][1], 0, diameter_set[i], length_set[i])) 46 path_list.append(L) 47 A = np.zeros((len(dot_set), len(dot_set))) # A是节点连接矩阵 48 for i in path_set: 49 A[i[0]][i[1]] = 1 50 A[i[1]][i[0]] = 1 51 #print(A) 52 53 #**************************************************************************** 54 # 假设某些管段流量，再求初始管段流量 55 v = len(dot_set) 56 e = len(path_set) 57 r = 2 + e - v 58 #print(v, e, r) 59 C = np.zeros((v, e)) 60 for i in range(len(path_set)): 61 m = 1 62 for j in path_set[i]: 63 C[j][i] = m 64 m = -1 65 #print(C) 66 #解方程，得出一组初始流量 67 init_q, useless1, u2, u3 = np.linalg.lstsq(C, dot_set_quantity, rcond=0) 68 #print(init_q) 69 70 # 搜索所有环算法*************************************************************** 71 # 找到所有从start到end的路径 72 def findAllPath(graph, start, end, path=[]): 73 path = path + [start] 74 if start == end: 75 return [path] 76 77 paths = [] # 存储所有路径 78 for node in graph[start]: 79 if node not in path: 80 newpaths = findAllPath(graph, node, end, path) 81 for newpath in newpaths: 82 paths.append(newpath) 83 return paths 84 85 allpath_list = [] 86 for i in range(len(dot_set)): 87 for j in range(i + 1, len(dot_set)): 88 if A[i][j] == 1: 89 allpath = findAllPath(graph, i, j) 90 allpath_list.append(allpath) 91 92 i = True 93 while i == True: 94 if len(allpath_list) != 1: 95 allpath_list[0].extend(allpath_list[1]) 96 del allpath_list[1] 97 else: 98 i = False 99 allpath_list = allpath_list[0] 100 101 for i in range(len(allpath_list) - 1): 102 if len(allpath_list[i]) == 2: 103 allpath_list[i] = [] 104 continue 105 for j in range(i + 1, len(allpath_list)): 106 if sorted(allpath_list[i]) == sorted(allpath_list[j]): 107 allpath_list[j] = [] 108 109 if len(allpath_list[-1]) == 2: 110 allpath_list[-1] = [] 111 m = 0 112 for i in allpath_list: 113 if i == []: 114 m += 1 115 for i in range(m): 116 allpath_list.remove([]) 117 #print(allpath_list)
back_cal.py
1 import time 2 time_1 = time.process_time() 3 4 import numpy as np 5 import pingcha 6 A = pingcha.C 7 #print(A) 8 q = pingcha.init_q 9 print('初始流量:\n',q) 10 Q = pingcha.dot_set_quantity 11 #print(Q) 12 #allpath_list = [[0, 2, 3, 1], [0, 2, 4, 3, 1], [2, 4, 3]] 13 allpath_list = pingcha.allpath_list 14 #allpath_list点转换成线编号,存入pathlist_p，用管段标号表示环路 15 pathlist_p = [] 16 for i in allpath_list: 17 l = [] 18 for j in range(len(i)): 19 a = [i[j-1], i[j]] 20 for n in range(len(pingcha.path_set)): 21 if a[1]==pingcha.path_set[n][0] and a[0]==pingcha.path_set[n][1]: 22 l.append(n) 23 break 24 elif a[0]==pingcha.path_set[n][0] and a[1]==pingcha.path_set[n][1]: 25 l.append(-n) 26 else: 27 None 28 pathlist_p.append(l) 29 #print(pathlist_p)#挑选出的闭合环路 30 31 l = [[0 for i in range(len(pingcha.path_set))] for j in range(len(pathlist_p))] 32 for i in range(len(pathlist_p)): 33 for j in pathlist_p[i]: 34 if j >= 0: 35 l[i][j] = 1 36 else: 37 l[i][-j] = -1 38 L = np.array(l) 39 #print(L) 40 s = [64*pingcha.length_set[i]/(3.1415926**2*100**2*pingcha.diameter_set[i]) for i in range(len(q))] 41 h = [s[i]*(q[i]**2) for i in range(len(q))] 42 #print('h:', h) 43 44 x = 0 45 t = 1#闭合环路的水头误差，当所有环路小于0.01m，即完成平差 46 while t > 0.01: 47 x+=1 48 closure_error_list = []#各环的闭合差组成的列表 49 for a in pathlist_p: 50 closure_error = 0#a环的闭合差 51 sum_sq = 0#环路sq之和 52 for b in a: 53 sum_sq += s[abs(b)]*abs(q[abs(b)]) 54 if b >= 0:#可能会有bug，0号管没法判定方向 55 closure_error += h[abs(b)] 56 else: 57 closure_error -= h[abs(b)] 58 closure_error_list.append(closure_error) 59 rivision_q = closure_error/(2*sum_sq)#校正流量 60 for b in a: 61 if (b>=0 and q[abs(b)]>0) or\ 62 (b<0 and q[abs(b)]<0): 63 q[abs(b)] -= rivision_q 64 elif (b<0 and q[abs(b)]>0) or\ 65 (b>=0 and q[abs(b)]<0): 66 q[abs(b)] += rivision_q 67 68 #根据经济流速选管径 69 t1 = 0 70 while True: 71 t1 += 1 72 if t1 == 4: 73 break 74 v = abs((q[abs(b)]*1000))/(3.1415926*(pingcha.diameter_set[abs(b)]/2)**2)#流速 75 if v<0.6: 76 if pingcha.diameter_set[abs(b)] <= 100: 77 break 78 else: 79 pingcha.diameter_set[abs(b)] -= 50 80 elif v>0.9: 81 if pingcha.diameter_set[abs(b)] >= 400: 82 break 83 else: 84 pingcha.diameter_set[abs(b)] += 50 85 else: 86 #print(pingcha.diameter_set[abs(b)]) 87 #print(v) 88 break 89 #print(pingcha.diameter_set[abs(b)]) 90 #print(v) 91 92 #print(rivision_q) 93 h = [s[i]*(q[i]**2) for i in range(len(q))] 94 t = max([abs(i) for i in closure_error_list]) 95 print('第', x,'次平差') 96 #print('h:', h) 97 print('管段流量:\n', q) 98 print('管径:', [i for i in pingcha.diameter_set]) 99 #print('closure_error_list:', closure_error_list) 100 print('环路最大水头闭合差:', t) 101 102 time_2 = time.process_time() 103 print('耗时:',time_2-time_1)
