随着热电联产系统及其他能源转换设备的应用范围不断扩大，多种能源系统的耦合和联系日益紧密。综合能源系统(integrated energy system，IES)打破了原有的冷/热/电/气供能系统独立设计、规划、运行和控制的模式，在规划、设计、建设和运行的过程中，对各类能源的分配、转化、存储、消费等环节进行有机协调与优化，可提高社会能源供应的安全性、灵活性、可靠性^[1]。然而，传统能源系统在建模与分析时，电/气/热系统一般分开处理，常忽略不同能源系统间的相互作用或做近似处理^[2]。因此如何准确有效地建立综合能源系统的多能流模型成为了一个需要解决的问题。

目前已有诸多学者对热力系统、燃气系统的建模工作进行了研究^[3,4,5]，对于更复杂的电－热－气互联系统，文献[6]与文献[7,8,9]分别采用2种不同的方法来进行多能潮流的计算；前者将系统的全部模型方程联立后实现一次性求解，而后者利用解耦的思想，将系统划分成多个子系统独立求解。无论是哪种方法，都需要采用牛顿拉夫逊法求解非线性方程组。以上文献的侧重点各有不同，但却难以得到一种较为通用的区域综合能源系统潮流的计算方法。

因此，本文采用解耦的思想，根据设备的工作模式总结出4种多能潮流的计算流程，让算法不仅能支持多种耦合设备参与，还可以适用于不同运行模式下的区域综合能源系统。

本文涉及的区域能源系统主要由电力网络、热力网络与天然气网络及不同类型的能源转换设备构成。由于允许对多个能源系统的运行计划进行全面优化，因此该系统比独立的能源系统更具效率且更加环保^[10]，它在能源传输、转换、配送等过程中起到了“承上启下”的关键作用。

区域电力系统主要由热电联产(combined heat and power，CHP)机组、燃气轮机组等设备供电，在并网运行模式下也可由大电网供电；区域热力系统由CHP机组、燃气锅炉、电热锅炉、热泵(heat pump，HP)等设备供热；区域天然气系统由气源供气，主要包括气井与储气站^[11]。表1以能源转换类型为依据对常见的耦合设备进行了分类，图1为典型区域综合能源系统的结构示意图。

区域电力系统的运行特性可由电力系统交流潮流模型来表示。经典的交流潮流模型是一个非线性模型，系统内各节点功率平衡方程如下所示：

