天然气管网规划中的混合整数规划模型与算法

Mixed-Integer Programming Models and Algorithms in Natural Gas Network Planning

摘要: 天然气管网规划问题极具挑战性，其核心是混合整数非线性规划(MINLP)特征。该特征主要源于复杂的网络拓扑、非线性的物理守恒关系，以及离散的管道建设与设备选型决策。本文针对该MINLP模型，采用分段线性近似技术，将其转化为可求解的混合整数线性规划(MILP)模型。本文为了加速MILP模型的求解效率，深入分析了管网的特殊拓扑结构和流量守恒特性，并基于此提出了强化的有效不等式(即加速割平面)。该不等式经验证可有效地应用于一般拓扑结构下的管网规划MILP模型。数值实验结果表明，所提出的不等式能够显著加强模型的线性松弛界，从而有效缩短MILP模型的求解时间。

Abstract: The natural gas network planning problem presents significant challenges, primarily characterized by its Mixed-Integer Nonlinear Programming (MINLP) nature. This complexity stems from intricate network topology, nonlinear physical conservation relationships, and the discrete decisions involved in pipeline construction and equipment selection. This paper addresses the MINLP model by employing a Piecewise Linearization technique to reformulate it into a solvable Mixed-Integer Linear Programming (MILP) model. To accelerate the solution efficiency of the resulting MILP model, we conduct an in-depth analysis of the pipeline network’s special topological structure and flow conservation properties. Based on this analysis, we propose a strengthened valid inequality (i.e., accelerating cutting plane). The inequality is demonstrated to be effectively applicable to MILP models for pipeline network planning under general topological structures. Numerical results confirm that the proposed inequality significantly tightens the linear relaxation bound of the model, consequently effectively reducing the MILP solution time.

1. 引言

天然气管网的规划布局与运行优化一直是能源领域的研究焦点和行业关注的热点。国际上，学者们围绕不同目标展开了深入探索。Misra等[1]专注于在维持管网运行压力的前提下实现运行成本的最小化，并对大型传输系统进行了分析。Alves等[2]则构建了管输收益最大化模型。Bazyar等[3]通过Torabi-Hassini方法，将经济、环境和社会目标集成到多元优化模型中，并采用多阶段方法求解。在国内，研究也取得了显著进展。陈志华[4]以天然气市场的社会福利最大化为目标，建立了长输管网的流量优化配置模型，并指出管网互联互通不足是制约效率的关键因素。

然而，伴随市场对天然气需求的持续增大，现有管网设施的输送能力已难以满足现阶段的需求。尤其是在我国天然气输送管道持续扩张并形成大型互联管网的背景下，如何协调气源、用户节点与现有拓扑结构，优化配置管网的流量流向，确定新建管线的选址与容量，成为管网设施建设中亟待解决的关键规划问题。该问题实质上是一个涉及连续变量和离散变量的混合整数非线性规划问题，求解难度极大。

针对该MINLP模型求解的挑战性，本文利用其结构特征，采用分段线性近似技术将其转化为可求解的混合整数线性规划(MILP)模型。此外，本文结合管网问题的物理意义和拓扑特性，提出了相应强化的加速割平面。数值结果表明，该方法在保证求解精度的前提下，能够显著提升求解效率，为天然气管网的规划决策提供了高效且可靠的优化工具。

2. 背景

天然气管网规划问题属于复杂的混合整数非线性规划(MINLP)范畴，即目标函数或约束中包含非线性项，且部分决策变量必须取整数值。MINLP问题的求解难度远高于线性规划(LP)和混合整数线性规划(MILP)，尤其当问题包含非凸函数时，寻找全局最优解变得极其困难[5] [6]。

2.1. MINLP问题求解算法的分类

在过去的几十年中，研究者们开发了多种针对MINLP的精确算法，这些方法主要根据其处理问题的凸性特征进行分类。

处理凸MINLP的经典算法，适用于目标函数和所有非线性约束函数均为凸函数的情况。

分支定界法(Branch-and-Bound, B&B)：这是最基础的框架，通过不断分支整数变量和求解连续松弛问题来收敛到全局最优解[5] [7]。

