文章编号：1000-3673（2016）10-2942-09    中图分类号：TM 72    文献标志码：A    学科代码：470·40

王英瑞，曾博，郭经，史佳琪，张建华

（新能源电力系统国家重点实验室(华北电力大学)，北京市 昌平区 102206）

Multi-Energy Flow Calculation Method for Integrated Energy System

Containing Electricity, Heat and Gas

WANG Yingrui, ZENG Bo, GUO Jing, SHI Jiaqi, ZHANG Jianhua

(State Key Laboratory of Alternative Electrical Power System With Renewable Energy Sources

(North China Electric Power University), Changping District, Beijing 102206, China)

Abstract: [1]In order to alleviate environmental pollution and improve utilization efficiency of end-used energy, integrated energy system (IES) becomes an important direction of energy structure adjustment in China. Under IES framework, traditional power flow calculation method, with power system as object, is difficult to meet requirements of analyzing coupled energy flow. This paper proposes an extended Newton-Raphson multi-energy flow calculation method, suitable for IES containing electricity, heat and gas. Firstly, models of electricity, heat and natural gas networks in IES are established. In view of complexity of traditional method to solve gas network including compressors, an improved practical method is put forward based on different control modes. On this basis, IES multi-energy flow model is built and Jacobi matrix reflecting multi-energy coupling relationship is derived considering grid-connected and island modes of power system in IES. Validity of the proposed method in multi-energy flow calculation and analysis of interacting characteristics is verified with numerical cases.

Key words: integrated energy system; combined heat and power; multi-energy flow; extended Newton-Raphson method; compressor

摘要：为缓解环境污染，提高终端用能效率，综合能源系统(integrated energy system，IES)已成为我国能源结构调整的重要方向。在IES框架下，传统以电力系统为对象的潮流计算方法将难以满足对互相耦合的能量流分析的需求。提出了一种适用于含电、热、气的IES的扩展Newton-Raphson多能流计算方法。首先，分别建立了IES中电力、热力和天然气网络的数学模型；针对现有方法对含压缩机的天然气网络建模的繁冗性问题，计及不同的控制方式，提出了改进的实用化处理方法；在此基础上，进一步构建了IES多能流计算模型，考虑IES中电力系统并网和孤岛2种运行方式，推导并得出了反映多能源耦合关系的雅可比矩阵。算例验证了该方法在不同运行场景下对多能流进行计算和互动特性分析的有效性。

关键词：综合能源系统；热电联产；多能流；扩展牛顿-拉夫逊法；压缩机

DOI：10.13335/j.1000-3673.pst.2016.10.004

0       引言

随着经济的快速发展，能源、环境问题日益突出，如何实现对能源的清洁高效利用成为近年来人们研究的重点[1]。综合能源系统(integrated energy system，IES)[2-3]集多种能源的生产、输送、分配、转换、存储和消费各环节于一体，能够实现对电、热、气等的综合管理和经济调度，为实现能源的综合利用提供了一种重要解决方案。在IES下，各类能源转换设备如热电联产机组(combined heat and power，CHP)[4]、电热锅炉和燃气锅炉等使电力、热力和天然之间紧密耦合，实现了多能源的互动及转换。如何对多能流(multi-energy flow，MEF)[5]的分布进行有效计算对于指导IES的投资规划和运行决策具有重要意义。

当前，针对MEF相关问题，国内外的研究主要集中在电/气或电/热2种能源网络的联合分析方面。文献[6]将现有的天然气水力计算方法和电力系统潮流计算方法相结合，提出了一种混合潮流的顺序求解法；文献[7-9]建立了能源集线器模型，分别采用集中式优化算法和分布式优化算法对电/气混合最优潮流进行求解；文献[10]从供能可靠性

的角度出发，考虑了天然气的传输延迟和可压缩性，研究了电/气耦合系统的短时优化运行问题；文献[11-12]对冷热电三联供系统进行了建模并研究了协调优化调度和节能问题；文献[13]建立了电/热耦合系统的稳态混合潮流模型，提出了统一计算和顺序计算两种混合潮流的计算方法；文献[14]建立了考虑电力网络和热力网络约束的电热能源集成系统优化模型，利用其研究了风电供热的效益；文献[15]对包含大容量储热的电-热联合系统的研究意义、应用前景和关键技术进行了综述。但上述研究仍未充分考虑MEF内部的耦合性。为弥补以上不足，文献[16-17]提出了以Energy hub为核心的电/热/气混合最优潮流的分解算法，但所用方法不能保证得到全局最优解；文献[18]考虑了时间尺度和网络约束，建立了区域综合能源系统的分层能量管理模型，但对热力部分的考虑仅限于可调节的热负荷。

