综合能源系统分析的统一能路理论(三):《稳态与动态潮流计算》(Python代码实现)

 👨‍🎓个人主页:研学社的博客 

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

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

⛳️座右铭:行百里者,半于九十。

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

目录

💥1 概述

📚2 运行结果

2.1 六节点热网动态潮流

2.2  六节点热网稳态潮流

 2.3 七节点气网动态潮流

2.4 七节点气网稳态潮流

🎉3 参考文献

🌈4 Python代码、数据、文章


💥1 概述

摘要:潮流计算作为能源网络分析的基础性应用,在各能源网络已形成成熟但不统一的计算模型与方法。为促进不同能源网络研究的学科融合,该文基于统一能路理论,针对天然气网络与供热网络,提出了相适应的潮流计算方法。在此基础上,补充了基值修正的迭代方法以提高潮流计算的精度,并应用了"边值–初值"等效的方法以在动态潮流计算中隐式地给定初始条件。基于统一能路的潮流计算方法不仅统一了不同能源网络的潮流计算,还统一了同一能源网络内的稳态潮流计算与动态潮流计算,奠定了多异质能流在多时间尺度上统一分析的应用基础。此外,相较以有限差分方法为代表的传统动态潮流计算方法,文中所提出的方法在计算性能上实现了显著的提升。

能源网络的潮流计算,指给定网络的结构、参数与边界条件,确定网络的运行状态,是能源网络

运行、规划的重要基础。根据网络运行状态是否发生变化,潮流计算可分为稳态潮流计算和动态潮流计算,表 1 给出了两者在相关属性上的区别能源网络的潮流计算,指给定网络的结构、参

数与边界条件,确定网络的运行状态,是能源网络运行、规划的重要基础。根据网络运行状态是否发生变化,潮流计算可分为稳态潮流计算和动态潮流计算,表 1 给出了两者在相关属性上的区别

作为能源网络分析的基础性应用,电力网络、天然气网络和供热网络已各自形成相对成熟的潮流计算模型与方法。文献[2]针对电力网络稳态潮流,介绍了高斯–塞德尔方法、牛顿–拉夫逊方法、快速分解法等潮流计算方法。文献 [3] 基 于Weymouth 方程和图论方法,建立了天然气网络的稳态潮流模型,并基于牛顿法提出了一种迭代计算方法,文献[4]将其推广到了多压力等级的天然气网络稳态潮流计算;考虑天然气的慢动态特性,文献[5]基于流体力学方程提出了天然气管道的动态仿真模型,并使用最小平方频谱方法进行了计算, 文献[6]将其从管道推广到网络,建立了天然气网络 的动态潮流模型,并在有限元计算软件 COMSOL Multiphysics 上实现了动态潮流计算。文献[7-8]建立了供热网络水力分析的稳态潮流模型,并应用改进平方根方法进行了稳态潮流计算,而文献[9]建立了供热网络热力分析的稳态潮流模型;考虑热力迁移的长时延效应,文献[10-11]基于热力学方程提出了供热网络热力分析的动态潮流模型,使用不同的有限差分格式进行了求解,并分析了差分稳定性。详细文章见第4部分。

本文包含“综合能源系统分析的统一能路理论(三):潮流计算”中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 参考文献

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

[1]陈彬彬,孙宏斌,吴文传,等.综合能源系统分析的统一能路理论(三): 稳态与动态潮流计算[J].中国电机工程学报, 2020, 40(15):11.

🌈4 Python代码、数据、文章

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值