(1){Pi=Ui∑j∈iUj(Gijcosθij+Bijsinθij)Qi=Ui∑j∈iUj(Gijsinθij−Bijcosθij)ΔPi=Pis−PiΔQi=Qis−Qi

式中：P_i和Q_i表示节点i的有功、无功功率；U_i为节点i的电压幅值；θ_ij＝θ_i－θ_j为节点i、j的电压相角差；Pis和Qis表示节点的注入有功功率与无功功率；G_ij和B_ij分别是节点导纳矩阵第i行、第j列元素的实部和虚部。

区域热力系统的建模，由水力模型与热工模型两部分组成^[3]。

热力网络的水力模型描述了管网中质流量与管道压降的关系。管网质流量需满足流量的连续性，管道的压降需满足环路压降方程，具体的水力模型如下：

(2){A1HmH=mqBHhf=0hf=KH|mH|mH

式中：A_1H表示热力网络的降阶节点－支路关联矩阵；m_H为热力管道质流量列向量；m_q表示热力网络源节点与负荷节点的注入质流量列向量；B_H表示热力网络的环路－支路关联矩阵；h_f表示热力管道压力损失列向量；K_H表示热力管道的阻力系数列向量，其与管道的长度，内径，粗糙度等硬件参数相关；|m_H|表示向量中每一个元素的绝对值。

热力网络中的每个节点均与3种温度相关联：供水温度T_s，出水温度T_o与回水温度T_r，温度模型主要用来计算负荷节点的供水温度与所有节点的回水温度。负荷节点的热功率通常可由热水流入节点的流量与其所释放的热量共同决定，如式(3)中1式；热水流经管道的热力损失可用式(3)中2式衡量；若有多条管道的热水汇入同一个节点，则称该节点为混合节点，混合节点的温度由从各个管道流入的热水流量与温度共同来决定，可用式(3)中3式表示。

(3){Φ=mqCp(Ts−To)Tend=(Tstart−Ta)e−λLHCpmH+Ta(∑j∈imj)Ti=∑j∈imjTj,end

式中：Φ为节点功率列向量；C_p表示水的比热容；T_s，T_o分别表示节点的供水温度列向量与回水温度列向量；T_start和T_end分别表示某一热水管道始节点水温与末节点水温；λ表示该管道单位长度的总传热系数；L_H表示该管道的长度；T_a为外界环境温度；m_j表示与节点i相连的管道j的质流量；T_j，end表示管道j末端温度；T_i表示热水混合后混合节点i温度。

区域天然气系统模型主要通过流量的连续性建立，即在任何节点的气负荷等于流入流出该节点支路气流量的代数和，可用矩阵形式表示为

(4)A1GQG=LG

式中：A_1G为天然气网络的降阶节点－支路关联矩阵；Q_G为天然气管道流量列向量；L_G为节点载荷列向量。

不同压强等级的天然气网络，其管道压降与气流量存在不同的关系，具体可将管网分为低压，中压与高压3种压强等级，采用不同的稳态流动方程^[11]。对于0～7 500 Pa的低压管网，可用Lacey方程描述：

(5){p1−p2=kq2k=11.7×103LD5

对于7.5×10⁴～7×10⁵ Pa的中压管网，可用polyflo方程描述：

(6){p12−p22=kq1.848k=27.24Lε2D4.848

对于7×10⁵ Pa以上的高压管网，可用Panhandle A方程描述：

(7){p12−p22=kq1.854k=18.43Lε2D4.854

式中：p₁，p₂分别表示管道的始末端的绝对压强；k表示管道的阻力系数；q表示管道的气流量；L，D分别表示管道的长度与内径；ε表示管道的效率因子，反映由于额外的摩擦以及牵制所引起的流量或压强损失^[11]。

区域燃气网络一般选用中压或高压流动方程，因此可将式(6)、式(7)写成如下向量表达：

(8)QG=(ΔpG/KG)1/m≜φ(ΔpG)

式中：Δp_G为管道始末节点压强的平方差；K_G为燃气管道的阻力系数列向量。

综合能源设备模型是建立各个能源网络之间耦合的关键环节，本文涉及到的能源耦合设备模型包括：燃气轮机组，CHP机组，热泵，压缩机，循环水泵等。

1)燃气轮机组。

燃气轮机组是耦合电力网络与燃气网络的关键部分，它在燃气网络中充当燃气负荷，在电网中充当发电设备。燃气轮机组的电出力与所消耗的燃气流量的关系^[4]：

(9)PGT=C1Qgas3+C2Qgas2+C3Qgas

式中：C₁，C₂，C₃根据不同的燃气轮机特征而定，简单起见，取C₁＝C₂＝0，考虑到理想状态下燃气轮机无其他损耗，所消耗的所有燃气均用来发电，则有：

(10)PGT=ηGTQgasqgas

式中：P_GT为机组的发电功率；η_GT为机组的发电效率；Q_gas为燃气耗量；q_gas表示单位体积的天然气的燃烧热值。

2) CHP机组。

CHP机组是一种既通过汽轮机发电生产电能，又通过发电后的蒸汽对热力用户供热的机组^[12]。本文主要考虑CHP机组的2种运行模式：背压模式与抽凝模式。背压运行模式下机组的发电量与发热量成正比，如式(11)中1式；抽凝运行模式下，随着抽气量的增加，机组热输出功率增大而电输出功率减小，如式(11)中2式，3种类型机组的电力输出功率、热力输出功率以及消耗燃气流量的关系如式(11)中3式^[13]。

(11){Pb=Φb/cmPex=−1ZΦex+ηexFinQgasCHP=ΦCHP+PCHPqgasηCHP

式中：P_b，Φ_b分别表示背压式CHP机组的电输出功率与热输出功率；c_m为机组的热电比；P_ex，Φ_ex分别表示抽凝式CHP机组的电输出功率与热输出功率；Z为描述CHP热输出功率与电输出变化比值的参数；η_ex表示全冷凝模式下发电功率的效率；F_in表示燃料输入率，本文中考虑为常数；QgasCHP为CHP机组所消耗的燃气；η_CHP为机组效率。

3)热泵。

热泵的一般模型如下：

(12)ΦHP=CCOPPHP

式中：Φ_HP，P_HP分别为热泵的热输出功率与电输入功率；C_COP为热泵的制冷系数。

4)压缩机。

在天然气网络的稳态模拟中，压缩机的约束条件必须考虑。本文选取“通过压缩机的出口压强为已知量”的假设进行分析。

压缩机根据消耗能源类型的不同可分为电驱动压缩机与燃气驱动压缩机，考虑绝热过程，2种压缩机的模型如下：

(13){PCP=qinαη(α−1)[(po/pin)(α−1)/α−1]QCP=PCP/qgas

式中：P_CP为电压缩机消耗电功率；q_in表示压缩机入口气流量；α为多变指数，一般取1.2～1.4；η表示压缩机的效率，一般取0.75～0.85；p_o，p_in表示压缩机的出口压强与入口压强。

5)循环水泵。

循环水泵需要克服环路的阻力损失，产生供回水管道之间的压差，使热水在系统内可以周而复始地循环。一般模型如下：

(14){Pp=mpgHp106ηpHp=2Δhmax+hc

式中：P_p表示循环水泵所消耗的电功率；m_p表示循环水泵所在节点的管道质流量；H_p表示循环水泵在管网中克服的水头压力；g表示重力加速度；η_p表示循环水泵的效率；Δh_max表示从循环水泵至整个供水网压力最低节点的压头损失；h_c表示热力用户允许的最小压头差。

各个能源网络模型之间的交替迭代求解是一种有效的多能潮流计算方法，由于每个网络模型在数学上均属于非线性方程组，故本文采用牛顿拉夫逊法计算每个模型的近似解。

区域电力网络的潮流求解方法已非常成熟，可用商业软件MATPOWER进行求解，此处不再赘述。

热力网络模型的求解方法大致可分为热力水工联立求解算法与热力水工分立求解算法^[3]。本文采用分立求解算法，具体步骤如下文所述。

1)根据热功率方程式(3)中1式，计算所有流入负荷节点的热水质流量m_q；

2)根据水力模型式(2)可列写牛顿法方程如下，计算管道质流量的修正量并更新m_H；

(15)ΔF(mH)=[ΔmΔh]=[A1HmH−mqBHKH|mH|mH]=0

3)令{T′s=Ts−TaT′r=Tr−Ta，则可根据热工模型式(3)中2，3两式得到关于T_s，T_r的线性方程。

4)求解并更新节点的供回水温度，将新的温度列向量代入步骤1)中进行新一轮的迭代。

