张 晴

(北京邮电大学 理学院，北京 100876)

随着天然气管道系统的规模越来越大，越来越复杂，导致这些设施的优化运行和规划的重要性也越来越大.此外，管网的投资成本和运行费用非常大，即使是系统利用率的小幅提高，也会涉及大量资金.

天然气行业的服务包括生产、输送和销售天然气.本研究的主要兴趣集中在通过管网输送天然气方面.在输送天然气的过程中，气体流经各种设备，包括管道及压缩机等管网元素，气体在管道中流动时，由于与管壁的摩擦，导致气体压力降低，不能稳定传输，因此需在管网的适当位置建立压气站给气体增压，使之稳定传输.天然气管网系统的总运行成本很大程度上取决于管网中压气站的运行成本.而压气站的运行成本，一般是以压气站消耗的燃料来衡量的.根据Luongo等人[1]的研究，压气站的运行成本占公司总运行预算的25%～50%.因此，传输网络的目标是在满足规定的输送流量和输送终端最低压力和温度要求的前提下，尽量减少压气站的总燃料消耗[2].

根据气体流量随时间的变化情况，管网系统可分为稳态和瞬态.当系统中表征气体流动的数值与时间无关时，系统被称为稳态.在此情况下，系统的约束条件，特别是描述气体流经管道的约束条件，可用代数非线性方程来描述.本文主要讨论稳态输气网络优化问题.

最近20年的工作旨在解决所有网络元素的非线性和离散方面.然而，由于目前的优化算法仍然无法解决真实网络的完整混合整数非线性规划(MINLP)模型，因此主要采用线性近似的方法，以获得混合整数线性规划(MILP)模型(参见Hacklnder[3]，2002年；Martin和Möller[4]，2005年；Martin等[5]，2006年；GeiBler 和Martin2015[6]，2015年；黄亚魁等[7]，2017年；Robert Burlacu等[8]，2019年).

1 管网模型

本节将介绍天然气运输网络优化中的一个基本问题.考虑一个天然气运输网络和预先已知的需要输送的天然气数量，即在每个入口，天然气可以进入网络；而在每个出口，天然气可以从网络中输出.此外，由于长距离输送管道上的摩擦引起的压降，必须考虑可增加气体压力的压缩机.上述内容导致了在找到满足所有物理和技术约束的可控元件的可接受配置的同时，最小化压缩机燃料消耗量的目标.由于气体物理学是通过非线性函数来描述的，考虑到可控元素的离散特性，该问题可被建模为MINLP.

1.1 稳态管网的MINLP模型

为了模型的简化，本文讨论的管网仅包含节点、管道和压缩机.

使用一个有向有限图G=(V,A)来模拟气体网络.节点集V由入口节点集Vs、出口节点集Vd和内部节点集Vi组成.弧的集合A由管道的集合Ap和压缩机的集合Ac组成.此外，弧a=(u,v)∈A上的质量流量用qa表示，节点v∈V处的压力和温度分别用pv和Tv来表示.考虑到一个规范的网络，对于每个节点v∈V，需给定压力变量pv和温度变量Tv的上下界约束：

(1)

(2)

此外，每个节点v∈V处需满足质量流量守恒：

(3)

其中：dv为v处的天然气需求量，满足dv≥0，若v∈Vd；dv≤0，若v∈Vs；dv=0，若v∈Vi.

1.1.1 管道

我国燃气管网中最常见的部件是管道.每一根管道a=(u,v)∈AP由其长度La、直径Da和粗糙度ka决定.管道内的气体流动可以用质量流量q、压力p、温度T、密度ρ和欧拉方程体系与具体状态方程来描述.在本文中，状态方程由真实气体的热力学标准方程给出：

ρRszT=p，

其中:Rs为通用气体常数，z为压缩系数.后者可用下式所示：

上式中的假临界温度pc和假临界温度Tc均为常数.通过一定的简化，管道上的压降可以用著名的Weymouth方程[9]以解析形式给出：

(4)

(5)

此外，管道上的流量上下界由以下约束给出：

(6)

在考虑焦耳-汤姆逊效应的影响下，管道上的温降可以用著名的Sukhov公式加上焦耳-汤姆逊效应的影响项[10]以解析形式给出：

∀a=(u,v)∈Ap，

(7)

