随着我国市场经济的快速发展，成品油的需求量也在大幅度增长，长距离、大流量、多卸油点的复杂成品油管道相继投产[1,2]。开泵方案的优劣直接影响管道运行与管理的安全性、经济性及稳定性[3,4]。为适应现场实际生产需求，输油泵机组存在串联、并联及混合连接等连接形式。针对水力系统复杂的成品油管道，时常辅助配置工作范围宽、适应性强的变频调速泵，增加了成品油管道开泵方案优化的难度[5,6]。

前人针对启泵方案优化问题进行了相关研究，目前最主要的研究方法为数学规划法。梁永图等[7]以管网运行周期内能耗费用最小为目标函数，建立了成品油管网优化运行数学模型，并采用分阶段动态规划方法进行模型求解。刘承婷等[8,9,10]以动态规划方法确定了最优启泵方案。以上方法对于长距离、多泵站管道求解速度较慢，且仅适用于配有定转速泵的优化问题。杨雪等[11,12,13,14]采用智能算法进行求解，在可接受的计算时间内得到了相对经济的开泵方案，但智能算法无法搜寻到全局最优解，且当问题规模扩大时收敛效果较差。

此外，也有学者通过建立数学模型求解特定问题。国外学者Cafaro等[15,16]通过采用泵特性方程内非线性项线性化处理的方式建立了混合整数线性规划(Mixed-Integer Nonlinear Programming，MILP)模型。Rejowski等[17]针对间歇输送模式，考虑了管道流量对摩阻损失及泵效率的影响，并采用连续时间变量建立混合整数非线性规划（Mixed Integer Nonlinear Programming，MINLP）模型。Bagirov[18]建立MINLP 模型，以二元变量表示泵启停状态，以连续变量表达泵运行时间，并开发出基于网格的组合搜索与Hooke-Jeeves 模式搜索方法的新算法用于模型求解。上述数学模型一般针对特定问题，不具有通用性，且与中国成品油管道的生产实际不够贴合。

以上研究都缺乏对泵机组连接形式以及泵调节方式的全面考虑，通常针对沿线泵站的定转速泵机组串联运行模式建立模型，无法实现对不同管道形式的统一求解，具有一定局限性。综上，基于已有批次输油计划，考虑泵机组并联运行模式以及变频泵的配置，以运行时间内管道全线运行能耗总和最小为目标建立MILP模型，用以求解成品油管道最优启泵方案。采用分段线性化方法实现非线性约束的线性化处理，并利用分支定界算法求解模型。

1 模型建立

1.1 模型假设

模型的假设条件：①管内油品为不可压缩流体；②不考虑混油长度，将相邻批次间视为混油界面；③不考虑压力和温度对油品密度及黏度造成的影响。

1.2 目标函数

模型基于离散时间表达[19]，以[Math Processing Error]T={1,2,...,tmax} 表示研究时间域内所有时间窗序号的集合，时间窗序号下标用tt 表示；以[Math Processing Error]I={1,2,...,imax} 表示管道沿线站场的编号，下标用ii 表示；以[Math Processing Error]K={1,2,...,Kmax} 表示站场ii 内所有输油泵编号的集合，下标用kk 表示（若第ii 站场非泵站，其[Math Processing Error]Kmax=0 ）。针对包含泵机组并联运行模式以及配备变频泵的成品油管道，以全线输油泵运行能耗最小为目标函数，其表达式为：

[Math Processing Error]minOF=∑t∑i∑kEPt,i,kτt （1）

式中：OFOF 为成品油管道全线输油泵总运行能耗，kW·h；[Math Processing Error]EPt,i,k 为第tt 时间窗内第i泵站的第k台泵的功率，kW；[Math Processing Error]τt 为第tt 时间窗的时长，h。

1.3 约束条件

在已知成品油管道批次输油计划的条件下，开泵方案优化约束主要分为6类：泵特性约束、泵转速约束、开泵数量约束、压力约束、流量约束及功率约束。

1.3.1泵特性约束

对于定转速泵机组以及额定转速下的变频泵机组，可采用最小二乘法拟合得到泵的流程-扬程特性曲线方程[20]。当变频泵转速发生改变时，其特性方程可由比例定律获得[21]。对于任意一台输油泵，其扬程需满足的约束条件为：

[Math Processing Error]Ht,i,k≤(1?BBPt,i,k)M+BVFi,k[(RVFt,i,kR0i,k)2(ai,k?bi,kQHt,i,k2?m)]+(1?BVFi,k)(ai,k?bi,kQHt,i,k2?m) （2）