具体流程框图如图2所示，其中C_s，C_r分别表示供水网与回水网的系数矩阵，Φ_salck表示热平衡节点的热出力。

区域燃气网络的节点压强与管道压差存在以下关系：

(16)−AGTPG=ΔPG

式中：A_G表示天然气网络的节点－支路关联矩阵；P_G表示燃气网络节点压强平方的列向量。

结合式(4)、式(8)与式(16)可得到燃气网络的牛顿方程如下：

(17)F(pG)=ΔpG=[A1Gφ(−AGTpG)−LG]=0

对于含有压缩机的燃气网络，可对压缩机所在的节点进行断口处理。由图3可知，入口气流量Q₁通过低压侧节点a经压缩机升压后从高压侧节点b流出，此时设出口气流量为Q₂。若压缩机为燃气驱动，其消耗的燃气来自于管道1的气流量，Q₁＝Q₂＋Q_cp；若压缩机为电力驱动，则Q₁＝Q₂。

设与入口节点相连的节点c的压强为p_link，压缩机出入口压强分别为p_inlet，p_outlet，则有如下关系：

(18)pinlet=f(plink,Q2,Qcp)

结合牛顿节点法与压缩机的处理方法，可得到如下区域燃气网络流量计算的流程图，见图4。

本文采用解耦的思想，在每次独立求解完某一系统潮流后更新相应的耦合设备，再根据这些耦合设备，完善下一个系统的源荷数据，从而使得潮流计算可以不断顺利进行，直至满足收敛判据。具体流程如图5所示，其中K_1Max表示最大的总迭代次数，K_2Max表示最大电－热迭代次数。

