CN108875269B - 一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法 - Google Patents

一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法 Download PDF

Info

Publication number
CN108875269B
CN108875269B CN201810745731.7A CN201810745731A CN108875269B CN 108875269 B CN108875269 B CN 108875269B CN 201810745731 A CN201810745731 A CN 201810745731A CN 108875269 B CN108875269 B CN 108875269B
Authority
CN
China
Prior art keywords
node
power
natural gas
pipeline
delta
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201810745731.7A
Other languages
English (en)
Other versions
CN108875269A (zh
Inventor
赵霞
冯欣欣
胡潇云
杨仑
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing University
Original Assignee
Chongqing University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chongqing University filed Critical Chongqing University
Priority to CN201810745731.7A priority Critical patent/CN108875269B/zh
Publication of CN108875269A publication Critical patent/CN108875269A/zh
Application granted granted Critical
Publication of CN108875269B publication Critical patent/CN108875269B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种考虑电力***多平衡机与天然气***慢动态特性的电‑气互联***多时段能流计算方法，主要步骤为：1)建立考虑多平衡机的电力***潮流模型。2)利用牛顿法求解电力***的潮流分布。3)建立天然气***暂态能流模型。4)利用牛顿法求解天然气***的能流分布。5)判断是否达到总仿真时间Tmax，若未达到总仿真时间，则将天然气***的能流计算结果作为下一时刻能流计算初值，重复步骤2)至步骤4)，若达到总仿真时间，则计算结束并输出电力***的潮流分布和天然气***的能流分布结果。本发明在电气互联***能流计算中同时考虑电力***多平衡机与天然气***慢动态特性，能够得到更准确的***运行状态。

Description