[Math Processing Error]Ht,i,k≥(BBPt,i,k?1)M+BVFi,k[(RVFt,i,kR0i,k)2(ai,k?bi,kQHt,i,k2?m)]+(1?BVFi,k)(ai,k?bi,kQHt,i,k2?m) （3）

[Math Processing Error]?BBPt,i,kM≤Ht,i,k≤BBPt,i,kM （4）

式中：[Math Processing Error]Ht,i,k 为第t时间窗第i泵站的第k台泵提供的扬程，mm ；[Math Processing Error]BBPt,i,k 为表征输油泵启动状态的二元变量，[Math Processing Error]BBPt,i,k=1 表示第t时间窗第i泵站的第k台泵处于启动状态，[Math Processing Error]BBPt,i,k=0 表示处于关闭状态；[Math Processing Error]BVFi,k 为表征输油泵类型的二元参数，[Math Processing Error]BVFi,k=1 表示第i泵站的第k台泵为变频泵，[Math Processing Error]BVFi,k=0 表示定频泵；[Math Processing Error]RVFt,i,k 为第t时间窗第i泵站的第k台泵为变频泵时其转速值，r/minr/min ；[Math Processing Error]R0i,k 为当第i泵站的第k台泵为变频泵时其额定转速值，r/minr/min ；[Math Processing Error]QHt,i,k 为第tt 时间窗第i泵站的第k台泵的过泵流量，[Math Processing Error]m3/h ；[Math Processing Error]ai,k 、[Math Processing Error]bi,k 分别为第i泵站的第k台泵的流量-扬程特性曲线拟合参数；M为极大值。

1.3.2泵转速约束

对于变频泵，其转速一般存在范围限制，其约束条件为：

[Math Processing Error]?BBPt,i,kM≤BVFi,kRVFt,i,k≤BBPt,i,kM （5）

[Math Processing Error](BBPt,i,k?1)M+RVmini,k≤BVFi,kRVFt,i,k≤(1?BBPt,i,k)M+RVmaxi,k （6）

式中：[Math Processing Error]RVmini,k 为当第i泵站第k台泵为变频泵时其转速下限，r/minr/min ；[Math Processing Error]RVmaxi,k 为当第i泵站第k台泵为变频泵时其转速上限，r/minr/min 。

1.3.3开泵数量约束

每个站场泵运行的数量小于该站场配备输油泵总数，其表达式为：

[Math Processing Error]NBPt,i=∑kBBPt,i,k （7）

[Math Processing Error]NBPt,i≤Ki （8）

式中：[Math Processing Error]NBPt,i 为第t时间窗第i泵站的开泵数量，台；[Math Processing Error]Ki 为第i泵站配备的输油泵总数，台。

1.3.4压力约束

管道全线压力必须满足能量守恒定律，已知首站进站压力，随后各站进站压力等于前一站的出站压力与站间管段压力损失之差[22]，其表达式为：

[Math Processing Error]PINt,1=PG （9）

[Math Processing Error]PINt,i+1=POUTt,i?PFt,i （10）

[Math Processing Error]POUTt,i=PINt,i+rt,ig[BSPi∑kHt,i,k+(1?BSPi)HBt,i] （11）

式中：[Math Processing Error]PINt,i 、[Math Processing Error]POUTt,i 分别为第t时间窗第i泵站的进站和出站压力，MPa；[Math Processing Error]PG 为首站进站压力，MPa；[Math Processing Error]PFt,i 为第t时间窗第i泵站前管段的压力损失，包括摩擦阻力损失和高程损失，MPa；[Math Processing Error]rt,i 为第tt 时间窗第i泵站的过站油品密度，[Math Processing Error]kg/m3 ；g为重力加速度，[Math Processing Error]m/s2 ；[Math Processing Error]BSPi 为表征第i泵站输油泵连接关系的二元参数，[Math Processing Error]BSPi=1 表示串联，[Math Processing Error]BSPi=0 表示并联；[Math Processing Error]HBt,i 为第t时间窗第i泵站是并联泵站时其扬程，m。

公式（11）中并联泵站的扬程的约束条件为：

[Math Processing Error]?BBLt,iM≤HBt,i≤BBLt,iM （12）

[Math Processing Error](BBLt,i?1)M+Ht,i,1≤HBt,i≤(1?BBLt,i)M+Ht,i,1 （13）

[Math Processing Error]NBPt,iKi≤BBLt,i≤NBPt,i （14）

式中：[Math Processing Error]BBLt,i 为表征第tt 时间窗第i泵站是否开泵的二元变量，[Math Processing Error]BBLt,i=1 表示有泵处于开启状态，[Math Processing Error]BBLt,i=0 表示无泵开启。