外逼近法(Outer-Approximation, OA)：由Duran和Grossmann提出，通过交替求解一个非线性规划(NLP)子问题和一个MILP主问题来迭代逼近最优解，在凸MINLP领域被广泛应用[6] [8]。

广义Benders分解法(Generalized Benders Decomposition, GBD)：适用于目标函数和约束中的非线性部分仅涉及连续变量的情况[5] [7]。

扩展割平面法(Extended Cutting Plane, ECP)：该方法主要适用于凸MINLP问题[5] [7]，其显著特征是不依赖于求解非线性规划(NLP)子问题，而是通过迭代，添加基于非线性约束和目标函数的线性化割平面，持续强化MILP松弛模型，直至收敛到最优解。

处理非凸MINLP的全局优化算法：天然气管网中的水力学约束(如流量与压降关系)通常是非凸的，使得问题更加复杂。主要的策略是构建原始问题的凸松弛模型，结合空间分支定界(Spatial B&B)框架进行全局搜索[6] [9]。该方法的核心在于开发有效的上/下估计函数(Under/Over-estimators)来紧致松弛界。例如通过使用多项式优化技术或半定规划(SDP)松弛等方法[6]。

2.2. 基于线性化和割平面的求解策略

求解MINLP的另一种主流且高效的策略是将其近似或精确地转化为MILP模型，从而利用成熟的MILP求解器(如Gurobi、CPLEX、SCIP、HIGHS、COPT)进行求解。对于复杂的非线性约束，分段线性近似(PLA)技术将非线性函数替换为一组线性约束和辅助整数变量，通过在函数的定义域内引入多个断点并使用特殊集(Special Ordered Sets, SOS)或凸组合模型(如$\lambda$ 模型)进行重构，可以将非凸或非线性的函数精确或近似地转化为MILP形式[10]。Vielma等人[10]对不同的分段线性近似建模技术及其对MILP求解性能的影响进行了深入研究和比较。这种方法牺牲了精度(在近似情况下)但极大地提高了求解效率，尤其适用于工程规划问题。本文即采用这种策略，将非线性的天然气规划模型转化为MILP模型进行求解。本文即采用这种策略，将非线性的天然气规划模型转化为MILP模型进行求解。

无论是由MINLP转化而来的MILP，还是原始的MILP模型，其求解效率的关键都在于线性规划(LP)松弛界的紧致程度。割平面方法是分支定界法中最重要的加速技术之一，通过迭代地向LP松弛中添加有效不等式(即割平面)，可以削去非整数解，同时不排除任何整数可行解，从而显著提升下界，减少搜索树的规模。针对特定问题的结构化割平面往往更有效。文献[7]中提到了用于求解凸MINLP的割平面方法。近年来，研究趋势集中于利用问题的特殊结构(如网络结构、组合结构)来推导出有效不等式。

然而，目前已有的有效不等式基本都只适用于特定的问题结构，为了考虑将有效不等式应用到本文所考虑的天然气管网规划问题中，我们需要根据问题的特殊结构进行构造。

本文的创新点正是结合天然气管网独特的网络拓扑结构和流量守恒特性，提出了有效不等式，以期克服PLA转化后模型规模大、LP松弛界较弱的问题，从而加速求解过程。

3. 分段线性近似与有效不等式推导

3.1. 问题建模

给定一个拓扑图$G\left(V,E\right)$ 表示天然气运输管网拓扑，$V$ 图中所有点的集合，表示需求点和供应点；$E$ 表示图中所有边的集合，表示运输天然气的管道。此处特别指出，为了将重点放在线性化模型与有效不等式的推导上面，对模型中参数的单位不做过多介绍。图中每条边上面有三个属性，分别是单价$u_{ij}^s,ij\in E,s\in \left\{1,2,3,4\right\}$ ，表示通过边$ij,ij\in E$ 运输的天然气处于区间$\left[a^s,a^{s+1}\right],s\in \left\{1,2,3,4\right\}$ 内的运输单价。长度$L_{ij},ij\in E$ 表示边的长度。容量上下限约束$\left[lo{w}_{ij},u{p}_{ij}\right]$ ，表示通过这条边的运输量必须处于区间$\left[lo{w}_{ij},u{p}_{ij}\right]$ 内。同时考虑到随着市场对于天然气需求量的增加，原始天然气管道总运输能力不足以满足市场总需求量，因此需要考虑在原始边$E$ 的位置新建管线。由此，为了使运输费用最小可以构建这样一个初始MINLP模型。