一种考虑电力***多平衡机与天然气***慢动态特性的电- 气互联***多时段能流计算方法
技术领域
本发明涉及电-气互联***能流计算技术领域，具体是一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法。
背景技术
能流计算作为电-气互联***基本运算，是电-气互联***运行和规划的基础。在实际***中，会出现负荷波动、断线故障等情况影响电-气互联***能流。因此有必要对电气互联***能流计算方法展开研究。
现有电-气互联***的能流计算方法中，在计算天然气***能流时，大部分采用天然气***的稳态模型而忽略了该***的慢动态特性。如2016年第184卷《Applied Energy》中“Steady-state analysis of the integrated natural gas and electric powersystem with bi-directional energy conversion”，采用天然气稳态模型计算所得的能流结果可能与实际***的运行结果有一定差异，可能导致操作人员做出不合实际或者次优的决策。另外，在电力***常规潮流计算方法中，线路开断、负荷波动等情况引起***中不平衡功率较大时，默认单个平衡节点吸收全网的不平衡功率是不恰当的，实际***中通常是多台发电机协作共同承担***的不平衡功率。如2017年第37卷第6期《电力自动化设备》中“电-气互联综合能源***多时段暂态能量流仿真”，采用单平衡机来承担***不平衡功率，由此进行潮流计算所得结果可能与***实际运行情况有较大偏差。
综上所述，为更准确地估计电-气互联***负荷波动、断线故障等状况发生后的网络状态，有必要充分考虑电力***多平衡机与天然气***慢动态特性来进行电气互联***能流计算。
发明内容
本发明的目的是解决现有技术中存在的问题。
为实现本发明目的而采用的技术方案是这样的，一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法，主要包括以下步骤：
1)输入电力***和天然气***的基本参数，并设定牛顿法收敛精度ε和总仿真时间Tmax。根据电力***和天然气***的基本参数、牛顿法收敛精度ε和总仿真时间Tmax建立考虑多平衡机的电力***潮流模型。
进一步，所述电力***与天然气***基本参数主要包括电力***网络参数和天然气***参数。
所述电力***参数主要包括电力***支路参数、变压器参数、发电机参数、负荷参数和节点类型。
所述天然气***参数主要包括气源参数、输气管道参数、各节点负荷参数和平衡节点压力初值。
进一步，建立考虑多平衡机的电力***潮流模型的主要步骤如下：
1.1)设定电力***。所述电力***中有Ne个节点。其中，Vδ节点数目为1，PV节点数目为r，PQ节点数目为(Ne-r-1)。电力***中有r+1台发电机组成多平衡机组。多平衡机组共同承担***的不平衡功率。
Vδ节点为的已知量为V和δ，有功功率和无功功率则为待求量。实际进行潮流计算时，总是把平衡节点与电压基准点合选成一个节点。一般选择电力***中主调频电厂的母线作为平衡节点。有时为了提高导纳矩阵算法的收敛性，也可以选择出线数目最多的发电厂母线作为平衡节点或者按其它原则选择平衡节点。
1.2)初始化电力***和天然气***的参数。
电力***中Vδ节点的有功功率初值设定原则为：不考虑功率损耗的情况下，电力***发电机输出有功功率等于负荷有功功率。
PQ节点电压的初值以平启动方法设置，即电压幅值标幺值初值设置为1，相角初值设置为0。
PV节点的电压相角初值设置为0。
天然气***中各节点的压力为天然气***稳态情况下各节点的压力值。天然气***管道分段后新增节点的压力初值设置为管道始末节点压力的等差数列值。管道的流量设置为稳态情况下管道中的流量。
1.3)建立考虑多平衡机的电力***潮流模型。
考虑多平衡机的潮流方程式如下所示：
Figure GDA0001769116170000031
式中，
Figure GDA0001769116170000032
表示电力***中除Vδ节点外其余节点编号所在集合。ζ表示电力***中PQ节点编号所在集合。PGi为节点i处发电机所发有功功率。PDi为节点i处负荷有功功率。Pical为计算得到节点i处传输线上的有功功率。QGj为节点j处发电机所发无功功率。QDj为节点j处负荷无功功率。Qjcal为计算得到节点j处传输线上的无功功率。αi为节点i处发电机的功率分配系数。ΔPΣ为***不平衡功率。
若节点i上无发电机，则PGi＝QGi＝0。
若节点i上无负荷，则PDi＝QDi＝0
ΔPi为电力***节点i处有功功率不平衡量。ΔQi为电力***节点i处无功功率不平衡量。
节点i处传输线上的有功功率Pical和节点i处传输线上的无功功率Qical分别如下所示：
Figure GDA0001769116170000033
式中，
Figure GDA0001769116170000034
为节点i的电压。
Figure GDA0001769116170000035
为节点j处电压的共轭值。
Figure GDA0001769116170000036
为节点i和节点j的互导纳矩阵Yij的共轭值。i、j为电力***节点。
功率分配系数αi如下所示：
Figure GDA0001769116170000037
若节点i处无发电机或发电机不参与承担***不平衡功率，则功率分配系数αi为0。
功率不平衡ΔPΣ如下所示：
Figure GDA0001769116170000038
式中，Ploss为***中的功率损耗。PGi为节点i处发电机所发有功功率。PDi为节点i处负荷有功功率。Ne为电力***节点总数。i为任意节点编号。
***中的功率损耗Ploss如下所示：
Figure GDA0001769116170000041
式中，Gij为节点i与节点j的互电导，δi为节点i处的相角。δj表示节点j处的相角。Ne为电力***节点总数。i，j为任意节点编号。Ui为节点i的电压幅值。
2)利用牛顿法求解电力***的潮流分布。
进一步，利用牛顿法求解电力***潮流的主要步骤如下：
2.1)计算电力***功率不平衡量矩阵ΔF。
电力***功率不平衡量矩阵ΔF如下所示：
Figure GDA0001769116170000042
式中，n＝Ne-1。Ne为电力***节点总数。△P为电力***有功功率不平衡量，△Q为电力***无功功率不平衡量。电力***功率不平衡量矩阵ΔF的前2(n-r)个元素为根据PQ节点计算的功率不平衡量，后r个元素为根据PV节点计算的功率不平衡量。
2.2)计算电力***的雅可比矩阵J。
雅可比矩阵J如下所示：
Figure GDA0001769116170000043
式中，电力***雅可比矩阵维数为(2n-r)×(2n-r)。
雅可比矩阵中的元素Hij、Jij、Nij和Lij分别如公式(8)所示：
Figure GDA0001769116170000044
式中，Hij和Jij分别为有功功率不平衡量和无功功率不平衡量对电压相角不平衡量求偏导所得值，而Nij和Lij分别为有功功率不平衡量和无功功率不平衡量对电压幅值不平衡量求偏导所得值。ΔPi为电力***节点i处有功功率不平衡量。ΔQi为电力***节点i处无功功率不平衡量。δi为节点i处的相角。Uj为节点j的电压幅值。
2.3)计算各节点电压的修正量。
根据功率不平衡量矩阵△F和雅可比矩阵J，计算各节点电压的修正量，从而得到电压修正量矩阵ΔX。所述电压修正量矩阵ΔX中包括电压相角修正量和电压幅值修正量。电压修正量矩阵ΔX如下所示：
ΔX＝J-1ΔF。 (9)
式中，J-1为雅可比矩阵J的逆矩阵。△F为功率不平衡量矩阵。
Figure GDA0001769116170000051
式中，上标T表示转置。n＝Ne-1。Ne为电力***节点总数。功率不平衡量矩阵ΔF的前2(n-r)个元素为根据PQ节点计算的功率不平衡量，后r个元素为根据PV节点计算的功率不平衡量。δi为节点i处的相角。ΔU1为节点i的电压幅值修正量。Δδi为节点i的电压相角修正量。
根据电压修正量更新节点电压向量X：
X(k+1)＝X(k)+ΔX(k)。 (11)
式中，k为迭代次数。当k＝1时，X(k)中各元素为电力***各状态变量初值。
2.4)收敛性判断
当max(|ΔX|)≥ε，则返回步骤2.1)继续迭代计算。
当不平衡量ΔX满足max(|ΔX|)<ε，则计算结束。
3)根据电力***潮流结果和电-气互联***的耦合方式，确定天然气***中天然气的产生量和消耗量，从而建立天然气***暂态能流模型。
进一步，根据电-气耦合设备计算等效天然气负荷量的主要步骤如下：
I)确定电-气互联***的电气耦合方式，将天然气转化为电能。所述电气耦合方式为利用燃气机组进行耦合。
II)设定电力***节点e处接有燃气机组，则燃气机组输出的有功功率Pe和燃气机组所消耗天然气的质量流量GfGPG,e的关系如下所示：
Figure GDA0001769116170000061
式中，ηGPG为燃气机组的效率。LHV为天然气低热值。ρ0为天然气密度。
进一步，建立天然气***暂态能流模型的主要步骤如下：
3.1)设定天然气***。天然气***中有NG个节点，压力已知节点数目为Nt，管道数目为NL，管道分段数为Np，压缩机数目为Nc。
3.2)对天然气***中的天然气管道进行等长度分段。
设定所述天然气***中包括m条管道。每条管道分为Np段。
记第L管道的首末节点的压力分别为p1、p2。记第L管道分段后从首端至末端Np+1个节点的压力分别为：
Figure GDA0001769116170000062
3.3)建立天然气***中燃气轮机驱动压缩机的模型
压缩机消耗的能量BHP如下所示：
Figure GDA0001769116170000063
式中，pout为压缩机出口压力。pin为压缩机入口压力。Gfcom为流经压缩机的流量。B为和天然气温度、压缩机效率、天然气热值有关的常量。z为与天然气热值和天然气压缩因子有关的常量。
当压缩机由天然气驱动时，压缩机所消耗的天然气流量τ如下所示：
Figure GDA0001769116170000064
式中，α、β和γ均为能量转换过程中的效率系数。BHP为压缩机消耗的能量。
压缩机压缩比Rc如下所示：
Figure GDA0001769116170000065
式中，pout为压缩机出口压力。pin为压缩机入口压力。
3.4)建立天然气***暂态能流方程
3.4.1)于任意一条管道，均列写四类非线性方程，则第L条管道的各节点功率不平衡量矩阵ΔF分别如下所示：
Figure GDA0001769116170000071
式中，
Figure GDA0001769116170000072
为暂态常数。
Figure GDA0001769116170000073
表示管道L第l个节点处的压力。
Figure GDA0001769116170000074
为管道L第l个节点处的压力初值，即该节点上一时刻的压力。
Figure GDA0001769116170000075
为管道L第l+1个节点处的流量。
Figure GDA0001769116170000076
为管道L第l-1个节点处的流量。Δt为时间步长。Δl为空间步长。
Figure GDA0001769116170000077
式中，
Figure GDA0001769116170000078
为暂态常数。
Figure GDA0001769116170000079
为管道L第l+1个节点处的压力。
Figure GDA00017691161700000710
为管道L第l+1个节点处的压力。Gfl L为管道L第l个节点处的流量。
Figure GDA00017691161700000711
式中，
Figure GDA00017691161700000712
为暂态常数。
Figure GDA00017691161700000713
为管道L第l+1个节点处的压力。
Figure GDA00017691161700000714
为管道L第l个节点处的压力。Gfl L为管道L第l个节点处的流量。
Figure GDA00017691161700000715
式中，
Figure GDA00017691161700000716
表示管道L第l个节点处的压力。Gfl L为管道L第l个节点处的流量。
Figure GDA00017691161700000717
为暂态常数。
Figure GDA00017691161700000718
为管道L第l+1个节点处的压力。
暂态常数
Figure GDA00017691161700000719
和暂态常数
Figure GDA00017691161700000720
分别如下所示：
Figure GDA00017691161700000721
式中，Z为管道中天然气的压缩系数。R为天然气的气体常数。
Figure GDA00017691161700000722
为天然气管道L的横截面积。T为天然气温度。fc L为管道L的摩擦因子。DL为天然气管道L的内径。
3.4.2)除天然气***压力已知节点之外的其他节点m的流量不平衡量ΔGm如下所示：
Figure GDA0001769116170000081
式中，Gfv→m为节点v流向节点m的流量。Gfm-source为气源注入节点m的流量。Sc为符号函数。当m为压缩机出口节点时，Sc为1。当m为压缩机入口节点时，Sc为-1。Gfm-load为节点m处气负荷消耗的流量。NG为压力天然气***节点个数。Nt为天然气***压力已知节点个数。τ为压缩机所消耗的天然气流量。Gfcom为流经压缩机的天然气流量。
压缩机压力约束方程ΔRc如下：
Figure GDA0001769116170000082
式中，Rc为压缩机c的压缩比。pout,c为压缩机c出口压力。pin,c为压缩机c入口压力。Nc为压缩机个数。
4)利用牛顿法求解天然气***的能流分布。
利用牛顿法求解天然气***能流的主要步骤如下：
4.1)计算天然气***的不平衡量ΔFG。
天然气***的不平衡量ΔFG如下：
Figure GDA0001769116170000083
式中，a∈{1,2,…,NL}。
4.2)计算雅可比矩阵JG。
Figure GDA0001769116170000091
式中，b∈{1,2,…,NG+1}。
雅可比矩阵JG为维数为[2NpNL+(NG-Nt+Nc)]×[2NpNL+(NG-Nt+Nc)]的方阵。
4.3)计算修正量
根据不平衡量△FG和雅克比矩阵JG，计算待求变量的修正量ΔXG。所述待求变量为天然气***能流，主要包括天然气***压力和流量分布。
修正量ΔXG如下所示：
Figure GDA0001769116170000092
式中，JG为雅可比矩阵。ΔFG为天然气***的不平衡量。
Figure GDA0001769116170000093
式中，Δpm为***节点压力修正量。
Figure GDA0001769116170000094
为***新增节点压力修正量。
Figure GDA0001769116170000095
为管道各处流量修正量。ΔGfcom,c为压缩机耗气量修正量。
根据公式(26)，对XG进行更新，即：
XG (k+1)＝XG (k)+ΔXG (k)。 (27)
式中，k表示迭代次数。ΔXG为修正量。XG (k)为天然气***能流。
4.4)收敛性判断
当max(|ΔXG|)≥ε，则返回步骤4.1)继续迭代计算。
当修正量ΔXG满足max(|ΔXG|)<ε，则计算结束。
5)判断是否达到总仿真时间Tmax，若未达到总仿真时间，则将天然气***的能流计算结果作为下一时刻能流计算初值，重复步骤2)至步骤4)，若达到总仿真时间，则计算结束并输出电力***的潮流分布和天然气***的能流分布结果。
本发明的技术效果是毋庸置疑的。本发明针对现有电-气互联***能流计算方法的不足，充分考虑电力***多平衡机与天然气***慢动态特性，提供一种电-气互联***多时段能流计算方法。在负荷波动、断线故障等造成***较大不平衡功率的情况下，本发明在电气互联***能流计算中同时考虑电力***多平衡机与天然气***慢动态特性，能够得到更准确的***运行状态。
附图说明
图1为方法流程框图；
图2为管道分段示意图；
图3为算例***结构图；
图4为两种方法计算电-气互联***能流所得到的天然气***节点5压力仿真结果图。
具体实施方式
下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。
实施例1：
一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法，主要包括以下步骤：
1)输入电力***和天然气***的基本参数，并设定牛顿法收敛精度ε和总仿真时间Tmax。根据电力***和天然气***的基本参数、牛顿法收敛精度ε和总仿真时间Tmax建立考虑多平衡机的电力***潮流模型。
进一步，所述电力***与天然气***基本参数主要包括电力***网络参数和天然气***参数。
所述电力***参数主要包括电力***支路参数、变压器参数、发电机参数、负荷参数和节点类型。
所述天然气***参数主要包括气源参数、输气管道参数、各节点负荷参数和平衡节点压力初值。
进一步，建立考虑多平衡机的电力***潮流模型的主要步骤如下：
1.1)设定电力***。所述电力***中有Ne个节点。其中，Vδ节点数目为1，PV节点数目为r，PQ节点数目为(Ne-r-1)。电力***中有r+1台发电机组成多平衡机组。多平衡机组共同承担***的不平衡功率。
1.2)初始化电力***和天然气***的参数。
电力***中Vδ节点的有功功率初值设定原则为：不考虑功率损耗的情况下，电力***发电机输出有功功率等于负荷有功功率。
PQ节点电压的初值以平启动方法设置，即电压幅值标幺值初值设置为1，相角初值设置为0。
PV节点的电压相角初值设置为0。
天然气***中各节点的压力为天然气***稳态情况下各节点的压力值。天然气***管道分段后新增节点的压力初值设置为管道始末节点压力的等差数列值。管道的流量设置为稳态情况下管道中的流量。
1.3)建立考虑多平衡机的电力***潮流模型。
考虑多平衡机的潮流方程式如下所示：
Figure GDA0001769116170000111
式中，
Figure GDA0001769116170000112
表示电力***中除Vδ节点外其余节点编号所在集合。ζ表示电力***中PQ节点编号所在集合。PGi为节点i处发电机所发有功功率。PDi为节点i处负荷有功功率。Pical为计算得到节点i处传输线上的有功功率。QGj为节点j处发电机所发无功功率。QDj为节点j处负荷无功功率。Qjcal为计算得到节点j处传输线上的无功功率。αi为节点i处发电机的功率分配系数。ΔPΣ为功率不平衡量。
若节点i上无发电机，则PGi＝QGi＝0。
若节点i上无负荷，则PDi＝QDi＝0
ΔPi为电力***节点i处有功功率不平衡量。ΔQi为电力***节点i处无功功率不平衡量。
节点i处传输线上的有功功率Pical和节点i处传输线上的无功功率Qical分别如下所示：
Figure GDA0001769116170000121
式中，
Figure GDA0001769116170000122
为节点i的电压。
Figure GDA0001769116170000123
为节点j处电压的共轭值。
Figure GDA0001769116170000124
为节点i和节点j的互导纳矩阵Yij的共轭值。i、j为电力***节点。
功率分配系数αi如下所示：
Figure GDA0001769116170000125
若节点i处无发电机或发电机不参与承担***不平衡功率，则功率分配系数αi为0。
***不平衡功率ΔPΣ如下所示：
Figure GDA0001769116170000126
式中，Ploss为***中的功率损耗。PGi为节点i处发电机所发有功功率。PDi为节点i处负荷有功功率。Ne为电力***节点总数。i为任意节点编号。
***中的功率损耗Ploss如下所示：
Figure GDA0001769116170000127
式中，Gij为节点i与节点j的互电导，δi为节点i处的相角。δj表示节点j处的相角。Ne为电力***节点总数。i，j为任意节点编号。Ui为节点i的电压幅值。
2)利用牛顿法求解电力***的潮流分布。
进一步，利用牛顿法求解电力***潮流的主要步骤如下：
2.1)计算电力***功率不平衡量矩阵ΔF。
电力***功率不平衡量矩阵ΔF如下所示：
Figure GDA0001769116170000128
式中，n＝Ne-1。Ne为电力***节点总数。△P为电力***有功功率不平衡量，△Q为电力***无功功率不平衡量。电力***功率不平衡量矩阵ΔF的前2(n-r)个元素为根据PQ节点计算的功率不平衡量，后r个元素为根据PV节点计算的功率不平衡量。
2.2)计算电力***的雅可比矩阵J。
雅可比矩阵J如下所示：
Figure GDA0001769116170000131
式中，电力***雅可比矩阵维数为(2n-r)×(2n-r)。
雅可比矩阵中的元素Hij、Jij、Nij和Lij分别如公式(8)所示：
Figure GDA0001769116170000132
式中，Hij为有功功率不平衡量对电压相角不平衡量求偏导所得值。
Jij为无功功率不平衡量对电压相角不平衡量求偏导所得值。
Nij,为有功功率不平衡量对电压幅值不平衡量求偏导所得值。
Lij为无功功率不平衡量对电压幅值不平衡量求偏导所得值。
ΔPi为电力***节点i处有功功率不平衡量。ΔQi为电力***节点i处无功功率不平衡量。δi为节点i处的相角。Uj为节点j的电压幅值。
2.3)计算各节点电压的修正量。
根据功率不平衡量矩阵△F和雅可比矩阵J，计算各节点电压的修正量，从而得到电压修正量矩阵ΔX。所述电压修正量矩阵ΔX中包括电压相角修正量和电压幅值修正量。电压修正量矩阵ΔX如下所示：
ΔX＝J-1ΔF。 (9)
式中，J-1为雅可比矩阵J的逆矩阵。△F为功率不平衡量矩阵。
Figure GDA0001769116170000141
式中，上标T表示转置。n＝Ne-1。Ne为电力***节点总数。功率不平衡量矩阵ΔF的前2(n-r)个元素为根据PQ节点计算的功率不平衡量，后r个元素为根据PV节点计算的功率不平衡量。δi为节点i处的相角。ΔU1为节点i的电压幅值修正量。Δδi为节点i的电压相角修正量。
根据电压修正量更新节点电压向量X：
X(k+1)＝X(k)+ΔX(k)。 (11)
式中，k为迭代次数。当k＝1时，X(k)中各元素为电力***各状态变量初值。
2.4)收敛性判断
当max(|ΔX|)≥ε，则返回步骤2.1)继续迭代计算。
当不平衡量ΔX满足max(|ΔX|)<ε，则计算结束。
3)根据电力***潮流结果和电-气互联***的耦合方式，确定天然气***中天然气的产生量和消耗量，从而建立天然气***暂态能流模型。
进一步，根据电-气耦合设备计算等效天然气负荷量的主要步骤如下：
I)确定电-气互联***的电气耦合方式，将天然气转化为电能。所述电气耦合方式为利用燃气机组进行耦合。
II)设定电力***节点e处接有燃气机组，则燃气机组输出的有功功率Pe和燃气机组所消耗天然气的质量流量GfGPG,e的关系如下所示：
Figure GDA0001769116170000142
式中，ηGPG为燃气机组的效率。LHV为天然气低热值。ρ0为天然气密度。
进一步，建立天然气***暂态能流模型的主要步骤如下：
3.1)设定天然气***。天然气***中有NG个节点，压力已知节点数目为Nt，管道数目为NL，管道分段数为Np，压缩机数目为Nc。
3.2)对天然气***中的天然气管道进行分段。
设定所述天然气***中包括m条管道。每条管道分为Np段。
记第L管道的首末节点的压力分别为p1、p2。记第L管道分段后从首端至末端Np+1个节点的压力分别为：
Figure GDA0001769116170000151
3.3)建立天然气***中燃气轮机驱动压缩机的模型
压缩机消耗的能量BHP如下所示：
Figure GDA0001769116170000152
式中，pout为压缩机出口压力。pin为压缩机入口压力。Gfcom为流经压缩机的流量。B为和天然气温度、压缩机效率、天然气热值有关的常量。z为与天然气热值和天然气压缩因子有关的常量。
当压缩机由天然气驱动时，压缩机所消耗的天然气流量τ如下所示：
Figure GDA0001769116170000153
式中，α、β和γ均为能量转换过程中的效率系数。BHP为压缩机消耗的能量。
压缩机压缩比Rc如下所示：
Figure GDA0001769116170000154
式中，pout为压缩机出口压力。pin为压缩机入口压力。
3.4)建立天然气***暂态能流方程
3.4.1)对于任意一条管道，均列写四类非线性方程，进而需要计算ΔF1 L、
Figure GDA0001769116170000155
四组不平衡量。第L条管道的各节点功率不平衡量矩阵ΔF分别如下所示：
Figure GDA0001769116170000156
式中，
Figure GDA0001769116170000157
为暂态常数。
Figure GDA0001769116170000158
表示管道L第l个节点处的压力。
Figure GDA0001769116170000159
为管道L第l个节点处的压力初值，即该节点上一时刻的压力。
Figure GDA00017691161700001510
为管道L第l+1个节点处的流量。
Figure GDA00017691161700001511
为管道L第l-1个节点处的流量。Δt为时间步长。Δl为空间步长。
Figure GDA0001769116170000161
式中，
Figure GDA0001769116170000162
为暂态常数。
Figure GDA0001769116170000163
为管道L第l+1个节点处的压力。
Figure GDA0001769116170000164
为管道L第l+1个节点处的压力。Gfl L为管道L第l个节点处的流量。
Figure GDA0001769116170000165
式中，
Figure GDA0001769116170000166
为暂态常数。
Figure GDA0001769116170000167
为管道L第l+1个节点处的压力。
Figure GDA0001769116170000168
为管道L第l个节点处的压力。Gfl L为管道L第l个节点处的流量。
Figure GDA0001769116170000169
式中，
Figure GDA00017691161700001610
表示管道L第l个节点处的压力。Gfl L为管道L第l个节点处的流量。
Figure GDA00017691161700001611
为暂态常数。
Figure GDA00017691161700001612
为管道L第l+1个节点处的压力。
暂态常数
Figure GDA00017691161700001613
和暂态常数
Figure GDA00017691161700001614
分别如下所示：
Figure GDA00017691161700001615
式中，Z为管道中天然气的压缩系数。R为天然气的气体常数。
Figure GDA00017691161700001616
为天然气管道L的横截面积。T为天然气温度。fc L为管道L的摩擦因子。DL为天然气管道L的内径。
3.4.2)除天然气***压力已知节点之外的其他节点m的流量不平衡量ΔGm如下所示：
Figure GDA00017691161700001617
式中，Gfv→m为节点v流向节点m的流量。Gfm-source为气源注入节点m的流量。Sc为符号函数。当m为压缩机出口节点时，Sc为1。当m为压缩机入口节点时，Sc为-1。Gfm-load为节点m处气负荷消耗的流量。NG为压力天然气***节点个数。Nt为天然气***压力已知节点个数。τ为压缩机所消耗的天然气流量。Gfcom为流经压缩机的天然气流量。
压缩机压力约束方程ΔRc如下：
Figure GDA0001769116170000171
式中，Rc为压缩机c的压缩比。pout,c为压缩机c出口压力。pin,c为压缩机c入口压力。Nc为压缩机个数。
4)利用牛顿法求解天然气***的能流分布。
利用牛顿法求解天然气***能流的主要步骤如下：
4.1)计算天然气***的不平衡量ΔFG。
天然气***的不平衡量ΔFG如下：
Figure GDA0001769116170000172
式中，a∈{1,2,…,NL}。
4.2)计算雅可比矩阵JG。
Figure GDA0001769116170000173
式中，b∈{1,2,…,NG+1}。
雅可比矩阵JG为维数为[2NpNL+(NG-Nt+Nc)]×[2NpNL+(NG-Nt+Nc)]的方阵。
4.3)计算修正量
根据不平衡量△FG和雅克比矩阵JG，计算待求变量的修正量ΔXG。所述待求变量为天然气***能流，主要包括天然气***各处的压力和流量分布。
修正量ΔXG如下所示：
Figure GDA0001769116170000181
式中，JG为雅可比矩阵。ΔFG为天然气***的不平衡量。
Figure GDA0001769116170000182
式中，Δpm为***节点压力修正量。
Figure GDA0001769116170000183
为***新增节点压力修正量。
Figure GDA0001769116170000184
为管道各处流量修正量。ΔGfcom,c为压缩机耗气量修正量。
根据公式(26)，对XG进行更新，即：
Figure GDA0001769116170000185
式中，k表示迭代次数。ΔXG为修正量。XG (k)为天然气***能流。
4.4)收敛性判断
当max(|ΔXG|)≥ε，则返回步骤4.1)继续迭代计算。
当修正量ΔXG满足max(|ΔXG|)<ε，则计算结束。
5)判断是否达到总仿真时间Tmax，若未达到总仿真时间，则将天然气***的能流计算结果作为下一时刻能流计算初值，重复步骤2)至步骤4)，若达到总仿真时间，则计算结束并输出电力***的潮流分布和天然气***的能流分布结果。
实施例2：
一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法的仿真实验，主要包括以下步骤：
1)输入基础数据及初始化
1.1)输入基础数据
本实施例的电-气互联***为IEEE14节点电力***与6节点天然气***联合而成的电-气互联***
输入电力***网络参数，包括：支路参数、变压器参数、发电机参数、负荷参数、节点类型。天然气***参数：气源参数、输气管道参数、各节点负荷参数、平衡节点压力初值。设置牛顿法收敛精度。电力***网络参数采用IEEE14节点电力网络数据，天然气***网络参数采用2017年第6卷32期《IEEE TRANSACTIONS ON POWER SYSTEMS》中“An EquivalentModel of Gas Networks for Dynamic Analysis of Gas-Electricity Systems”一文中的6节点天然气***的数据。
1.2)参数初始化
电力***中，节点1为Vδ节点，节点2、3、6、8为PV节点，其余节点为PQ节点，PQ节点的电压的初值采用平启动的方法选取，即电压幅值标幺值为1，相角为0.PV节点的电压相角设置为0。
Vδ节点的有功功率初值按照“不考虑功率损耗的情况下，***发电机输出有功功率等于负荷有功功率”的原则来取，在本算例中计算得到输出有功为219MW。
天然气***中各节点压力的初值设置为算例***稳态条件下各节点压力，节点1到6的压力分别为p1＝4MPa,p2＝4.13MPa,p3＝4.04MPa,p4＝2.80MPa,p5＝3.58MPa,p6＝3.06MPa,管道分段后产生的新增节点压力设置为管道首末端压力的等差数列值，各支路流量设置为稳态条件下各管道的流量，管道1到6的流量分别为Gf1＝-4.6,Gf2＝17.6,Gf3＝10.6,Gf4＝34.4,Gf5＝-8.4,Gf6＝26.4,单位为kg/s.新增节点处的流量初值取为与所处管道流量一致。
2)建立考虑多平衡机的潮流模型。
初始情况下该***的功率损耗：
Figure GDA0001769116170000191
计算***中的不平衡功率：
Figure GDA0001769116170000192
设置节点1处发电机承担不平衡功率的份额为60％，节点3发电机承担30％，节点6发电机承担10％，节点2与节点8所承担不平衡功率的份额为0。根据式(1)列写考虑多平衡机的潮流方程，如节点3的潮流方程可列写为下式：
Figure GDA0001769116170000193
3)牛顿法求解考虑多平衡机的电力***潮流
3.1)计算功率不平衡量
通过上面列写的潮流方程，将电力***网络参数以及各状态变量取值代入潮流方程，可以计算得到功率不平衡量ΔF：
Figure GDA0001769116170000201
3.2)计算雅可比矩阵
计算得到雅可比矩阵如下：
Figure GDA0001769116170000202
3.3)计算电压修正量
根据上述步骤获得的不平衡量△F和雅克比矩阵J，以第一次迭代后的结果为例，计算出各节点电压的修正量：
Figure GDA0001769116170000203
根据上述结果更新节点电压向量X：
Figure GDA0001769116170000204
3.4)收敛性判断
当修正量ΔX满足max(|ΔX|)<ε，则计算结束。
当max(|ΔX|)≥ε，则返回步骤3.1继续迭代计算。
根据前面步骤，计算可知当迭代4次后，满足收敛条件，此时max(|△X|)＝9.77E-10＜ε。
3.5)计算发电机功率及***功率损耗
根据公式(1)计算得到各发电机节点所发出的功率，根据式5可计算得到***的功率损耗。计算得到节点3与节点6处发电机所发出的有功功率为：PG3＝3.83MW。PG4＝1.28MW。
4)电气耦合设备计算天然气产生量与消耗量
本算例中采用燃气机组作为电气耦合设备，算例***中，电力***节点3与节点6处发电机为燃气机组，分别于天然气***节点4与节点5相连。燃气机组的效率ηGPG取80％，天然气低热值LHV取37.26MJ/m3，天然气密度ρ0取0.7174kg/m3，则根据式(12)计算得到电力***节点3与节点6处发电机所需天然气流量为：GfGPG3＝0.0921kg/s。GfGPG6＝0.0307kg/s。
天然气***中节点5与节点6处原负荷加上燃气机组所需的天然气量为电气互联***中天然气***节点5与节点6处的实际负荷。
5)天然气***暂态能流
5.1)天然气***中管道分段
管道分段数越多，差分方程与原偏微分方程结果越相近，但相应的计算量越大，计算时间越长。本算例中选取管道分段数为10.
5.2)天然气***暂态能流方程
根据公式(16)至公式(22)可以列写一组非线性方程。本算例中共有125个方程与125个待求状态变量。
6)牛顿法求解天然气***能流
6.1)计算不平衡量
将步骤(1)中输入的天然气***网络参数以及设置的各状态变量初值代入式23，可以计算得到不平衡量：
Figure GDA0001769116170000211
6.2)计算雅可比矩阵
计算得到雅可比矩阵如下
Figure GDA0001769116170000212
注意，该雅可比矩阵为125阶方阵，且不是01矩阵，上式仅仅展示了雅可比矩阵中部分元素的值。
6.3)计算修正量
根据上述步骤获得的不平衡量△FG和雅克比矩阵JG，以第一次迭代后的结果为例，计算出各状态变量的修正量：
Figure GDA0001769116170000213
根据上述结果更新状态变量：
Figure GDA0001769116170000214
6.4)收敛性判断
当修正量ΔXG满足max(|ΔXG|)<ε，则计算结束。
当max(|ΔXG|)≥ε，则返回步骤6.1)继续迭代计算。
根据前面步骤，计算可知当迭代3次后，满足收敛条件，此时max(|ΔXG|)＝1.55E-9＜ε。
7)多时段电-气互联***能流计算
当仿真时间t<Tmax时，执行t＝t+Δt，将能流计算结果作为下一次能流计算中各状态变量的初值，并返回执行步骤2)～6)，计算下一时刻电气互联***能流。
当仿真时间t>Tmax时，能流计算完毕，结束运算。
实验效果：
以附图3所示的以IEEE14节点电力***与6节点天然气***联合而成的电-气互联***为仿真对象，设计以下仿真算例并通过三种方法对比验证本发明方法的优越性和有效性。算例：初始时刻时，电气互联***处于稳态，此时电力***出现负荷波动，波动后电力***负荷数据如表1所示，考虑多平衡机时各发电机承担份额如表2所示。计算该负荷波动后电气互联***能流。方法1：考虑电力***多平衡机但忽略天然气***慢动态特性对电气互联***进行能流计算。方法2：考虑电力***单平衡机与天然气***慢动态特性对电气互联***进行能流计算。方法3：考虑电力***多平衡机与天然气***慢动态特性进行能流计算。
表1负荷波动后电力***各节点负荷
节点编号 有功功率(MW) 无功功率(MVar)
1 0 0
2 21.70 12.70
3 99.20 19.00
4 52.80 -3.90
5 7.60 1.60
6 11.20 7.50
7 0 0
8 0 0
9 29.50 16.60
10 9.00 5.80
11 3.50 1.80
12 6.10 1.60
13 13.50 5.80
14 14.90 5.00
表2考虑多平衡机时各发电机的承担份额
发电机节点编号 份额
1 60％
2 0
3 30％
6 10％
8 0
图4中的两条曲线分别为方法1和方法3计算得到负荷波动后天然气***节点5处压力的变化曲线。方法1中忽略了天然气***的慢动态特性，在负荷波动后，节点5处的压力立即达到了稳态值3.13MPa。由图可见在前300min，方法1与方法3的能流计算结果相差较大，忽略天然气***的慢动态特性计算能流所得结果与***实际状态会有较大差异，因此，有必要考虑天然气***慢动态特性进行能流计算。
表3方法3计算负荷波动前后各发电机输出有功功率及***有功损耗
发电机节点编号 负荷波动前有功功率(MW) 负荷波动后有功功率(MW)
1 226.66 233.07
2 40.00 40.00
3 3.83 7.04
6 1.28 2.35
8 0 0
通过方法3计算可得，负荷波动前***的有功功率损耗为12.76MW，负荷波动后***的有功功率损耗为13.45MW，负荷波动造成***的功率缺额为10MW，因此，***波动后总的不平衡功率为10.69MW。表3为负荷波动前后通过方法3计算得到电力***各发电机节点所发出的有功功率。可见节点1处发电机承担了60％的不平衡功率，增发了6.41MW的有功功率，节点3与节点6分别增发了不平衡功率30％与10％数值的有功功率，可见各发电机严格按照份额来承担***不平衡功率，进而验证了本方法的正确性。
表4方法2与方法3计算所得发电机输出有功功率
Figure GDA0001769116170000231
表4为方法2与方法3计算所得负荷波动后电力***发电机输出的有功功率，方法2中节点1处的发电机承担了全网的功率损耗以及负荷波动所造成的***不平衡功率，而方法3中由节点1、节点3、节点6处的发电机共同承担，两种方法计算得到各发电机的输出功率差异较大，实际***中的不平衡功率是多台发电机共同承担的，通过方法3计算所得到的能流结果更符合实际情况。因此能流计算时有必要充分考虑电力***多平衡机机制。
由此可见，在负荷波动、断线故障等造成***较大不平衡功率的情况下，电气互联***能流计算中有必要同时考虑电力***多平衡机与天然气***慢动态特性，得到更准确的***运行状态。