管道沿线各站进、出站压力必须满足压力限制约束，其表达式为：

[Math Processing Error]PImini≤PINt,i≤PImaxi （15）

[Math Processing Error]POmini≤POUTt,i≤POmaxi （16）

式中：[Math Processing Error]PImini 、[Math Processing Error]PImaxi 分别为第i泵站的进站压力下限和上限，MPa；[Math Processing Error]POmini 、[Math Processing Error]POmaxi 分别为第i泵站的出站压力下限及上限，MPa。

管道中的高点需要满足高点压力下限约束，其表达式为：

[Math Processing Error]PHt,x≥PHmin （17）

式中：[Math Processing Error]PHt,x 为第t个时间窗高点x的压力，MPa；[Math Processing Error]PHmin 为全线压力低限，MPa。

管道中的低点压力不得超过该处管道最大承压，其表达式为：

[Math Processing Error]PLt,y≤PLmax （18）

式中：[Math Processing Error]PLt,y 为第t个时间窗低点y的压力，MPa；[Math Processing Error]PLmax 为该点处管道最大承压，MPa。

1.3.5流量约束

过泵流量与输油泵机组连接形式有关，需要满足范围约束。对于串联泵机组，过泵流量为该泵站的出站管段流量；对于并联泵机组，生产实际中常配置相同型号，其过泵流量为出站管段流量与当前启泵台数之比，其表达式为：

[Math Processing Error]?BBPt,i,kM≤QHt,i,k≤BBPt,i,kM （19）

[Math Processing Error](BBPt,i,k?1)M+QHmini,k≤QHt,i,k≤(1?BBPt,i,k)M+QHmaxi,k （20）

[Math Processing Error]QPt,i+(BSPiBBPt,i,k?1)M≤QHt,i,k≤(1?BSPiBBPt,i,k)M+QPt,i （21）

[Math Processing Error]QPt,iNBPt,i+[(1?BSPi)BBPt,i,k?1]M≤QHt,i,k≤[1?(1?BSPi)BBPt,i,k]M+QPt,iNBPt,i （22）

式中：[Math Processing Error]QHmini,k 、[Math Processing Error]QHmaxi,k 分别为第i泵站的第k台泵的过泵流量下限和上限，[Math Processing Error]m3/h ；[Math Processing Error]QPt,i 为第tt 时间窗第i泵站出站流量。

1.3.6 功率约束

输油泵运行功率与过泵流量、油品密度、扬程以及泵效率有关，其约束条件为：

[Math Processing Error]?BBPt,i,kM≤EPt,i,k≤BBPt,i,kM （23）

[Math Processing Error](BBPt,i,k?1)M+QHt,i,krt,igHt,i,kηPt,i,k≤EPt,i,k≤(1?BBPt,i,k)M+QHt,i,krt,igHt,i,kηPt,i,k （24）

式中：[Math Processing Error]ηPt,i,k 为第t时间窗第i泵站的第k台泵的效率。

2 模型求解

2.1 分段线性化

上述成品油管道开泵方案统一优化模型存在非线性约束，数学本质为MINLP模型，对其直接求解难度大，求解效率低[23]。因此，引入分段线性化策略，将MINLP模型转换为MILP模型进行求解，保证计算精度的同时，提高了求解效率[24]。

采用的分段线性化方法有两种（图1）：第一种采用一系列线性分段一次函数拟合连续非线性函数（图1a）；第二种采用一系列分段区间中值代替该区间内非线性函数值（图1b，m为列宾宗公式系数）。分段数可根据实际问题规模及函数非线性程度加以选取，所划区间越密集则求解精度越高，但同时问题规模也将急剧增大，因此须合理选取分段数。另外，区间不必等间隔划分，针对非线性程度较强的部分可适当加密区间。

图1 分段线性化方法示意图

以公式（2）为例进行说明，该式中变频泵转速平方项与流量幂次项呈现复杂的非线性耦合关系。针对转速平方项，将转速分为r个区间，采用第一种分段线性化方法，得到公式（25）。针对流量幂次项，将过泵流量分为s个区间，采用第二种分段线性化方法，得到公式（26），其表达式分别为：

[Math Processing Error]RVFt,i,k2|r≈A1rRVFt,i,k+B1r （25）

[Math Processing Error]QHt,i,k2?m|s≈(QHmins+QHmaxs2)2?m （26）