1.1.2 压缩机

为了补偿气网中的压降，管网中安装了将入口压力提高到较高出口压力的压缩机.每台压缩机a=(u,v)∈Ac由其绝热效率ηad,a和能耗成本系数ca指定.本文仅考虑涡轮压缩机，其可行工作范围由其特性图来描述，见图1中的灰色区域[11].

图1 涡轮压缩机的特性图Figure 1 Characteristic diagram of turbo compressor

其中：Qa代表体积流量，与qa有如下关系：

(8)

其中：ρ0为标况下的气体密度.此外，压缩机的绝热过程会导致绝热焓Had,a的变化：

(9)

其中：κ是绝热指数，假定为常数.此时压缩系数za可近似为[10]：

(10)

压缩机将进口压力pu提高到较高的出口压力pv所消耗的功率Na为：

(11)

当压缩机关闭时，没有气体流经压缩机，其入口和出口的压力和温度是解耦的，显然也不消耗功率.当压缩机运行时，其在消耗功率的同时也增加了沿弧线方向流动的气体的压力和温度，即qa≥0，pv≥pu，Tv≥Tu，逆压缩机弧线方向的流动是不可能的.引入二元变量sa来表示压缩机的不同状态：sa=1表示压缩机a是开的，反之亦然.除了能耗约束(12)外，模型中还加入了以下约束条件：

(12)

(13)

(14)

(15)

puybmin≤pv≤puybmax，∀a=(u,v)∈Ac，

(16)

sa∈{0,1}，a∈Ac.

(17)

1.2 MINLP问题的总结

管网的运行成本主要为压缩机的燃料消耗，因此，本文假设每台压缩机a∈Ac的能耗成本系数ca≥0.考虑到前两小节中描述的燃气网络中的不同元素，可得到下述MINLP问题：

s.t.(1)～(17).

2 网络约简、线性化及算法框架

本节将介绍一个用于构建MINLP模型的MILP近似的方法及算法框架，以解决第二节中的问题.本节讨论的内容是文献[7]中的方法，为了内容的完整性，在此节进行简要回顾.

2.1 网络约简

用虚拟压缩机代替顺序连接的压缩机，该虚拟压缩机的行为表现为其组件以最佳方式运行[12].此处使用的约简方法为只替换顺序相连的压缩机站，如图2转变为图3所示.

图2 循环网络Figure 2 Cyclic network

图3 约简网络Figure 3 Reduced network

2.2 线性化-凸组合法

对原始网络进行约简后的简化网络的规模虽然会减少，但由于目标函数和某些约束的非凸非线性，导致直接调用MINLP求解器得到的解的质量不会很高.其中一种处理该问题的方法是基于建立混合整数线性模型的思想.用到的线性化方法是凸组合法[13]，其具体形式为：

λj≥0,j=1,…,m;zi∈{0,1},i=1,…,n.

其中S:={S1,…,Sn}为非线性函数y(x)的可行域剖分成的单纯形的集合，m是S中的顶点个数.此外，引入了辅助变量λ和z及其相应的线性约束，其中z是0-1变量.

2.3 算法框架

求解稳态燃气网络的基于线性化的算法如下所述：

算法3.1 优化燃气网络的基于线性化的算法

步骤1 网络约简

1)将网络分解为仅包含顺序压缩器的小型子网站；

2)用虚拟压缩器替换每个子网并获得简化的网络；

步骤2 建立管网的MINLP模型

步骤3 线性化

1)交替固定流量和其他变量，并选择离散点；

2)优化每个子网以获得每个离散点的能耗和出口温度；

3)使用可行网格点划分可行区域；

4)建立MILP模型；

步骤4 解决MILP模型

步骤5 如果通过步骤4获得的解满足一定条件，则停止算法；若不满足条件，转到步骤3.

3 等价替换—连续化方法

模型中的整数变量的存在加大了问题的求解难度，而已知的0-1规划的求解方法可以分成离散方法和连续化方法两类.

离散方法是直接解整数规划的解法，连续化方法就是把问题中的整数变量变为连续变量之后进行求解，最简单的就是把0-1规划中的0-1整数变量松弛到[0,1]区间上[14].本节的内容就是通过在模型中添加一个关于0-1变量的非线性等式约束，从而把原问题中的一部分0-1变量变为连续变量，并证明了进行替换后的问题与原问题的等价性，与此同时利用罚函数将约束吸收到目标函数.该替换可以降低问题求解难度，并提高问题求解速度.