$\text{min}{\displaystyle \sum }_{i,j}{\displaystyle \sum }_s{u}_{ij}^s\ast L_{ij}\ast x_{ij}^s+{\displaystyle \sum }_{i,j}{\displaystyle \sum }_s{u}_{ij}^s\ast L_{ij}\ast y_{ij}^s+C\ast {\displaystyle \sum }_{i,j}{L}_{ij}\ast {\left(y_{ij}\right)}^p$

$\begin{array}{c}s.t.{\displaystyle \sum }_k\left(x_{ki}+y_{ki}\right)-{\displaystyle \sum }_j\left(x_{ij}+y_{ij}\right)=b_i,i\in V\end{array}$ (1)

$\begin{array}{c}lo{w}_{ij}\le \left|x_{ij}-x_{ji}\right|\le u{p}_{ij},ij\in E\end{array}$ (2)

$\begin{array}{c}{z}_{ij}^s\ast a^s\le x_{ij}^s\le z_{ij}^s\ast a^{s+1},s\in \left\{1,2,3,4\right\}\end{array}$ (3)

$\begin{array}{c}z{\text{'}}_{ij}^s\ast a^s\le y_{ij}^s\le z{\text{'}}_{ij}^s\ast a^{s+1},s\in \left\{1,2,3,4\right\}\end{array}$ (4)

$\begin{array}{c}{\displaystyle \sum }_s{z}_{ij}^s=1,{\displaystyle \sum }_{s}z{\text{'}}_{ij}^s=1,z_{ij}^s,z{\text{'}}_{ij}^s\in \left\{0,1\right\}\end{array}$ (5)

$\begin{array}{c}{x}_{ij}={\displaystyle \sum }_s{x}_{ij}^s,y_{ij}={\displaystyle \sum }_s{y}_{ij}^s\end{array}$ (6)

$\begin{array}{c}0\le x_{ij},y_{ij}\end{array}$ (7)

参数说明，此处仅对当前模型中出现的变量结合实际进行介绍，具体可参考后文的变量说明表：

$x_{ij}$ ：从节点i到j的流量。

$y_{ij}$ ：从节点i到j的扩容流量。

$C$ ：扩容成本常系数。

$b_i$ ：流守恒给定的供给量或者需求量。

$x_{ij}^s$ ：从节点i到j处于区间$\left\{0\right\}\cup \left[a^s,a^{s+1}\right]$ 的流量。

$y_{ij}^s$ ：从节点i到j处于区间$\left\{0\right\}\cup \left[a^s,a^{s+1}\right]$ 的扩容流量。

$z_{ij}^s$ ：取值为1时，表示$x_{ij}^s\in \left[a^s,a^{s+1}\right]$ ；取值为0时，表示$x_{ij}^s=0$ 。