式中：[Math Processing Error]A1r 、[Math Processing Error]B1r 分别为第r个转速区间的一次函数表达式常系数；[Math Processing Error]QHmins 、[Math Processing Error]QHmaxs 分别为第s个流量区间的下限和上限；m为列宾宗公式系数，此处为0.25。

2.2 求解方法

分支定界算法以其灵活且便于使用计算机求解的优势成为解决MILP问题的重要方法。针对线性化后的MILP数学模型，采用以分支定界算法为基本求解原理的大规模数学规划求解器Gurobi进行求解，最终得到全局最优开泵方案。

3 实际算例

以中国某条成品油管道为研究对象，在管道批次注入顺序、初始时刻管道内油品状态以及中间各卸油站下载计划已知的基础上，可求解出研究时间域内管道整体的最优开泵方案。已知管道全长1 247.0 km，全线各类站场共16座，包括首站1座（以IS表示）、末站1座（以TS表示）、独立清管站1座、分输泵站3座（以P表示）及分输站10座（以D表示），管内顺序输送0#柴油、90#汽油、93#汽油。已知各分输泵站基本信息数据（表1），首站泵机组并联运行，其余3座分输泵站泵机组均串联运行，并通过配备的调速电机控制进出站压力，系统结构较为复杂。

表1 各分输泵站基本信息

已知该算例的输油计划批次运移图（图2）。其中，0#柴油以红色表示，90#汽油以绿色表示，93#汽油以蓝色表示。图中横坐标代表时间，左侧纵坐标轴代表管道的初始状态，右侧纵坐标代表站场体积坐标，黑线表示批次界面运移过程。横向矩形段表示各卸油站的卸油操作，其宽度表示下载流量大小，起、终点对应到横轴时间轴分别代表操作的开始、结束时刻。研究时间域长为190 h，初始状态管内油品顺序为90#汽油-0#柴油-93#汽油-90#汽油，0时刻开始首站依次输入0#柴油-90#汽油。

图2 输油计划批次运移图

利用上文提出的优化模型，基于批次运行计划对开泵方案进行优化。研究时间域分为15个时间窗，求解得到各时间窗起始时刻开泵方案，通过输油泵的启停及变频泵转速变化调节出站压力。使用Intel(R) Xeon(R) CPU E5-2640 v3@2.60GHz处理器电脑，并运用MATLAB R2016a对模型进行编程，得到各泵站优化开泵方案（表2）及各泵站进、出站压力曲线（图3）。分析可知：求解所得开泵方案下各站场在任意时刻进、出站压力均在允许范围内。首站配置并联定频泵，由批次运移图可知研究时间域内首站定流量注入，因此首站出站压力仅与过泵油品类型有关。在135 h时由于首站换油导致出站压力发生改变，与图3a相符。95.61 h至112.39 h之间，P3泵站不分输且末站下载流量增大，因此该时间窗内P3泵站和P4泵站需要提供更高的出站压力，增大变频泵转速，与图3c、图3d及表2相符。112.39 h后，P3泵站下载油品变更为93#汽油，其后管段水力损失降低，因此P3泵站内变频泵转速降低。综上，所得启泵方案可随全线水力工况变化进行相应调整，以降低管道运行能耗。

表2优化开泵方案结果 导出到

注：配泵方案中“11”代表该泵站内两种型号泵机组均处于开启状态；“10”代表该泵站内第一种型号泵机组开启，第二种型号泵机组关闭；“1”代表该站内泵机组处于开启状态。

图3 泵站进出站压力曲线

采用梁永图等[4]提出的动态规划算法求解，计算耗时为312.03 s，管道运行总能耗为8.121×105kW·h；采用上文的优化模型求解，求解耗时213.76 s，管道运行总能耗为7.708×105kW·h，相比动态规划法能耗降低约5%。可见，采用上文所提出的MILP优化模型进行求解能够降低管道运行能耗，具有较好的求解效果。

4 结论

（1）针对包含泵机组并联运行模式以及配备变频泵的成品油管道，考虑沿线节点压力约束、过泵流量约束等约束条件，基于已有批次输油计划，建立了以运行时间域内管道全线泵机组运行能耗总和最小为目标的成品油管道开泵方案统一优化MILP模型。

（2）采用分段线性化方法对模型中的非线性约束进行线性化处理，并利用分支定界算法对模型进行求解。在合理分阶段的情况下，求解效率较高，并可求得问题的全局最优解。

（3）该模型具有通用性，针对不同类型的泵机组连接模型以及输油泵调节形式的管道均可使用，通过调整初始参数设置即可求解，因此对其他实际运行的成品油管道开泵方案优化也具有指导意义。