Claims (6)

1.一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法，其特征在于，包括以下步骤：
1)输入所述电力***和天然气***的参数，并设定总仿真时间Tmax；建立考虑多平衡机的电力***潮流模型；
建立考虑多平衡机的电力***潮流模型的步骤如下：
1.1)设定电力***；所述电力***中有Ne个节点；其中，Vδ节点数目为1，PV节点数目为r，PQ节点数目为Ne-r-1；电力***中有r+1台发电机组成多平衡机组；多平衡机组共同承担***的不平衡功率；
2)初始化电力***和天然气***的参数；
电力***中Vδ节点的有功功率初值设定原则为：不考虑功率损耗的情况下，电力***发电机输出有功功率等于负荷有功功率；
PQ节点电压的初值以平启动方法设置，即电压幅值标幺值初值设置为1，相角初值设置为0；
PV节点的电压相角初值设置为0；
天然气***中各节点的压力为天然气***稳态情况下各节点的压力值；天然气***管道分段后新增节点的压力初值设置为管道始末节点压力的等差数列值；管道的流量设置为稳态情况下管道中的流量；
3)建立考虑多平衡机的电力***潮流模型；
考虑多平衡机的潮流方程式如下所示：
Figure FDA0003585555830000011
式中，
Figure FDA0003585555830000012
表示电力***中除Vδ节点外其余节点编号所在集合；ζ表示电力***中PQ节点编号所在集合；PGi为节点i处发电机所发有功功率；PDi为节点i处负荷有功功率；Pical为计算得到节点i处传输线上的有功功率；QGj为节点j处发电机所发无功功率；QDj为节点j处负荷无功功率；Qjcal为计算得到节点j处传输线上的无功功率；αi为节点i处发电机的功率分配系数；ΔPΣ为***不平衡功率；
若节点i上无发电机则PGi＝QGi＝0；
若节点i上无负荷则PDi＝QDi＝0；
ΔPi为电力***节点i处有功功率不平衡量；ΔQi为电力***节点i处无功功率不平衡量；
节点i处传输线上的有功功率Pical和节点i处传输线上的无功功率Qical分别如下所示：
Figure FDA0003585555830000021
式中，
Figure FDA0003585555830000022
为节点i的电压；
Figure FDA0003585555830000023
为节点j处电压的共轭值；
Figure FDA0003585555830000024
为节点i和节点j的互导纳矩阵Yij的共轭值；i、j为电力***节点；
功率分配系数αi如下所示：
Figure FDA0003585555830000025
若节点i处无发电机或发电机不参与承担***不平衡功率，则功率分配系数αi为0；
功率不平衡ΔPΣ如下所示：
Figure FDA0003585555830000026
式中，Ploss为***中的功率损耗；PGi为节点i处发电机所发有功功率；PDi为节点i处负荷有功功率；Ne为电力***节点总数；i为任意节点编号；
***中的功率损耗Ploss如下所示：
Figure FDA0003585555830000027
式中，Gij为节点i与节点j的互电导，δi为节点i处的相角；δj表示节点j处的相角；Ne为电力***节点总数；i,j为任意节点编号；Ui为节点i的电压幅值；
2)利用牛顿法求解电力***的潮流分布；
3)确定电-气耦合设备中天然气的产生量和消耗量，从而建立天然气***暂态能流模型；
4)利用牛顿法求解天然气***的能流分布；
5)判断是否达到总仿真时间Tmax，若未达到总仿真时间，则将天然气***的能流计算结果作为下一时刻能流计算初值，重复步骤2)至步骤4)，若达到总仿真时间，则计算结束并输出电力***的潮流分布和天然气***的能流分布结果。
2.根据权利要求1所述的一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法，其特征在于：所述电力***与天然气***参数包括电力***网络参数和天然气***参数；
所述电力***网络参数包括电力***支路参数、变压器参数、发电机参数、负荷参数和节点类型；
所述天然气***参数包括气源参数、输气管道参数、各节点负荷参数和平衡节点压力初值。
3.根据权利要求1所述的一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法，其特征在于：利用牛顿法求解电力***潮流的步骤如下：
1)计算电力***功率不平衡量矩阵ΔF；
电力***功率不平衡量矩阵ΔF如下所示：
Figure FDA0003585555830000032
式中，n＝Ne-1；Ne为电力***节点总数；△P为电力***有功功率不平衡量，△Q为电力***无功功率不平衡量；电力***功率不平衡量矩阵ΔF的前2(n-r)个元素为根据PQ节点计算的功率不平衡量，后r个元素为根据PV节点计算的功率不平衡量；
2)计算电力***的雅可比矩阵J；
雅可比矩阵J如下所示：
Figure FDA0003585555830000031
式中，电力***雅可比矩阵维数为(2n-r)×(2n-r)；
雅可比矩阵中的元素Hij、Jij、Nij和Lij分别如公式(8)所示：
Figure FDA0003585555830000041
式中，Hij和Jij分别为有功功率不平衡量和无功功率不平衡量对电压相角不平衡量求偏导所得值；Nij和Lij分别为有功功率不平衡量和无功功率不平衡量对电压幅值不平衡量求偏导所得值；ΔPi为电力***节点i处有功功率不平衡量；ΔQi为电力***节点i处无功功率不平衡量；δi为节点i处的相角；Uj为节点j的电压幅值；
3)计算各节点电压的修正量；
根据功率不平衡量矩阵△F和雅可比矩阵J，计算各节点电压的修正量，从而得到电压修正量矩阵ΔX；所述电压修正量矩阵ΔX中包括电压相角修正量和电压幅值修正量；电压修正量矩阵ΔX如下所示：
ΔX＝J-1ΔF； (9)
式中，J-1为雅可比矩阵J的逆矩阵；△F为功率不平衡量矩阵；
Figure FDA0003585555830000042
式中，上标T表示转置；n＝Ne-1；Ne为电力***节点总数；功率不平衡量矩阵ΔF的前2(n-r)个元素为根据PQ节点计算的功率不平衡量，后r个元素为根据PV节点计算的功率不平衡量；δi为节点i处的相角；ΔUi为节点i的电压幅值修正量；Δδi为节点i的电压相角修正量；
根据电压修正量更新节点电压向量X：
X(k+1)＝X(k)+ΔX(k)； (11)
式中，k为迭代次数；当k＝1时，X(k)中各元素为电力***各状态变量初值；
4)收敛性判断
当max(|ΔX|)≥ε，则返回步骤1)继续迭代计算；ε为牛顿法收敛精度；
当不平衡量ΔX满足max(|ΔX|)<ε，则计算结束。
4.根据权利要求1或2所述的一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法：根据电-气耦合设备计算等效天然气负荷量的步骤如下：
1)确定电-气互联***的电气耦合方式，将天然气转化为电能；所述电气耦合方式为利用燃气机组进行耦合；
2)设定电力***节点e处接有燃气机组，则燃气机组输出的有功功率Pe和燃气机组所消耗天然气的质量流量GfGPG,e的关系如下所示：
Figure FDA0003585555830000051
式中，ηGPG为燃气机组的效率；LHV为天然气低热值；ρ0为天然气密度。
5.根据权利要求1所述的一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法，其特征在于：建立天然气***暂态能流模型的步骤如下：
1)设定天然气***；天然气***中有NG个节点，压力已知的节点数目为Nt，管道数目为NL，管道分段数为Np，压缩机数目为Nc；
2)对天然气***中的天然气管道进行分段；
设定所述天然气***中包括m条管道；每条管道分为等长度的Np段；
记第L管道的首末节点的压力分别为p1、p2；记第L管道分段后从首端至末端Np+1个节点的压力分别为：
Figure FDA0003585555830000052
3)建立天然气***中燃气轮机驱动压缩机的模型
压缩机消耗的能量BHP如下所示：
Figure FDA0003585555830000053
式中，pout为压缩机出口压力；pin为压缩机入口压力；Gfcom为流经压缩机的流量；B为和天然气温度、压缩机效率、天然气热值有关的常量；z为与天然气热值和天然气压缩因子有关的常量；
当压缩机由天然气驱动时，压缩机所消耗的天然气流量τ如下所示：
Figure FDA0003585555830000061
式中，α、β和γ均为能量转换过程中的效率系数；BHP为压缩机消耗的能量；
压缩机压缩比Rc如下所示：
Figure FDA0003585555830000062
式中，pout为压缩机出口压力；pin为压缩机入口压力；
4)建立天然气***暂态能流方程
4.1)对于任意一条管道，均列写四类非线性方程，则第L条管道的各节点功率不平衡量矩阵ΔF分别如下所示：
Figure FDA0003585555830000063
式中，
Figure FDA0003585555830000064
为暂态常数；
Figure FDA0003585555830000065
表示管道L第l个节点处的压力；
Figure FDA0003585555830000066
为管道L第l个节点处的压力初值，即该节点上一时刻的压力；
Figure FDA0003585555830000067
为管道L第l+1个节点处的流量；
Figure FDA0003585555830000068
为管道L第l-1个节点处的流量；Δt为时间步长；Δl为空间步长；
Figure FDA0003585555830000069
式中，
Figure FDA00035855558300000610
为暂态常数；
Figure FDA00035855558300000611
为管道L第l+1个节点处的压力；
Figure FDA00035855558300000612
为管道L第l-1个节点处的压力；
Figure FDA00035855558300000613
为管道L第l个节点处的流量；
Figure FDA00035855558300000614
式中，
Figure FDA00035855558300000615
为暂态常数；
Figure FDA00035855558300000616
为管道L第l+1个节点处的压力；
Figure FDA00035855558300000617
为管道L第l个节点处的压力；Gfl L为管道L第l个节点处的流量；
Figure FDA00035855558300000618
式中，
Figure FDA00035855558300000619
表示管道L第l个节点处的压力；Gfl L为管道L第l个节点处的流量；
Figure FDA00035855558300000620
为暂态常数；
Figure FDA00035855558300000621
为管道L第l+1个节点处的压力；
暂态常数
Figure FDA00035855558300000622
和暂态常数
Figure FDA00035855558300000623
分别如下所示：
Figure FDA0003585555830000071
式中，Z为管道中天然气的压缩系数；R为天然气的气体常数；
Figure FDA0003585555830000072
为天然气管道L的横截面积；T为天然气温度；
Figure FDA0003585555830000073
为管道L的摩擦因子；DL为天然气管道L的内径；
4.2)除天然气***压力已知节点之外的其他节点m的流量不平衡量ΔGm如下所示：
Figure FDA0003585555830000074
式中，Gfv→m为节点v流向节点m的流量；Gfm-source为气源注入节点m的流量；Sc为符号函数；当m为压缩机出口节点时，Sc为1；当m为压缩机入口节点时，Sc为-1；Gfm-load为节点m处气负荷消耗的流量；NG为压力天然气***节点个数；Nt为天然气***压力已知节点个数；τ为压缩机所消耗的天然气流量；Gfcom为流经压缩机的天然气流量；
压缩机压力约束方程ΔRc如下：
Figure FDA0003585555830000075
式中，Rc为压缩机c的压缩比；pout,c为压缩机c出口压力；pin,c为压缩机c入口压力；Nc为压缩机个数。
6.根据权利要求1或5所述的一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法，其特征在于：利用牛顿法求解天然气***能流的步骤如下：
1)计算天然气***的不平衡量ΔFG；
天然气***的功率不平衡量ΔFG如下：
Figure FDA0003585555830000076
式中，a∈{1,2,…,NL}；ΔGm为除天然气***压力已知节点之外的其他节点m的流量不平衡量；
2)计算雅可比矩阵JG；
Figure FDA0003585555830000081
式中，b∈{1,2,…,NG+1}；
雅可比矩阵JG为维数为[2NpNL+(NG-Nt+Nc)]×[2NpNL+(NG-Nt+Nc)]的方阵；天然气***中有NG个节点，压力已知的节点数目为Nt，管道数目为NL，管道分段数为Np，压缩机数目为Nc；
3)计算修正量
根据不平衡量△FG和雅克比矩阵JG，计算待求变量的修正量ΔXG；所述待求变量为天然气***能流，包括天然气***的压力和流量分布；
修正量ΔXG如下所示：
Figure FDA0003585555830000082
式中，JG为雅可比矩阵；ΔFG为天然气***的不平衡量；
Figure FDA0003585555830000083
式中，Δpm为***节点压力修正量；
Figure FDA0003585555830000084
为***新增节点压力修正量；
Figure FDA0003585555830000085
为管道各处流量修正量；ΔGfcom,c为压缩机耗气量修正量；
根据公式(26)，对XG进行更新，即：
XG (k+1)＝XG (k)+ΔXG (k)； (27)
式中，k表示迭代次数；ΔXG为修正量；XG (k)为天然气***能流；
4)收敛性判断
当max(|ΔXG|)≥ε，则返回步骤1)继续迭代计算；ε为牛顿法收敛精度；
当修正量ΔXG满足max(|ΔXG|)<ε，则计算结束。
CN201810745731.7A 2018-07-09 2018-07-09 一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法 Expired - Fee Related CN108875269B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810745731.7A CN108875269B (zh) 2018-07-09 2018-07-09 一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810745731.7A CN108875269B (zh) 2018-07-09 2018-07-09 一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法