$z{\text{'}}_{ij}^s$ ：取值为1时，表示$y_{ij}^s\in \left[a^s,a^{s+1}\right]$ ；取值为0时，表示$y_{ij}^s=0$ 。

$p$ ：非线性项指数，$0<p<1$ 。

决策变量：

$x_{ij}$ ，$y_{ij}$ ，$x_{ij}^s$ ，$y_{ij}^s$ ，$z_{ij}^s$ ，$z{\text{'}}_{ij}^s$ 。其中$x_{ij}$ ，$y_{ij}$ ，$x_{ij}^s$ ，$y_{ij}^s$ 为连续变量，$z_{ij}^s$ ，$z{\text{'}}_{ij}^s$ 为0-1二元变量。

模型实际含义如下：

目标函数：极小化运输成本。成本计算由三部分组成，${\displaystyle \sum }_{i,j}{\displaystyle \sum }_s{u}_{ij}^s*L_{ij}*x_{ij}^s$ 表示通过已有管线进行天然气运输的总成本。${\displaystyle \sum }_{i,j}{\displaystyle \sum }_s{u}_{ij}^s*L_{ij}*y_{ij}^s$ 表示通过新建管线进行运输的总成本。这两项都是线性函数。$C*{\displaystyle \sum }_{i,j}{L}_{ij}*{\left(y_{ij}\right)}^p$ 表示新建管线的建设成本是非线性函数，且由于包含${\left(y_{ij}\right)}^{0.244}$ ，是一个凹函数。因此，需要求的问题是一个凹的混合整数非线性规划问题，求解难度极大，即使是具有特殊性的MINLP也通常被证明是NP-hard问题[6] [10]。

约束(1)：流守恒约束，满足各节点的供给量或者需求量约束。

约束(2)：对应拓扑图中边的容量约束。

约束(3)到(6)：通过0-1变量限制流量$x_{ij}$ ，扩容流量$y_{ij}$ 的取值区间。

3.2. 分段线性近似模型

首先模型中的绝对值约束(2)可以通过引入0-1变量将其线性化等效于约束(8) (9) (10)。

$\begin{array}{c}lo{w}_{ij}\ast {\alpha }_{ij}\le x_{ij}\le u{p}_{ij}\ast {\alpha }_{ij}\end{array}$ (8)

$\begin{array}{c}lo{w}_{ij}\ast {\beta }_{ij}\le x_{ji}\le u{p}_{ij}\ast {\beta }_{ij}\end{array}$ (9)

$\begin{array}{c}{\alpha }_{ij}+{\beta }_{ij}\le 1, {\alpha }_{ij},{\beta }_{ij}\in \left\{0,1\right\}\end{array}$ (10)

引入0-1变量将原始模型中的绝对值约束改写为线性约束之后，整个模型仅目标函数中的最后一项包含${\left(y_{ij}\right)}^p$ 是非线性的，而所有约束都是线性的。

由于原始目标函数是对${\left(y_{ij}\right)}^p$ 的非线性项求和，可以将每一个包含${\left(y_{ij}\right)}^p$ 的求和项单独进行分段线性近似如下：

${\left(y_{ij}\right)}^p\approx {\displaystyle \sum }_m{e}^m*y_{ij}{}^m+f^m*t_{ij}^m$

$\begin{array}{c}{\displaystyle \sum }_m{t}_{ij}^m=1,t_{ij}^m\in \left\{0,1\right\}\end{array}$ (11)

$\begin{array}{c}{l}_m\ast t_{ij}^m\le y_{ij}{}^m\le r_m\ast t_{ij}^m\end{array}$ (12)

$\begin{array}{c}{y}_{ij}={\displaystyle \sum }_s{y}_{ij}^s={\displaystyle \sum }_m{y}_{ij}^m\end{array}$ (13)

线性化引入新的0-1变量$t_{ij}^m$ 表示在分段区间$\left[l_m,r_m\right]$ 内使用$e^m\ast y_{ij}{}^m+f^m\ast t_{ij}^m$ 替代原始${\left(y_{ij}\right)}^p$ 的取值。$e^m,f^m$ 为常数，分别表示区间$\left[l_m,r_m\right]$ 内线性近似函数的斜率，截距。通过不断增大$m$ 的取值，增加分段数来拟合原始非线性项，能够不断地逼近原始目标函数，但是随着$m$ 的增大，0-1变量数逐渐增加，模型求解效率急剧下降。通过对原始模型进行线性化，可以得到新的MILP模型目标函数。

$\text{min} {\displaystyle \sum }_{i,j}{\displaystyle \sum }_s{u}_{ij}^s\ast L_{ij}\ast x_{ij}^s+{\displaystyle \sum }_{i,j}{\displaystyle \sum }_s{u}_{ij}^s*L_{ij}*y_{ij}^s+C\ast {\displaystyle \sum }_{i,j}\left(L_{ij}\ast {\displaystyle \sum }_m{e}^m\ast y_{ij}{}^m+f^m\ast t_{ij}^m\right)$

约束条件由(1) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13)组成。至此，MINLP问题被转换成了MILP问题。

3.3. 有效不等式推导

构造有效不等式：

${\displaystyle \sum }_{ij\in E^{\prime }}{y}_{ij}\ge all\_req-maxflow$

$all\_req$ 表示已知的总需求量，$maxflow$ 表示通过图论中经典算法——最小费用最大流算法(MCMF)计算出的原始拓扑网络在不新建或扩容管线的情况下所能输送的天然气总量。

$E^{\prime }$ 表示通过MCMF算法确定的所有有效流量路径中，流向指向需求点的管线集合。

不等式构造逻辑：我们约束集合$E^{\prime }$ 中所有管线的总扩容量不小于总需求量$all\_req$ 与最大输送总量 $maxflow$ 的差值$all\_req-maxflow$ 。该不等式的引入旨在保证模型具有可行解空间，并且不会排除任何一个最优可行解(即该不等式是有效的)。

为了更好地理解不等式中$all\_req$ ，$maxflow$ ，$E^{\prime }$ 的实际意义。消除$E^{\prime }$ 的模糊性，避免在不同拓扑上构造不等式产生歧义。

通过下面两幅图对比，从实际应用角度出发能更加容易理解其实际含义。

图1中每条边上的第一个数字表示网络流通过边的花费，第二个数字表示边的最大容量。粉色节点A表示供应点，黑色节点B，C表示中间节点需求量为0，蓝色节点D，E表示需求点。假设需要将A的总供应量15全部运往D，E，表示此示例中$all\_req=15$ 。通过MCMF算法可以得到最大流的网络流路径，如下图所示。

Figure 1. Initial Network

图1. 原始网络

Figure 2. MCMF Network

图2. 最小费用最大流网络

图2中红色边为最大流经过的边，并且通过箭头表明了流向。边上的绿色数字表示通过边的流量。统计发现，目前的网络最大运输能力为12，表示此示例中$maxflow=12$ 。因此需要新建边，按照前文对于新建边集合$E^{\prime }$ 的描述：$E^{\prime }$ 表示通过MCMF算法确定的所有有效流量路径中，流向指向需求点的管线集合。也就是与需求点D，E相连，并且有流量流向D，E的边的集合。由此可以得到在此示例中$E^{\prime }$ 表示为边BD，BE，CE的集合。综上所述，详细阐述了如何准确获取不等式中的各项参数，构造有效不等式。通过将不等式添加到MILP模型中可以加速问题求解。

3.4. 变量说明表

前文建模过程中与数学模型中所提到的所有变量符号含义说明如表1所示。

Table 1. Variables Description

表1. 变量说明表

变量符号	含义
$G$	给定拓扑图
$V$	拓扑图中所有节点集合
$E$	拓扑图中所有边集合
$i$	拓扑图中节点，$i\in V$
$j$	拓扑图中节点，$j\in V$
$s$	单价变化分段数，$s\in \left\{1,2,3,4\right\}$
$u_{ij}^s$	边$ij$ 运输单价，取值由$s$ 决定
$L_{ij}$	边$ij$ 长度，1000~2000之间的随机值
$C$	新建管线成本计算常系数
$p$	非线性项指数，$0<p<1$ 随机值
$b_i$	节点$i$ 的供应量$b_i<0,$ 需求量$b_i\ge 0$
$lo{w}_{ij}$	原始管线运输量下界约束
$u{p}_{ij}$	原始管线运输量上界约束
$x_{ij}$	从$i->j$ 的原始管线运输量，连续变量
$y_{ij}$	从$i->j$ 的新建管线运输量，连续变量
$x_{ij}^s$	$i->j$ 原始管线运输量，由$z_{ij}^s$ 决定，连续变量
$y_{ij}^s$	$i->j$ 新建管线运输量，由${z^{\prime }}_{ij}^s$ 决定，连续变量
$z_{ij}^s$	0-1变量，决定$x_{ij}^s$ 取值区间
${z^{\prime }}_{ij}^s$	0-1变量，决定$y_{ij}^s$ 取值区间
$a^s$	不同取值区间端点值
${\alpha }_{ij}$	0-1变量，取1表示流向$i->j$
${\beta }_{ij}$	0-1变量，取1表示流向$j->i$
$m$	分段线性近似的分段数，取2
$e^m$	分段线性近似的线性函数斜率
$f^m$	分段线性近似的线性函数截距
$t_{ij}^m$	0-1变量，决定$y_{ij}^m$ 的取值区间
$y_{ij}^m$	分段线性近似函数的近似变量，连续变量
$l_m$	分段区间的左端点
$r_m$	分段区间的右端点
$E^{\prime }$	有效不等式边集
$all\_req$	总需求量
$maxflow$	原始网络最大运输量

4. 数值实验

本次数值实验在Windows 11操作系统环境下进行，硬件配置为CPU R5-3550H和16 GB RAM。我们采用Python 3.10编程语言，并结合Pulp库构建MILP模型。模型求解由SCIP(9.1.0)完成。实验设计旨在对比两种情景下的计算效率：基准模型(不添加有效不等式)与强化模型(添加本文提出的有效不等式)。通过比较两种模型设置下的求解时间以验证所提出的有效不等式在加速MILP求解方面的有效性。

4.1. 测试数据

随机生成50，60，70，80，90节点的网络拓扑图，边数为节点数的3倍。供应节点数为节点数的1/10，需求节点数处于节点数的3/10~8/10之间，剩余节点为中间点。保证供应点的总供应量大于总需求量；保证原始随机网络的最大流小于总需求量。每个节点规模生成20组测试数据。分段设置在100，150，200，分成两段[100, 150]、[150, 200]进行线性近似。最终需要求解的0-1变量数为边数的12倍。由于0-1变量数对问题求解效率的影响极大，而不同网络规模的拓扑中，0-1变量数最终可以通过节点总数的固定倍数表示，因此，后续在实验结果中通过节点规模来标识实验组别。

4.2. 数值结果

表2呈现了不同网络规模下基准模型(未添加不等式)与强化模型(添加本文提出的有效不等式)的平均求解时间对比。其中，节点规模代表随机生成网络中的节点总数。数值实验结果表明，本文提出的有效不等式能够显著提升MILP模型的求解效率。在所有测试规模下，强化模型的平均求解时间均低于基准模型。特别地，在节点规模为50的情况下，求解时间的缩短效果最为显著，求解时间从172秒降至133秒，加速比约为22.7% (即求解时间降低到原始的77.3%)。整体而言，该不等式能够在不同网络规模下将平均求解时间缩短3.8%至22.7%。这有力地证明了该不等式组能够有效加强模型的线性松弛界，从而显著加速一般拓扑结构下管网规划MILP模型的求解过程。

Table 2. Average solution times (Unit: seconds)

表2. 平均求解时间表(单位：秒)

节点规模	50	60	70	80	90
基准模型	172	292	285	401	738
强化模型	133	281	270	341	656

5. 结论

本文聚焦于天然气管网规划这一具有挑战性的决策问题，该问题被建模为复杂的混合整数非线性规划模型，其难度主要源于非线性的物理守恒关系和离散的规划决策(新建管线选址与容量)。本文主要针对天然气管网规划的MINLP模型，本文成功应用了分段线性近似(Piecewise Linearization)技术，将其高效地转化为可由常数求解器处理的混合整数线性规划模型。这一转化策略为大规模管网规划提供了可行的求解路径。为克服线性近似模型可能存在的线性松弛界较弱、求解效率低的问题，本文深入分析了管网的特殊拓扑结构和流量守恒特性，并基于此提出了强化的有效不等式(即加速割平面)，与现有的有效不等式相比，本文提出的有效不等式充分结合了天然气运输网络的特性，更加具有针对性。数值实验结果有力地证明了该加速割平面的有效性。通过对比基准模型(未添加不等式)和强化模型(添加不等式)的平均求解时间，结果显示，所提出的不等式能够显著加强模型的线性松弛界，从而有效缩短MILP模型的求解时间。在本次实验数据中，求解时间最多可缩短22.7%，证明了所提不等式在提升管网规划求解效率方面的价值。

综上，本文的方法在保证规划精度的前提下，极大地提高了天然气管网规划问题的求解效率，为能源基础设施的优化决策提供了高效且可靠的工具。