【python】EI顶刊复现：综合能源系统分析的统一能路理论(三)：稳态与动态潮流计算程序代码！

适用平台：python 3.8；模块：pandas、numpy、scipy、matplotlib2

程序基于统一能路理论，针对天然气网络和供热网络，借鉴电力系统潮流计算方法，提出了（7节点）气网-（6节点）热网的稳态-动态潮流计算方法，奠定了多能流在对时间尺度上统一分析的基础。程序中算例丰富、注释清晰、干货满满，可扩展性和创新性很高！下面对文章和程序做简要介绍！

程序创新点：

1）借鉴电力系统潮流计算方法，提出了气网-热网的稳态-动态潮流计算方法;

2）应用了“边值–初值”等效的方法在动态潮流计算中隐式地给定初始条件；

3）统一了同一能源网络内的稳态潮流计算与动态潮流计算。

主要工作：

能源网络的潮流计算，即给定网络的结构、参数和边界条件，确定网络的运行状态，是能源网络运行、规划的重要基础。目前，电力网络、供热网络、天然气网络已各自形成相对成熟的潮流计算模型与算法。但是，这些模型与算法相互迥异，在综合能源系统分析中形成了天然的学科壁垒，给相关理论与应用的发展带来了挑战。本文旨在提出针对综合能源系统分析的统一能路理论，基于前两部分中提出的包含气路、水路、热路在内的统一能路模型，针对供热网络和天然气网络的稳态潮流计算和动态潮流计算，发展出基于统一能路的潮流计算方法。补充了基值修正的迭代方法以提高潮流计算的精度，并应用了“边值–初值”等效的方法以在动态潮流计算中隐式地给定初始条件。基于统一能路的潮流计算方法不仅统一了不同能源网络的潮流计算，还统一了同一能源网络内的稳态潮流计算与动态潮流计算，奠定了多异质能流在多时间尺度上统一分析的应用基础。此外，相较以有限差分方法为代表的传统动态潮流计算方法，文中所提出的方法在计算性能上实现了显著的提升。

（1）统一能路理论：

统一能路理论是一种基于能量流、物质流和信息流的能源系统分析方法。该理论将能源系统中的设备、技术和过程看作一个整体，通过分析能量转化、传递和利用过程中的各种因素，实现能源系统的优化设计和调度。在统一能路理论中，电热气耦合系统可看作一个复杂的能源网络。系统中存在电能、热能、燃气等多种能源形式的转化和利用。各种能源之间相互、相互影响，形成一个多层次、多过程的能源流动结构。

（2）稳态-动态潮流区别

（3）气网节点分类

（4）热网节点分类

（5）基于统一能路的动态潮流计算

从数学上说，动态潮流计算的本质是求解一组刻画网络状态变量在时间、空间上分布的时空偏微分-代数方程组。在之前的工作中，已分别将时域中微分方程刻画的天然气网络动态潮流和供热网络动态潮流等价为频域中若干不同频率的正弦稳态的叠加，对于每一个频率下的正弦稳态过程，由一组相量描述的代数方程刻画。因此，基于统一能路的动态潮流计算，归结为若干组复数代数方程的求解；每组复数代数方程求解得到响应相量，通过傅里叶反变换恢复成时域中对应频率的正弦稳态响应；叠加所有正弦响应，即得到原能源网络的时域动态潮流。规程如下图所示：

程序结果：

（1）6节点热网稳态潮流

、

（2）6节点热网动态潮流

（1）7节点气网稳态潮流

（2）7节点气网动态潮流

部分程序：

@contextmanager
def context(event):
    t0 = time.time()
    print('[{}] {} starts ...'.format(time.strftime('%Y-%m-%d %H:%M:%S'), event))
    yield
    print('[{}] {} ends ...'.format(time.strftime('%Y-%m-%d %H:%M:%S'), event))
    print('[{}] {} runs for {:.2f} s'.format(time.strftime('%Y-%m-%d %H:%M:%S'), event, time.time()-t0))


with context('数据读取与处理'):
    tb1 = pd.read_excel('./6节点热网动态data.xls', sheet_name='Node').fillna(0)
    tb2 = pd.read_excel('./6节点热网动态data.xls', sheet_name='Branch')
    tb3 = pd.read_excel('./6节点热网动态data.xls', sheet_name='Device', header=None, index_col=0)
    tb4 = pd.read_excel('./6节点热网动态data.xls', sheet_name='Dynamic')
    # 水力参数
    L = tb2['length'].values * 1e3
    D = tb2['diameter'].values
    lam = tb2['fraction'].values
    npipes, nnodes = len(tb2), len(tb1)
    As = np.array([np.pi*d**2/4 for d in D])
    mb = np.ones(npipes) * 50  # 基值平启动
    rho = 1000
    Ah = np.zeros([nnodes, npipes], dtype=np.int32)
    for i,row in tb2.iterrows():
        Ah[row['from node']-1, i] = 1
        Ah[row['to node']-1, i] = -1
    fix_p = np.where(tb1.type1.values=='定压力')[0]
    fix_G = np.where(tb1.type1.values=='定注入')[0]
    # 热力参数
    c = 4200
    miu = tb2.disspation.values


with context('稳态水力计算'):
    err = []  # 失配误差记录
    mbs = [mb.copy()]  # 流量基值的迭代过程记录
    for itera in range(100):  # 最大迭代次数
        # 更新支路参数
        R = [lam[i]*mb[i]/rho/As[i]**2/D[i]*L[i] for i in range(npipes)]
        E = [-lam[i]*mb[i]**2/2/rho/As[i]**2/D[i]*L[i] for i in range(npipes)]
        # 追加各支路阀、泵的参数
        for i,row in tb2.iterrows():
            if row.pump > 0:
                kp1, kp2, kp3, w = tb3.loc['pump-%d'%int(row.pump),:]
                R[i] += -(2*kp1*mb[i]+kp2*w)
                E[i] += (kp1*mb[i]**2-kp3*w**2)
            if row.valve > 0:
                kv, _, _, _ = tb3.loc['valve-%d'%int(row.valve),:]
                R[i] += 2*kv*mb[i]
                E[i] -= -kv*mb[i]**2
        E = np.array(E).reshape([-1,1])
        yb = np.diag([1/Ri for Ri in R])
        Y = np.matmul(np.matmul(Ah, yb), Ah.T)
        Ygg = Y[fix_G][:,fix_G]
        Ygp = Y[fix_G][:,fix_p]
        Ypg = Y[fix_p][:,fix_G]
        Ypp = Y[fix_p][:,fix_p]
        pp = tb1['pressure(MPa)'].values[fix_p].reshape([1,1]) * 1e6
        G = tb1['injection(kg/s)'].values.reshape([-1,1]) + np.matmul(np.matmul(Ah, yb), E)
        Gg = G[fix_G,:]

欢迎感兴趣的小伙伴关注，小编会不定期更新高质量的学习资料、文章和程序代码，为您的科研加油助力！