基于混合整数线性优化算法的天然气管网优化

Natural Gas Pipeline Network Optimization Based on a Mixed Integer Linear Optimization Algorithm

摘要: 随着天然气需求的增加和管网结构的复杂化，如何高效地优化天然气管网的设计和运行成为了一个重要的研究课题。本文提出了一种基于混合整数线性优化算法的天然气管网优化方法，该方法通过建立天然气管网的数学模型，以用户满意度最大和管输能耗最小为优化目标，以提高管网的输送效率并降低能耗。结果表明，该模型能够有效提升管网的运行效率，减少能耗，并且具有较强的适应性。

Abstract: With the increasing demand for natural gas and the increasing complexity of pipeline networks, efficiently optimizing the design and operation of natural gas pipeline networks has become an important research topic. This paper proposes a natural gas pipeline network optimization method based on a mixed integer linear optimization algorithm. This method establishes a mathematical model of the natural gas pipeline network and aims to maximize user satisfaction and minimize pipeline energy consumption, thereby improving pipeline transmission efficiency and reducing energy consumption. Results show that the model can effectively improve pipeline network operational efficiency, reduce energy consumption, and exhibit strong adaptability.

1. 引言

天然气作为一种重要的能源，广泛应用于工业、家庭以及发电等领域。随着全球能源需求的不断增长和城市化进程的加快，天然气管网的规模和复杂性不断提升。天然气管网的优化设计，不仅涉及管道布局、流量分配、节点压力控制等多个因素，还需要考虑能源消耗、系统安全性和用户满意度等方面。因此，如何高效地优化天然气管网，成为了当前面临的重要课题。

在天然气管网的优化过程中，传统的优化方法，如线性规划、非线性规划等，虽然在一定程度上可以提供优化方案，但在面对管网复杂约束条件时，往往存在较大的局限性。这些方法难以同时处理多种非线性、离散化的约束条件，且求解过程往往需要较长时间，无法满足现代天然气管网对实时性和高效性的需求。

近年来，随着计算机技术的不断发展，许多学者[1]-[4]提出了不同的天然气管网优化方法。冷建成[5]对原油多级星式集输管网建立了非线性混合整数规划模型，并采用神经网络方法求解。Cafaro和Grossmann [6]提出了大规模非凸混合整数非线性规划，开发了一种分级解耦方法来确定长期规划策略，包括钻井计划、管道参数以及工艺处理设备的规模，以使项目的净现值最大化。徐宁等[7]综合考虑输送模式、体积、运输设备、供需能力及站场建设方案等约束，以天然气供应链平均日运行费用最小为目标函数，建立了混合整数线性规划模型。Tan和Barton [8]提出了一种多周期混合整数线性规划优化模型，该模型可以为利用橇装设备将伴生气或滞留气货币化的决策者确定最佳的动态分配和运营决策。然而，国内外主要致力于降低能耗、运行费用与提高管网输送能力，以系统运行费用最低为目标函数，对压力、流量、温度及设备参数等主要变量进行决策。尚无直接适用于天然气管网适应性分析的管网运行优化模型、算法和软件。

为了更好地应对这一挑战，本文提出了一种基于混合整数线性优化算法的天然气管网优化方法，并以最大化用户满意度和最小化管输能耗为优化目标，力求在多约束条件下实现管网设计和运行的最优配置。

2. 基于混合整数线性优化算法的天然气管网优化模型

2.1. 用户满意度

在天然气行业中，用户满意度定义为用户对于所接受管网系统供气服务的满意程度，用于表征管网系统的供气服务能力，可由下式表示：

$用户满意度\left(S\right)=供气数据/用户期望$ (1)

管网的供气服务通常包含以下几个方面：气量、压力、气质、供气策略、应急情况供气、价格等。为了衡量管网的用户满意度，建立了以下评价体系及指标[9]-[11]，如下图1。

Figure 1. User satisfaction evaluation system and indicators

图1. 用户满意度评价体系及指标

