💥1 概述
本文包含“综合能源系统分析的统一能路理论（三）：潮流计算”中7节点天然气网络和6节点供热网络的数据与Python源码
注意
1. 第三方package要求：pandas、numpy、scipy、matplotlib
2. python版本要求：python 3.x
3. 给出的源码仅对目标算例负责。举例：7节点气网算例中的压气机扬程保持恒定，故未对其进行傅里叶分解；基值修正部分，标准的做法应对计算得到的流量取绝对值再修正（因为基值是非负的），由于算例中给出的正方向恰好为实际正方向，故源码中未取绝对值。如需扩充、修改算例内容，请自行修改代码。
4. 天然气网络的算例中，RT用声速（340m/s）的平方代替。这是一个常数，具体取值不影响算法本身。
📚2 运行结果
2.1 六节点热网动态潮流
2.2 六节点热网稳态潮流
2.3 七节点气网动态潮流
2.4 七节点气网稳态潮流
部分代码：
with context('读取数据与处理'): pipe_table = pd.read_excel('./7节点气网稳态data.xls', sheet_name='Branch') node_table = pd.read_excel('./7节点气网稳态data.xls', sheet_name='Node') numpipe = len(pipe_table) # 支路数 numnode = len(node_table) # 节点数 l = pipe_table['长度(km)'].values * 1e3 # 长度，m d = pipe_table['管径(mm)'].values / 1e3 # 管径，m lam = pipe_table['粗糙度'].values # 摩擦系数 cp = pipe_table['压气机(MPa)'].values * 1e6 # 支路增压，Pa c = 340 # 声速 Apipe = np.pi*d**2/4 # 管道截面积 v = np.ones(numpipe)*5 # 流速基值 with context('基于能路的潮流计算'): MaxIter = 100 # 最大迭代次数 err = [] # 误差记录 vb = [v.copy()] # 基值记录 for itera in range(MaxIter): # 支路参数 Rg = [lam[i]*v[i]/Apipe[i]/d[i] for i in range(numpipe)] Lg = [1/Apipe[i] for i in range(numpipe)] Cg = [Apipe[i]/c**2 for i in range(numpipe)] Ug = [-lam[i]*v[i]**2/2/c**2/d[i] for i in range(numpipe)] # 支路导纳矩阵 Yb1, Yb2, Zb, Ub = [], [], [], [] f = 0 # 稳态，只取零频率分量 for i in range(numpipe): Z, Y = Rg[i], 0 za = np.cosh(np.sqrt(Ug[i]**2+4*Z*Y)/2*l[i]) - Ug[i]/np.sqrt(Ug[i]**2+4*Z*Y)*np.sinh(np.sqrt(Ug[i]**2+4*Z*Y)/2*l[i]) za = za*np.exp(-Ug[i]*l[i]/2) zb = -2*Z/np.sqrt(Ug[i]**2+4*Z*Y)*np.sinh(np.sqrt(Ug[i]**2+4*Z*Y)/2*l[i]) zb = zb*np.exp(-Ug[i]*l[i]/2) zc = -2*Y/np.sqrt(Ug[i]**2+4*Z*Y)*np.sinh(np.sqrt(Ug[i]**2+4*Z*Y)/2*l[i]) zc = zc*np.exp(-Ug[i]*l[i]/2) zd = np.cosh(np.sqrt(Ug[i]**2+4*Z*Y)/2*l[i]) + Ug[i]/np.sqrt(Ug[i]**2+4*Z*Y)*np.sinh(np.sqrt(Ug[i]**2+4*Z*Y)/2*l[i]) zd = zd*np.exp(-Ug[i]*l[i]/2) Yb1.append((za*zd-zb*zc-za)/zb) # 稳态计算中，接地支路不起作用 Yb2.append((1-zd)/zb) # 稳态计算中，接地支路不起作用 Zb.append(-zb) Ub.append(1-za*zd+zb*zc) yb, ub = np.diag(1/np.array(Zb)), np.diag(Ub) # 节点-支路关联矩阵 A = np.zeros([numnode, numpipe]) Ap = np.zeros([numnode, numpipe]) for row in pipe_table.iterrows(): A[int(row[1][1])-1, row[0]] = 1 A[int(row[1][2])-1, row[0]] = -1 Ap[int(row[1][1])-1, row[0]] = 1 Ap[int(row[1][2])-1, row[0]] = 0 # 节点导纳矩阵 Yg_ = np.matmul(np.matmul(A, yb), A.T) - np.matmul(np.matmul(np.matmul(A, yb), ub), Ap.T) # 节点分类 fix_G = node_table[node_table['节点类型']=='定注入'].index.values fix_p = node_table[node_table['节点类型']=='定压力'].index.values Yg_11 = Yg_[fix_G][:,fix_G] Yg_12 = Yg_[fix_G][:,fix_p] Yg_21 = Yg_[fix_p][:,fix_G] Yg_22 = Yg_[fix_p][:,fix_p] assert np.linalg.cond(Yg_11)<1e5 # 确认矩阵不奇异 # 形成广义节点注入向量(给定) 与 节点压力向量(给定) Gn_1 = node_table[node_table['节点类型']=='定注入']['注入 (kg/s)'].values.reshape([-1,1]) # kg/s Gn_1 -= np.matmul(np.matmul(A[fix_G,:], yb), cp.reshape([-1,1])) pn2 = node_table[node_table['节点类型']=='定压力']['气压 (MPa)'].values.reshape([-1,1]) * 1e6 # Pa # 求解零频率网络方程 pn1 = np.matmul(np.linalg.inv(Yg_11), (Gn_1 - np.matmul(Yg_12, pn2))).real Gn_2 = (np.matmul(Yg_21, pn1) + np.matmul(Yg_22, pn2)) Gn_2 += np.matmul(np.matmul(A[fix_p,:], yb), cp.reshape([-1,1])) # 计算失配误差 p = [] pn1, pn2 = pn1.reshape(-1).tolist(), pn2.reshape(-1).tolist() for node in node_table['节点类型'].values: p.append(pn1.pop(0) if node=='定注入' else pn2.pop(0)) p = np.array(p).reshape([-1,1]) I = np.matmul(A.T, p).reshape(-1) + cp.reshape(-1) - (np.array(Ub).reshape([-1,1])*np.matmul(Ap.T, p)).reshape(-1) for i in range(numpipe): I[i] = abs(I[i]*np.diag(yb)[i]/Apipe[i]/(np.matmul(Ap.T, p).reshape(-1)[i])*c**2) err.append(np.linalg.norm(I-v)) print('第%d次迭代，失配误差为%.5f'%(itera+1, err[-1]))
🎉3 参考文献
部分理论来源于网络，如有侵权请联系删除。