上述研究表明，IES下的MEF计算已受到国内外学者的广泛关注，但是仍存在以下问题：1）当前研究大多从优化角度侧重于分析多种能源之间的分配关系，而计算模型对电/热/气之间耦合关系及互动机制的考虑仍不够完善；2）现有MEF计算方法对天然气系统中压缩机的处理方法[19]过程繁琐，若将其直接应用于实际大规模系统中，会使计算复杂度大大增加。

鉴于上述问题，本文对含电、热、气的IES的MEF计算方法进行了研究。首先建立了IES中各子系统及耦合环节的数学模型，并提出了一种针对含压缩机的天然气管道的实用化处理流程；在此基础上推导出MEF的计算模型，进而提出了一种扩展的Newton-Raphson电/热/气MEF计算方法，考虑IES中电力系统的运行方式推导出了反映多能源耦合关系的雅可比矩阵，最后给出了MEF计算流程。算例应用验证了本文所提算法在不同应用场景下的有效性。

1  含电、热、气的综合能源系统

含电、热、气的IES由分布于一个区域内的电力系统、热力系统、天然气系统和作为耦合环节的多能源转换设备组成，如图1所示。电力网络将电能从电源侧传输到负荷侧；热力系统由热源、供热网络、回热网络和热负荷构成，高温热水将热量从热源经供热管道输送到热负荷处，通过散热器将热能传输给用户后转变为低温热水，再经回热管道流回热源；天然气系统由气源、供气管道、压缩机和负荷组成，由燃气轮机驱动或电动机驱动的压缩机用来提高供气压力以保证供气量。CHP机组、电锅

图1  含电、热、气的IES

Fig. 1  IES containing electricity, heat and gas

炉和燃气锅炉等设备能够根据需求实现不同能源之间的转换。此外，为提高经济性和运行灵活性，系统中还安装有电储能、热储能和储气设备。在IES框架下对多种能源进行综合调度和协调优化，对满足人类日益增长的用能需求并减少碳排放具有重要意义，因而具有良好普适性和广阔的推广前景，各国在进行能源领域政策的制定时对其重视程度也逐渐提高。

在进行电/热/气MEF计算时，系统的节点分类及对应的已知量和待求量如表1所示。

表1  IES中各节点的分类

Tab. 1  Classification of nodes in IES

系统

变量

已知

未知

电力 系统

电压相角θ

电力平衡节点

PV节点和PQ节点

电压幅值U/kV

PV节点和平衡节点

PQ节点

有功功率P/MW

PV节点和PQ节点

平衡节点

无功功率Q/Mvar

PQ节点

PV节点和PQ节点

热力 系统

水头压力H/m

指定某一节点

除指定节点之外的节点

热功率Φ/MW

除平衡节点之外的节点

平衡节点

供热温度Ts/℃

热源节点

热负荷节点

输出温度To/℃

热负荷节点

—

回热温度Tr/℃

—

所有节点

热水流量m/kg×s-1

—

所有热力管道

天然

气系

统

节点压力p/bar

天然气平衡节点

除平衡节点之外的节点

节点天然气负荷L/m3×h-1

除平衡节点之外的节点

天然气平衡节点

天然气流量f/m3×h-1

—

所有供气管道

2  系统建模

2.1  电力系统模型

与传统的电力系统潮流模型相同，IES中的电力系统用经典的交流潮流模型描述，其节点的功率表达式[20]如下：

                                             (1)

式中：P、Q为节点的有功功率和无功功率；Y为节点导纳矩阵；为节点电压相量。

2.2  热力系统模型

1）水力模型。

热水在网络中流动应满足网络基本定律[13]：各管道的流量在各节点处应满足流量连续性方程，即节点处注入流量等于流出的流量；在一个由管道组成的闭合回路中，水在各管道中流动的压头损失之和为0，即

                                                     (2)

式中：As为供热网络的节点-支路关联矩阵；m为各管道流量；mq为各节点流出的流量；Bh供热管网的回路-支路关联矩阵；hf为压头损失向量，其计算方法为

                                                   (3)

式中K为管道的阻力系数矩阵。

2）热力模型。

对于每一个热负荷节点，供热温度Ts表示热水注入负荷节点之前的温度，输出温度To表示热水流出负荷节点时的温度，回热温度Tr表示从负荷节点流出的热水在管道节点处与其他管道的水混合之后的温度。热力模型[21]为

                                           (4)

                                 (5)

                                     (6)