本文采取用气量满意度进行计算，第i个用户满意度通过式(2)进行计算。

$RE{L}_i=\frac{Q_i^{dem}}{qdmi{n}_i}$ (2)

式(2)中：$RE{L}_i$ 是第i个用户的满意度；$Q_i^{dem}$ 是天然气的实际获得量；$qdmi{n}_i$ 是天然气的用气需求量。

考虑到天然气管网通过含有多个用户，且各用户重要度存在差异，则系统用户满意度通过式(3)进行计算。

$REL={\displaystyle \sum _{i\in N_a}{\lambda }_i\cdot RE{L}_i}$ (3)

式(3)中：$REL$ 是用户的满意度；${\lambda }_i$ 表示第i个用户在整体中的权重；$RE{L}_i$ 是第i个用户的满意度。

2.2. 管输能耗

管输能耗计算公式通过式(4)进行计算。

$Z_{\min }={\displaystyle \sum _{i=1}^n{Q}_i{L}_i{F}_i}$ (4)

式(4)中：$Z_{\min }$ 为目标函数管道的能耗量总和，kgce；n为管道数量；$Q_i$ 为第i条管道的输气量，m³；$L_i$ 为第i条管道的长度，km；$F_i$ 为第i条管道的单位周转量能耗(简称“生产单耗”)，kgce/(m³∙km)。

目标函数中关键在于确定单位周转量能耗计算公式。

在优化模型求解过程中获得各条管道的输量$Q_i$ 与生产单耗$F_i$ 的拟合关系，是进行能耗计算与优化的关键环节。通过对各管道输量与历史能耗数据的收集、整理，拟合得到单耗与周转量的线性、分段线性以及非线性函数关系。已知周转量，通过计算可以获得输量，由此找到了周转量与生产单耗的关系即找到了输量与生产单耗之间的关系。以管道生产管理系统的数据库为基础，研究得出了一种能耗历史数据的处理方法，通过拟合获得天然气主干管道的生产单耗与周转量的关系。

在输气管道的能耗评价中，通常采用单耗作为评价指标，该指标可由耗气量、耗电量以及周转量计算得到，其计算公式通过式(5)进行计算。

$F=\frac{W_g\alpha +W_e\beta }{G}$ (5)

式(5)中：$F$ 为单耗，kgce/(m³∙km)；$W_g$ 为耗气量，m³；$W_e$ 耗电量，kW∙h；$\alpha$ 为天然气折标煤系数，取1.33 kgce/m³；$\beta$ 为电折标煤系数，取0.1229 kgce/(kW∙h)；$G$ 为管道周转量，m³∙km。

生产单耗–周转量关系拟合步骤如下：

步骤1：压缩机启机增压情况确认。根据能耗数据反映的启机情况及管道的基础信息判断管道的启机类型，分为常年启机增压管道、间歇启机增压管道、非增压管道。

步骤2：月平均单耗值计算。按照月份划分管道各月的平均单耗值，分析对比不同月份单耗差异。

步骤3：按照周转量划分台阶，计算管道的平均单耗值。

步骤4：对单耗–周转量关系进行拟合。通过对管道每个月不同周转量台阶单耗数据的分析，确定单耗–周转量关系拟合的月份划分及周转量范围划分的方法，开展一次和二次函数关系拟合。

下图2展示了某管道系统周转量与生产单耗的拟合结果。

Figure 2. Fitting results of turnover and unit consumption

图2. 周转量与单耗拟合结果示意图

2.3. 约束条件

天然气管网流量分配优化模型中需要满足气源进气量约束、管网流量约束、管网压力约束以及管网水力约束。下面分开介绍各项约束条件：

(1) 管网流量约束：

a. 进入节点的流量之和等于离开的流量之和。

${\displaystyle \sum _{\left(i,j\right)\in E}{x}_{ij}\left(t\right)}-{\displaystyle \sum _{\left(j,l\right)\in E}{x}_{jl}}\left(t\right)=0\forall i,j\in V$ (6)

