💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

2.1 六节点热网动态潮流

2.2  六节点热网稳态潮流

 2.3 七节点气网动态潮流

2.4 七节点气网稳态潮流

🎉3 参考文献

🌈4 Python代码、数据

---

💥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  # 长度，md = pipe_table['管径(mm)'].values / 1e3  # 管径，mlam = pipe_table['粗糙度'].values  # 摩擦系数cp = pipe_table['压气机(MPa)'].values * 1e6  # 支路增压，Pac = 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], 0za = 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]] = 1A[int(row[1][2])-1, row[0]] = -1Ap[int(row[1][1])-1, row[0]] = 1Ap[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.valuesfix_p = node_table[node_table['节点类型']=='定压力'].index.valuesYg_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/sGn_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))).realGn_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 参考文献

部分理论来源于网络，如有侵权请联系删除。

🌈4 Python代码、数据