在市政給水管網當中，管網平差的目的是在已知節點流量、管段長度的情況下，求得各管段流量和對應的經濟管徑，本科生學習階段了解并掌握管網平差原理及方法是必不可少的環節，
在下面的程式當中，將利用哈代克羅斯法求得管段流量，
程式分為三個.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)
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/2037.html
標籤：其他