式中，$x_{ij}\left(t\right)$ 表示t时刻，第i节点到第j节点的流量，10⁴ Nm³/day；$\left(i,j\right)$ 和$\left(j,l\right)$ 分别表示节点i到节点j之间的管道和节点j到节点l之间的管道。

b. 双向管道的流向约束，表示管道中最多有一个管道流向。

$y_{ij}\left(t\right)+y_{ji}\left(t\right)\le 1\forall \left(i,j\right)\in E_b$ (7)

$y_{ij}\left(t\right)\in \left\{0,1\right\}\forall \left(i,j\right)\in E_b$ (8)

式中，$E_b$ 表示双向管道的集合，$y_{ij}\left(t\right)$ 为控制t时刻双向管道流向的二元决策变量。

c. 气源点供气量之和等于供给需求点的气量之和。

${\displaystyle \sum _{j\in S}{x}_{sj}\left(t\right)}={\displaystyle \sum _{i\in D}{x}_{id}\left(t\right)}$ (9)

式中，$x_{sj}\left(t\right)$ 表示t时刻，第j个气源点供给的气量，10⁴ Nm³/day。

d. 气源点供给流量小于气源能力上限，供给需求点的流量小于需求点的需求量，网络中任何管道的流量不会超过管道的管输能力，且流量不能为负，上述约束如下式所示。

$x_{ij}\left(t\right)\le y_{ij}\left(t\right)C_{ij}\left(t\right)\forall \left(i,j\right)\in E_b$ (10)

$x_{ij}\left(t\right)\le C_{ij}\left(t\right)\forall \left(i,j\right)\in E-E_b$ (11)

式中，$C_{ij}\left(t\right)$ 为第t天管道$\left(i,j\right)$ 的管输能力，10⁴ Nm³/day，其为能力矩阵C的一个元素。此处将气源能力、需求点的需求量以及管道的管输能力统一采用能力矩阵C表示。

(2) 管网压力约束：

a. 管网各节点压力约束。

$P_{i,\min }\le P_i\le P_{i,\max }$ (12)

式中，$P_i$ 为第i个节点的压力，MPa，$P_{i,\min }$ 和$P_{i,\max }$ 分别为第i个节点压力的上下限，MPa。

需要指出的是，压缩机站和调节阀上下游节点除需满足式(12)所示节点压力约束外，还需要满足其上下游压力约束。

b. 压缩机站上下游压力约束。

$P_{c_i,}{}_{lower}\le P_{c_i,1}\le P_{c_i,2}\le P_{c_i,}{}_{upper}$ (13)

式中，$P_{c_i,1}$ 、$P_{c_i,2}$ 、$P_{c_i,}{}_{lower}$ 、$P_{c_i,}{}_{upper}$ 分别为第i个调节阀上游压力、下游压力。上游压力限制和下游压力限制，MPa。

(3) 管网水力约束：

$y_{ij}\left(p_i^2\left(t\right)-p_j^2\left(t\right)\right)=y_{ij}\left(\frac{\lambda Z{\Delta }_{*}TL}{C_0^2{D}^5}{x}_{ij}^2\left(t\right)\right)\forall \left(i,j\right)\in E_b$ (14)

$p_i^2\left(t\right)-p_j^2\left(t\right)=\frac{\lambda Z{\Delta }_{*}TL}{C_0^2{D}^5}{x}_{ij}^2\left(t\right)\forall \left(i,j\right)\in E-E_b$ (15)

式中，$p_i\left(t\right)$ 和$p_j\left(t\right)$ 分别为第t时刻第i和j节点的压力，MPa，λ为水力摩阻系数、Z为天然气在管输条件下的压缩因子，${\Delta }_{*}$ 天然气的相对密度，T为输气温度，K，L输气管道计算段长度，km，D为输气管内径，m，$C_0$ 为常数，其值随各参数单位而定，上述参数均为节点i到节点j之间的管道的参数。

