管網平差的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)