把1.1.2节中出现的变量sa通过如下所述的连续化方法变为[0,1]区间上的连续变量，对∀j均满足下式：

证明：首先证明左式推右式：

由sa,j=0/1⟹sa,j(1-sa,j)=0⟹∑jαjsa,j(1-sa,j)=0，且满足0≤sa,j≤1.

现在证明右式推左式：

∵0≤sa,j≤1，故0≤sa,j(1-sa,j)≤0.25，

又∵∑jαjsa,j(1-sa,j)=0⟹对∀j，都有αj>0，故对∀j,sa,j(1-sa,j)=0均成立，⟹sa,j=0/1，证毕.

因此，可以把模型中的整数变量连续化为[0,1]区间上的连续变量，进而减小问题的求解难度，并加快问题求解速度.

4 数值结果

将前两节中的算法应用于我国西部天然气管网系统，报告数值结果.

如图4所示，该测试管网由27个压缩机站，32条管道，64个节点和2个循环组成[7].输气量为1.4×106m3/d，压缩机出入口的压力范围均是[8,12](MPa)，温度取值范围为[0,50](℃).在图4中，竖着的黑色梯形块符号代表的是压缩机，箭头表示的是气体的流动方向.图5是相应的约简管网.

图4 测试管网Figure 4 Test network

图5 约简管网Figure 5 Reduced network

建立约简管网的MINLP模型，采用算法3.1求解该模型，并将此结果与使用等价连续化的结果进行比较.为了表达的便利，算法3.1称为方法1，对整数变量进行等价连续化的算法3.1称为方法2.

本节中介绍的所有计算都是在Intel Core i7，1.80 GHz的处理器上完成的.在每一次迭代中，均需选择离散点，并对每一子网进行优化.测试中使用DP优化子网，使用建模语言YALMIP[15]为约简管网构建MILP模型.作为MILP求解器，本文使用CPLEX12.8.0[16]，其中所有参数都设置为默认值.

可看到可行的解决方案是通过三次迭代获得的.如前所述，简化网络仅包含虚拟压缩机.因此，需测量耗气量与出口温度的准确性.当误差小于给定的公差ε时，则接受该解.在实验中，设置误差精度ε=5%，超过两次迭代后，精度取其一半.表1中给出两种方法得到的解.

表1 MILP约简网络解的质量比较Table 1 Quality comparison of MILP approximate network solutions

从表1中可看出，第一次迭代时，除了VCS5子网的能耗误差稍大于给定误差外，其余子网的误差均达到给定精度，再根据表中的运行时间，表明两种方法均可在较短时间内得到可行解，这对于实际管网优化问题来说十分重要，而且方法2的用时比方法1的用时短.第三次迭代后，耗气量和出口温度均满足给定的误差精度，说明两种方法都可以提供令人满意的解决方案.此外，方法2获得的耗气量较低，为22 711.94 m3/d，比方法1的耗气量降低了133.85 m3/d.两者均提供了一种子网VCS2和VCS3关闭的情况下的运行方案.

跟踪约简管网的解可得到原始网络的解，表2显示出相应结果，其中粗体说明对应的压缩机是关闭的，可看到在满足管网各物理和技术约束的条件下，两种方法均关闭了15个压缩机，但每个压缩机的出入口压力、温度不尽相同.换句话说，通过对管道中的气体进行适当的增加或减少气体压力或者温度可达到降低管网系统耗气量的效果.由此可见，管网中各变量值的微小改变可带来非常可观的经济效益.

表2 MILP原始网络解的比较Table 2 Comparison of MILP original network solutions

综上所述，上节提出的方法可有效减少管网的耗气量，运行时间也有所减少.

5 结 语

本文提出了针对模型中整数变量的等价替换形式，可进一步减少问题求解时间，甚至可影响CPLEX的分支过程得出更优的解.数值结果也证明了本文方法的有效性.下一步希望能找到更多的燃气管网实例来研究该方法有效性，同时进一步提高算法的求解效率.然而，实际生产中的管网系统的变量随时间是不断变化的(即关于时间的函数)，所以对非稳态的天然气管网运行优化的研究值得进一步探索.