由于管网水力约束为非线性约束，且方程组庞大、未知项众多，直接采用优化模型求解比较困难。为了简化约束便于模型求解，尝试将水力约束分段线性化处理[12] (图3)，具体流程如下所示：

步骤一：针对某一个具体的普通管道，根据其设定的天然气输量范围设定均匀划分为10段，$x_{ij`1}^k,k=1,2,\cdots ,11$ ；并根据起点压力的范围，分别取7个起点压力$p_i^m,m=1,2,\cdots ,7$ 。

步骤二：在每个流量范围内，再均匀划分为5段，共取6个流量点$x_{ij}^{k_l},l=1,2,\cdots ,6$ 。

步骤三：根据管输工艺计算公式求解出10 × 6 × 7 = 420组末端压力值$p_j^w$ 。w对应一组k_l，m。

步骤四：在每一个流量范围内，由6 × 7组起点压力、终点压力和管道内天然气流量($p_i^m,p_j^w,x_{ij}^{k_l}$ )，根据线性回归公式$f\left(x\right)=Ax+b$ ，拟合出一个管道系数，其中$f\left(x\right)=p_i^2-p_j^2$ ，$x=x_{ij}$ 。最终得到每段流量范围的线性函数，如下图所示。由于式中还包含平方项，在建立模型时另定义$P_i=p_i^2$ ，表示压力的平方。

Figure 3. Schematic diagram of piecewise linearization

图3. 分段线性化示意图

对于优化模型中，还需要将分段线性约束转换为连续的线性约束，此时采用0-1变量将分段线性函数转化为线性约束。对于一个n段线性函数$f\left(x\right)$ 的分点为$b^1\le b^2\le \cdots \le b^n\le b^{n+1}$ ，引入连续变量$w^k$ 和0-1变量$z^k$ ，将x和$f\left(x\right)$ 表示为：

$x={\displaystyle \sum _{k=1}^{n+1}{w}^k{b}^k}$ (16)

$f\left(x\right)={\displaystyle \sum _{k=1}^{n+1}{w}^{k}f\left(b^k\right)}$ (17)

连续变量$w^k$ 和0-1变量$z^k$ 满足下列约束：

$\{\begin{array}{l}{w}^1\le z^1,w^2\le z^1+z^2,\cdots ,w^n\le z^{n-1}+z^n,w^{n+1}\le z^n\\ {\displaystyle \sum _{k=1}^{n+1}{w}^k=1}\\ {\displaystyle \sum _{k=1}^n{z}^k=1}\end{array}$ (18)

因此，通过上述方法，可以将非线性的管网水力约束$p_i^2\left(t\right)-p_j^2\left(t\right)=\frac{\lambda Z{\Delta }_{*}TL}{C_0^2{D}^5}{x}_{ij}^2\left(t\right)$ 均可转化为线性约束，如下所示：

$x_{ij}\left(t\right)={\displaystyle \sum _{k=1}^{11}{w}_{ij}^k{b}_{ij}^k}$ (19)

$P_i-P_j={\displaystyle \sum _{k=1}^{11}{w}_{ij}^{k}f\left(b_{ij}^k\right)}$ (20)

$\{\begin{array}{l}{w}^1\le z^1,w^2\le z^1+z^2,\cdots ,w^{10}\le z^9+z^{10},w^{11}\le z^{10}\\ {\displaystyle \sum _{k=1}^{11}{w}^k=1}\\ {\displaystyle \sum _{k=1}^{10}{z}^k=1}\end{array}$ (21)

式中，各参数含义如上所述，下标ij代表节点i至节点j的管道，$b_{ij}^k=x_{ij}^k,f\left(b_{ij}^k\right)=f\left(x_{ij}^k\right),k=1,\cdots ,11$ 。

决策变量：管道流量x_ij、双向管道的流向y_ij以及管道两端节点压力P_io。

3. 算例应用

3.1. 管网概况

为验证开发方法的可行性，以神府区块管网作为算例进行分析。本文将测试的天然气管道系统视为有向加权图，$G=\left(V,E\right)$ ，其中V代表图中的节点，E代表代表图中的弧，其拓扑结构如图4所示。天然气管网中气源视为源节点(Source)，各需求点视为汇节点(Sink)，压缩机站视为一般节点，天然气管道则视为包含管输能力约束的网络弧。

Figure 4. Pipeline network topology of the shenfu block

图4. 神府区块管网拓扑结构图

3.2. 基础参数

根据现场调研可获得优化模型的详细数据如下：

已建神府南区增压站压缩机出口能够提供的压力为4.5 MPa~5.5 MPa，外输神木首站压力为5.0 MPa。

神府中区增压站压缩机进口压力0.3 MPa~1.2 MPa，压缩机出口压力5.5 MPa~6.7 MPa。

拟建的府谷1#集气站和黄甫1#集气站，无相关压缩机参数，长输管道沿线站场压缩机为了便于采购、调配，通常设计参数统一一致，因此2座站场压缩机参数与神府中区增压站一致：进口压力0.3 MPa~1.2 MPa，压缩机出口压力5.5 MPa~6.7 MPa。

府谷压气站接收压力3.5 MPa~4.5 MPa。神府中区增压站进站压力0.7 MPa，集气站在不增压的情况下外输压力不大于0.9 MPa。

3.3. 算例一

此算例模拟了高需求工况下的管网运行优化方案，设定神府南区增压站进口压力为5.0 MPa，流量为11.8 × 10⁸ Nm³/a，府谷分输站进站压力为4.5 MPa，流量为0.2 × 10⁸ Nm³/a，黄甫1#集气站出站压力为5.95 MPa，流量为5 × 10⁸ Nm³/a，府谷1#集气站出站压力为5.90 MPa，流量为3.62 × 10⁸ Nm³/a，神府中区增压站出站压力为5.77 MPa，流量为3.45 × 10⁸ Nm³/a，以此作为可变参数的输入。优化结果包括各管道的输量、节点的压力、资源侧的外输量和分输站的供气量。该算例的输出结果如下表1~4所示，其中，输量的单位，万方/day。

Table 1. Output information of pipelines

表1. 管道输出信息

管道编号	输量	管输下限	管输上限	是否卡边	负荷率
1#管道	140.8571	0	142.8571	0	98.60%
2#管道	140.8571	0	142.8571	0	98.60%
3#管道	244.2857	0	285.7143	0	85.50%
4#管道	244.2857	0	285.7143	0	85.50%
5#管道	5.7143	0	114.2857	0	5.00%
6#管道	337.1429	0	371.4286	0	90.77%
7#管道	337.1429	0	371.4286	0	90.77%
8#管道	337.1429	0	371.4286	0	90.77%
9#管道	337.1429	0	371.4286	0	90.77%

Table 2. Output information on the resource side

表2. 资源侧输出信息

上载点	外输量	外输能力	外输量/外输能力
黄甫1#集气站	140.8571	142.8571	98.60%
府谷1#集气站	103.4286	103.4286	100%
神府中区集气站	98.5714	98.5714	100%

Table 3. Output information of nodes

表3. 节点输出信息

节点	压力(MPa)
1#节点(黄甫外输压力)	5.9500
2#节点	5.9267
3#节点(府谷外输压力)	5.9001
4#节点	5.8315
5#节点(神府中区外输压力)	5.7700
6#节点(府谷分输站进站压力)	5.7699
7#节点	5.5813
8#节点	5.4317
9#节点	5.2361
10#节点(神府南区增压站进站压力)	5.0763

Table 4. Distribution station output information

表4. 分输站输出信息

需求点	供气量	需求	用户满意度
府谷站外输需求	5.7143	5.7143	100.00%
神府南区外输需求	337.1429	337.1429	100.00%

3.4. 算例二

此算例模拟了低需求工况下的管网运行优化方案，设定神府南区增压站进口压力为5.0 MPa，流量为3.01 × 10⁸ Nm³/a，府谷分输站进站压力为4.5 MPa，流量为2.12 × 10⁸ Nm³/a，府谷1#集气站出站压力为5.06 MPa，流量为2.30 × 10⁸ Nm³/a，神府中区增压站出站压力为5.05 MPa，流量为3.05 × 10⁸ Nm³/a，以此作为可变参数的输入。优化结果包括各管道的输量、节点的压力、资源侧的外输量和分输站的供气量。该算例的输出结果如下表5~8所示，其中，输量的单位，万方/day。

Table 5. Output information of pipelines

表5. 管道输出信息

管道编号	输量	管输下限	管输上限	是否卡边	负荷率
1#管道	0.0000	0	142.8571	0	0.00%
2#管道	0.0000	0	142.8571	0	0.00%
3#管道	59.4286	0	285.7143	0	20.80%
4#管道	59.4286	0	285.7143	0	20.80%
5#管道	60.5714	0	114.2857	0	53.00%
6#管道	86.0000	0	371.4286	0	23.15%
7#管道	86.0000	0	371.4286	0	23.15%
8#管道	86.0000	0	371.4286	0	23.15%
9#管道	86.0000	0	371.4286	0	23.15%

Table 6. Output information on the resource side

表6. 资源侧输出信息

上载点	外输量	外输能力	外输量/外输能力
黄甫1#集气站	-	-	-
府谷1#集气站	59.4286	65.7143	90.43%
神府中区集气站	87.1429	87.1429	100%

Table 7. Output information of nodes

表7. 节点输出信息

节点	压力(MPa)
1#节点(黄甫外输压力)	—
2#节点	—
3#节点(府谷外输压力)	5.0600
4#节点	5.0549

续表

5#节点(神府中区外输压力)	5.0501
6#节点(府谷分输站进站压力)	5.0386
7#节点	5.0363
8#节点	5.0257
9#节点	5.0107
10#节点(神府南区增压站进站压力)	5.0000

Table 8. Distribution station output information

表8. 分输站输出信息

需求点	供气量	需求	用户满意度
府谷站外输需求	60.5714	60.5714	100.00%
神府南区外输需求	86.0000	86.0000	100.00%

4. 结论

本文针对天然气管网运行优化中存在的复杂约束与高效调度问题，提出了一种基于混合整数线性优化算法的优化方法，具体结论如下：

(1) 优化模型建立。通过建立以用户满意度最大化和管输能耗最小化为目标的数学模型，并结合分段线性化方法处理非线性水力约束，提升了模型的可解性和实用性。

(2) 算例分析验证。在调研基础上建立了管网运行优化模型，并对高需求与低需求两种典型工况进行了算例分析。结果表明：在高需求工况下，管网整体运行接近设计极限，各集气站均满负荷外输，系统虽紧张但仍能满足分输站需求；在低需求工况下，部分管道处于低负荷或闲置状态，黄甫集气站未启用，但府谷和神府中区仍能保证供气安全与用户满意度。研究结果说明，该优化模型能够准确反映不同工况下的运行特征与瓶颈，为后续管网的规划建设、调度优化及运行管理提供了有效的理论依据。

(3) 未来展望。随着天然气需求持续增长和管网运行环境日益复杂，未来研究可从以下几方面展开：一是进一步提升水力方程近似精度，减少分段线性化带来的误差；二是拓展至多工况和动态运行情境，增强模型的适应性和鲁棒性；三是引入不确定性建模与鲁棒优化，将气源波动、需求不确定性及价格变化纳入统一框架；四是结合大数据与智能算法，实现实时预测与自适应优化。上述方向将推动天然气管网优化方法向智能化、低碳化与可持续化发展，具有广阔的应用前景。

基金项目

重庆科技大学研究生科技创新项目“神府区块地面集输系统管网仿真模拟及运行优化”(YKJCX2420130)资助。