式(4)为热力网络的节点热功率Φ的表达式，Cp为水的比热容；式(5)表示管道末端温度Tend和始端温度Tstart的关系，Ta为环境温度，l为管道的热传导系数，L为管道的长度；式(6)表示热水在节点处混合前后的温度关系，mout、Tout和min、Tin分别为流出和流入的管道中水的流量和温度。

2.3  天然气系统模型

传统的方法在对含压缩机的天然气系统的状态进行求解时，往往采用对同时含管道模型和压缩机模型的非线性方程组进行整体求解的方法[19]，由于压缩机模型的特殊性，这样会带来计算量大和计算时间过长等问题。为此，本文将含压缩机的管道从网络中分离出来，对不含压缩机的网络和含压缩机的管道2部分分别建模，并提出了一种含压缩机的管道的流量计算方法，通过该方法计算出的管道流量可直接作为其两端网络节点的等效负荷，则在网络计算中不再涉及压缩机的具体数学模型，从而大大简化计算。

1）不含压缩机的天然气网络模型。

该部分模型用于描述天然气网络中节点压力和管道流量所满足的关系。

天然气管道r的稳态流量fr可表示[19,22-23]为

                        (7)

式中：Kr为管道常数；f为函数符号；Δpr2 =pi2 –pj2定义为管道r的压力降；sij用于表征天然气的流动方向，当pi>pj时取+1，否则取-1。

天然气网络的流量连续性方程为

                                                         (8)

式中：Ag为去掉含压缩机管道的天然气网络的节点-支路关联矩阵；f为各管道天然气流量；L为各节点流出的流量。

记Πi = pi2，ΔΠr = Δpr2，则天然气管道的压力降向量ΔΠ可表示为

                                                   (9)

2）含压缩机的管道流量计算方法。

含燃气轮机驱动的压缩机的管道如图2所示，其中，fcom为流过压缩机的流量，fcp为压缩机消耗的天然气量，fmi为压缩机入口管道的流量，fon为出口管道的流量。

图2  含燃气轮机驱动的压缩机的管道

Fig. 2  Pipeline with compressor driven by gas turbine

设压缩机出口处为节点o，入口处为节点i，则含压缩机的管道的数学模型可表示为

                            (10)

式中：kcp为压缩比；Kmi、Kon为入口管道和出口管道的管道常数；pm、pn、pi、po为图2中4个节点的压力；qgas为天然气热值；Tgas为天然气温度；a为多变指数。

压缩机通常有以下4种控制模式：1）出口压力已知；2）入口压力已知；3）通过压缩机的流量已知；4）压缩比已知。本文提出了一种实用化的含压缩机的管道流量计算方法，其基本思路如下：

1）模式1。po已知，利用po、pn及式(10)中第一行计算fon，将其作为fmi的初值，结合pm计算pi和kcp，利用第2行计算fcp，从而得到fmi的新值，反复计算至fmi变化量的绝对值|Δfmi|小于给定值ε。

2）模式2。pi已知，利用式(10)中第四行计算fmi作为fon的初值，从而计算出po及kcp，利用第二行计算fcp并更新fon=fmi-fcp，直至fon变化量的绝对值|Δfon|小于ε。

3）模式3。fcom即fon已知，则可首先利用式(10)的第一行计算出po，从而转化为模式1。

4）模式4。kcp已知，首先给出po的初值，则pi=po/kcp，再用式(10)第1行和第4行分别计算fmi和fon，则fcp=fmi-fon，利用第2行计算fcom，即fon，进一步利用利用第1行计算po的新值，反复计算直至po变化量的绝对值|Δpo|小于ε。计算流程见图3。

图3  含压缩机的管道流量的计算流程

Fig. 3  Flow chart for mass rate calculating of the pipeline with compressor

利用以上计算方法得到的fmi和-fon可分别作为节点m和n的等效负荷feq。

2.4  CHP环节模型

IES中的多能源耦合环节包括CHP机组、电锅炉和燃气锅炉等，为简单起见，本文在此仅考虑含微燃机CHP机组的情况。其电功率PCHP、热功率ΦCHP和耗气量Fin的关系[24]为

                                             (11)

式中：cm为CHP机组的热电比例系数；ηe与所发出的电功率有关，本文假设其为常数。

3  综合能源系统MEF计算方法

3.1  MEF模型

综合上述，构建面向IES的MEF模型

      (12)