Publications (2)

Publication Number Publication Date
CN108875269A CN108875269A (zh) 2018-11-23
CN108875269B true CN108875269B (zh) 2022-06-10

Family

ID=64300265

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810745731.7A Expired - Fee Related CN108875269B (zh) 2018-07-09 2018-07-09 一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法

Country Status (1)

Country Link
CN (1) CN108875269B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3989118A4 (en) 2019-06-28 2022-06-29 SZ DJI Technology Co., Ltd. Target tracking method and system, readable storage medium and moving platform
CN111429302B (zh) * 2020-04-14 2023-07-25 南京工程学院 一种综合能源***稳态能流计算中天然气***初值计算方法
CN111400937B (zh) * 2020-04-29 2023-12-01 上海电力大学 一种综合能源***潮流计算方法
CN112329217B (zh) * 2020-10-23 2022-07-15 重庆大学 一种电-水综合能源***的准稳态能流分析方法
CN112818492B (zh) * 2021-02-02 2022-12-02 山东大学 计及能源集站状态变量的电-气耦合网络能量流求解方法
CN113410867B (zh) * 2021-06-30 2022-06-14 华南理工大学 一种基于全纯函数的气电联合网络动态仿真计算方法
CN114492958B (zh) * 2022-01-07 2024-06-04 中国能源建设集团广东省电力设计研究院有限公司 一种基于多能流计算的能源管理方法和装置

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101139981A (zh) * 2007-10-19 2008-03-12 无锡压缩机股份有限公司 压缩机往复惯性力平衡结构
CN101685968A (zh) * 2009-07-24 2010-03-31 重庆大学 配电网可靠性评估的故障扩散方法
CN101702521A (zh) * 2009-11-20 2010-05-05 河海大学 计及多平衡机影响的电力***状态估计方法
CN105631710A (zh) * 2015-12-31 2016-06-01 谢卓甫 一种利用管道气体压力形成气体市场价格的方法
CN106874554A (zh) * 2017-01-13 2017-06-20 重庆大学 一种天然气***能流计算的改进方法
CN107291990A (zh) * 2017-05-24 2017-10-24 河海大学 基于电‑气互联综合能源***暂态模型的能量流仿真方法
CN107368927A (zh) * 2017-08-01 2017-11-21 重庆大学 基于目标级联分析的电‑气能量流分布式协同优化计算方法
CN107947245A (zh) * 2017-11-20 2018-04-20 重庆大学 考虑天然气***约束的等值最优潮流模型构建方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9099866B2 (en) * 2009-09-01 2015-08-04 Aden Seaman Apparatus, methods and systems for parallel power flow calculation and power system simulation

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101139981A (zh) * 2007-10-19 2008-03-12 无锡压缩机股份有限公司 压缩机往复惯性力平衡结构
CN101685968A (zh) * 2009-07-24 2010-03-31 重庆大学 配电网可靠性评估的故障扩散方法
CN101702521A (zh) * 2009-11-20 2010-05-05 河海大学 计及多平衡机影响的电力***状态估计方法
CN105631710A (zh) * 2015-12-31 2016-06-01 谢卓甫 一种利用管道气体压力形成气体市场价格的方法
CN106874554A (zh) * 2017-01-13 2017-06-20 重庆大学 一种天然气***能流计算的改进方法
CN107291990A (zh) * 2017-05-24 2017-10-24 河海大学 基于电‑气互联综合能源***暂态模型的能量流仿真方法
CN107368927A (zh) * 2017-08-01 2017-11-21 重庆大学 基于目标级联分析的电‑气能量流分布式协同优化计算方法
CN107947245A (zh) * 2017-11-20 2018-04-20 重庆大学 考虑天然气***约束的等值最优潮流模型构建方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Influence of Power-to-Gas-Technology on Unit Commitment and Power System Operation;I. Marjanovic;《IEEE》;20071231;第536-540页 *
Multi-time period combined gas and electricity network optimisation;Modassar Chaudry 等;《Electric Power Systems Research》;20081231;第1265-1279页 *
Steady-state analysis of the integrated natural gas and electric power system with bi-directional energy conversion;Qing Zeng 等;《Applied Energy》;20160518;第1483-1492页 *
多平衡机潮流计算在调度员潮流中的应用;崔丽华 等;《广东电力》;20060930;第13-17页 *
电-气互联综合能源***多时段暂态能量流仿真;卫志农 等;《电力自动化设备》;20170630;第41-47页 *
获取电力***运行方式的多平衡机潮流模型;杜正春 等;《电力***自动化》;20140825;第41-47页 *