在使用该算法时，可根据电－热耦合设备在不同情况下的运行模式将多能潮流的计算顺序划分为4种类型，可用n_Flag表示。

当n_Flag＝0时，3种网络之间没有耦合，此时无需考虑多能潮流的计算顺序。

在第2，第3种情况下，综合能源系统电－热耦合设备的工作模式是统一的，即：“以热定电”或“以电定热”。n_Flag＝1时，系统中耦合设备的输出热功率均为已知量，热力系统的数据完备，可优先进行热力系统的流量计算；随后可根据1.2节中的耦合设备模型计算它们的电出力与电负荷，使得电力系统潮流计算得以进行。n_Flag＝2的情况与n_Flag＝1类似，不同之处在于需要先计算数据完善的电力系统潮流再计算热力系统的热工水力模型。在2个上层网络的输出功率均确定的情况下，它们在燃气网络中以载荷的形式完善燃气系统的计算数据，于是在最后可以求解燃气网络的节点压强与管道流量。

在第4种工作模式下，n_Flag＝3，IES中电－热耦合设备2种工作模式同时存在，这种情况在孤岛运行模式下很常见。CHP在满足自身电/热出力约束的同时还要满足系统的模型方程，此时在上层网络中无论先计算电力系统潮流还是热力系统潮流，均在某些节点含有未知的出力或负荷；所以需要给定相应的初值，在迭代计算过程中不断进行修正。

本文区域综合能源系统由IEEE 33节点电力网络^[14]，巴厘岛32节点单环路热力网络，36节点中压天然气网络^[15]组成。热平衡节点为H₁，电平衡节点为E₁，燃气网络的气源节点G₃₆通过固定压缩比为1.5的电压缩机升压至6×10⁵ Pa后，向其他负荷节点供气。为更好地说明能源子系统间的交互过程，假定热力网络不与上级网络相连。设置IES系统能流计算的收敛精度为10^－6，电热耦合与电气耦合的收敛精度为10^－3。系统的具体结构如图6所示。

在算例设置方面，为充分说明综合能源系统的交互过程以及能源耦合设备对多能潮流计算的影响，以下设置了2个算例来进行说明，其能源站的配置见附录表A1—A2部分。

算例1将综合能源系统设置为并网模式，电力网络的电平衡节点与大电网相连，主要为验证所提潮流计算方法的正确性。

由所提多能潮流计算的流程图可知，此算例中热电耦合设备的运行模式均为“以热定电”，说明需采用n_Flag＝2的计算流程，同时由于在能源站1存在电压缩机，说明计算要以电气耦合为收敛判据迭代求解。

选取热力管网中的主干线：H₁-H₂-H₅-H₁₁-H₁₃-H₁₄-H₁₉-H₂₂-H₂₅-H₂₈-H₃(网络中标蓝部分)。图7(a)中质流量柱状图中反映了各个管道热水质流量的大小，其中紫色的折线图表示主干线流量变化，可以看出从两热源主干道流出的热水流量不断减小，最终汇聚于H₂₂处；对比供回水温度模型可发现，H₂₂作为供水网的混合节点，其供水温度在主干线节点中最低；由于回水网与供水网的流向完全相反，H₂₂作为回水网的起点，沿着主干线流向其他节点时具有最高的回水温度。