式中：第1—2行分别表示电力系统的有功偏差和无功偏差，第3—6行分别表示热力系统的节点热功率偏差、供热网络回路压力降偏差、供热温度偏差和回热温度偏差，第7行表示天然气系统节点流量偏差；PSP、QSP、ΦSP和LSP为系统给定的有功功率、无功功率、热功率和天然气负荷；Asl、Agl分别为供热网络和天然气网络去掉压缩机支路后形成的降阶的关联矩阵；Cs、Cr为分别与供热网络、回热网络的结构和流量有关的矩阵，bs、br为分别与供热温度、输出温度有关的列向量，计算方法见文献[13]；x=[θ, U, m, Ts,load, Tr,,load, Π]T为系统状态量。

3.2  面向MEF的扩展Newton-Raphson算法

本文将Newton-Raphson法推广至MEF问题中，基于扩展Newton-Raphson法的MEF计算方法的迭代公式如下：

                                   (13)

令ΔFe=[ΔP,ΔQ] T、ΔFh=[(ΔΦ,Δp), (ΔTs,ΔTr)] T、ΔFg=Δf分别表示与电、热、气有关的偏差量， xe=[θ, U] T、xh=[m, (Ts,load, Tr,,load)] T、xg=Π分别表示与电、热、气有关的状态量，则J可表示为

  (14)

式中：对角块Jee、Jhh、Jgg分别表示单独的电、热、气系统自身潮流与自身状态量之间的关系，其具体表达式与传统的电力潮流、热力流和天然气流的计算所用的表达式相同；非对角块表示不同能源之间的耦合关系。限于篇幅，正文中仅给出非对角块表达式的推导方法，对角块的表达式见附录A。

IES的天然气网络中，由于平衡节点连接气源，天然气系统内部的状态发生变化时，其供需波动会由平衡节点供气量的变化来承担，不会对热力系统和电力系统产生影响，故式(14)中Jeg、Jhg均为0。

热力系统平衡节点处的热功率由工作于以热定电(following the thermal load，FTL)模式的CHP机组提供，当热力系统中的状态量发生变化时，平衡节点处热功率的波动会同时使该机组所发出的电功率和天然气耗量发生变化，故式(14)中Jeh、Jgh均为非零项。该CHP机组的热功率Φsource,i、电功率Psource,i和燃气耗量Fsource,i表示为

                         (15)

式中Asource,i为热力系统的节点关联矩阵中与该热源有关的行。

针对IES中的电力网络，需考虑2种运行模式：1）并网模式，当电力系统与外部大电网相连时，系统内电功率的波动由大电网来平衡，此时Jhe、Jge均为0；2）孤岛模式，当电力系统工作于孤岛模式时，考虑将电力系统平衡节点设在某一工作于以电定热(FEL)模式的CHP机组处，此时电力系统状态变化所导致的平衡节点处功率的波动会使该CHP的热功率和天然气消耗量均发生变化，故Jhe、Jge均不为0。孤岛模式电力系统平衡节点处的CHP出力和燃气耗量表示式为

                   (16)

根据以上分析，分别将式(15)(16)代入MEF偏差计算式并进行一阶求导，可得2种模式下式(14)的雅可比矩阵非对角块表达式如表2所示。

表2  雅可比矩阵非对角块的表达式

Tab. 2  Expressions for non-diagonal blocks of the Jacobi matrix

对角块

电力系统并网

电力系统孤岛

Jeh

diag{(Ts-To)}Asource,i/cm

diag{(Ts-To)}Asource,i/cm

Jeg

0

0

Jhg

0

0

Jhe

0

cm[Re{jUp Ypk* Uk*}，

Re{-jUp Ypk*(cosθpk-jsinθpk)}]

Jge

0

-cm[Re{jUp Ypk* Uk*}，

Re{-jUp Ypk*(cosθpk-jsinθpk)}]/ηe

Jgh

-diag{(Ts-To)}Asource,i/(cmηe)

-diag{(Ts-To)}Asource,i/(cmηe)

注：Ypk为节点导纳矩阵中的第p行第k列的元素；θpk为节点相角差；“*”表示取共轭。

3.3  计算流程

基于扩展Newton-Raphson的MEF计算的基本思路如下：计算开始时，首先读取IES中各子系统和耦合设备的参数，并输入各状态量的初始值；每次迭代求解之前，都需要根据CHP的工作模式，利用此时获得的温度、流量和电压值计算相应的热功率或电功率以及天然气耗量，为本次迭代计算ΔF做准备；一次迭代完成之后，除更新x之外，还应利用更新后的状态量计算热源的回热温度Tr,source，并计算含压缩机的天然气管道流量，以获得压缩机所在管道两端节点的等效负荷。当|Δx|