Also Published As

Publication number Publication date
CN108875269A (zh) 2018-11-23

Similar Documents

Publication Publication Date Title
CN108875269B (zh) 一种考虑电力***多平衡机与天然气***慢动态特性的电-气互联***多时段能流计算方法
CN107291990B (zh) 基于电-气互联综合能源***暂态模型的能量流仿真方法
CN111061986B (zh) 一种多运行模式的热电综合能源***潮流计算方法
CN107947245B (zh) 考虑天然气***约束的等值最优潮流模型构建方法
CN109636052B (zh) 一种气电联合***的协同规划方法
CN110601203A (zh) 一种电-气耦合***分段线性化最优潮流计算方法
CN107069741A (zh) 一种新型线性化潮流计算方法
CN112531716A (zh) 电-水互联***混合能流的统一标幺计算方法
CN110866213A (zh) 电-气综合能源***多网络稳态能量流分析方法及装置
CN111523247A (zh) 基于分立求解的综合能源***多能流联合计算方法及装置
CN111400873A (zh) 电热综合能源***二阶锥规划抗差状态估计方法和***
Larsen et al. Equivalent models for district heating systems
Yang et al. A dynamic linearization and convex relaxation-based approach for a natural gas optimal operation problem
CN115062555A (zh) 一种基于非平衡节点㶲的综合能源***㶲流直接计算方法
CN105302948B (zh) 单相流体网络建模方法
CN114221346A (zh) 一种综合能源***潮流计算方法
CN113627021A (zh) 一种基于序列凸规划的电气互联***最优能流计算方法
Bazrafshan et al. A fixed-point iteration for steady-state analysis of water distribution networks
CN103701125B (zh) 一种基于序列二次规划法的配电网柔性潮流算法
CN106410811B (zh) 首次迭代小阻抗支路端点改变雅可比矩阵的潮流计算方法
CN111881492B (zh) 一种基于cfd方法的汽轮机阀门配汽管理函数生成方法
CN112257355A (zh) 一种掺氢天然气的中低压配气管网建模方法
CN108964143B (zh) 一种电-气综合能源***的天然气网络静态等值模型
CN104598728B (zh) 一种计及频率变化的含风力发电的电力***状态估计方法
Xianxi et al. Modeling and simulation of steam pipeline network with multiple supply sources in iron& steel plants

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20220610