对燃气网络的分析类似，选取气流量主干道G₃₆-G₈-G₉-G₁₂-G₁₁-G₁₀-G₇-G₃₃-G₁-G₂₆-G₂₉-G₃₀(网络中标蓝部分)，由图7(b)中可看出从管道14到管道13的气流量有较大的跌落，其原因在于位于另一条干线的G₃₄具有较大的气负荷；随着干线气流量的不断减小，节点的压强也不断降低，符合实际运行的规律。

算例2将综合能源系统设置为孤岛模式，并且在算例1的基础上添加了更多的耦合设备与耦合节点。电力系统不与大电网相连，系统中电平衡节点与热平衡节点均通过CHP机组提供能量。同时注意到原先燃气网络中距离气源较远的节点压强偏低，若要求网络中所有节点的压强均在1.5×10⁵ Pa以上，可在G₃₃设置压缩机。

由于热源平衡节点的CHP工作模式为“以热定电”，电平衡节点的CHP工作模式为“以电定热”；由所提算法可知，此场景下需要同时满足热电耦合与电气耦合2种收敛判据迭代求解，因此可以选择抽凝式与背压式2种具有不同电热输出特性的CHP机组，以便能更好地满足电力负荷与热力负荷需求。位于G₃₃处的压缩机选择电力驱动，利用所提方法，将G₃₃处理为一对辅助节点，设置其出口节点G₃₇压强为2.1×10⁵ Pa。

图8给出了算例1、2中综合能源系统中热网的流量与温度结果。对比发现，算例2在孤岛模式下，由于原先由大电网输送至区域电力网络的电功率现在全部由能源站S2承担，为满足所有的电负荷需求，该站的背压式CHP具有较大的电/热输出功率，在热网水力模型中反映在H₃₁具有比热平衡节点H₁更高的管道质流量，如图8(a)所示；在温度模型中反映在最低的干线供水温度由H₂₂偏移至H₁₉，回水温度变化不大，如图8(b)所示。

图9给出了算例1、2中电网潮流的结果对比。算例1，2电负荷的变化主要来自能源站S2，它相对于其他节点的有功负荷相对较小；并且算例1、2中无功负荷未发生变化，所以电网中支路传输功率、节点电压与节点相位变化不大。

图10展示了两算例中气网主干线上气流量与节点压强的分布情况。其中孤岛运行模式下，能源站的电/热出力设备最终消耗的能量均来自于天然气，承担电平衡节点的S2在G₃₄需要更多的燃气负荷，导致气源节点和输送至G₃₄的管道9、12、14的燃气流量明显增加。由于压缩机的加入，使G₃₃的压强有了明显的升高，随后的干线节点压强接近2×10⁵ Pa，符合实际运输要求。

综上，在区域综合能源系统并网模式下，系统间的耦合程度较低，能源耦合设备对各个系统的影响程度相对较低；然而，在孤岛模式下，CHP机组不仅要满足自身的热电输出功率限制，还要满足区域负荷的需求，这会很大程度地影响互联后系统源节点的出力情况；又由于系统中的热能、电能很大程度上依赖于天然气供应，耦合节点的气负荷若始终维持较高的水平，可能会使气网管道长时间重载而造成安全隐患。

本文依据能源耦合设备的运行模式提出了一种多能潮流仿真计算方法。该算法适用于不同的综合能源系统运行场景，并且可以匹配多种不同类型的能源耦合设备。以IEEE 33节点的区域热电力系统、32节点的区域热力系统及36节点的区域天然气系统共同组成的区域综合能源系统为算例，在验证算法有效性的同时，探究了各个能源子系统间的交互过程与运行规律。