成品油管道分輸計劃優(yōu)化的并行模擬退火算法

陳海宏 ，左麗麗，吳長春，李清平

1 中國石油大學(xué)(北京)油氣管道輸送安全國家工程實(shí)驗室/城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗室，北京 102249

2 中海油研究總院深水工程國家重點(diǎn)實(shí)驗室，北京 100028

成品油管道運(yùn)行與管理的重要內(nèi)容之一是制定批次計劃。可行的批次計劃需滿足站場分輸/注入流量、站場分輸/注入條件、管段運(yùn)行流量、批量和批次次序等約束[1]。對于中石油現(xiàn)有成品油管道，借助于批次計劃模擬軟件，計劃編制人員手工制定一個月的管道批次計劃通常需耗時半天以上。為提高計劃編制人員的工作效率，許多研究人員提出了批次計劃優(yōu)化模型與算法，并開發(fā)了相應(yīng)的優(yōu)化軟件[2-9]，明顯縮短了計劃編制時間。針對成品油管道批次計劃優(yōu)化問題，文獻(xiàn)[10～13]采用離散時間方程建立了混合整數(shù)線性規(guī)劃(MILP)模型，而文獻(xiàn)[14～19]采用連續(xù)時間方程建立了批次計劃優(yōu)化的MILP模型。與離散模型相比，連續(xù)模型的規(guī)模更小，計算時間更短，優(yōu)化結(jié)果質(zhì)量更高[20-22]。當(dāng)問題規(guī)模較小時，通過CPLEX商業(yè)優(yōu)化軟件一般可在較短時間內(nèi)獲取MILP模型最優(yōu)解。隨著問題規(guī)模的擴(kuò)大，CPLEX軟件的計算時間將大幅增加，其往往只能獲取MILP模型較優(yōu)解，有時甚至無法獲取可行解，故CPLEX不適合于管道規(guī)模大、計劃周期長、批次數(shù)目多、油品種類多的情況。為保證在較短時間內(nèi)產(chǎn)生較優(yōu)解，文獻(xiàn)[11]提出采用模擬退火(SA)算法優(yōu)化單點(diǎn)注入、多點(diǎn)分輸成品油管道分輸計劃。然而，對于大規(guī)模問題，SA算法應(yīng)用效果依舊不佳。為改善SA算法的優(yōu)化效果，實(shí)現(xiàn)SA算法并行計算，本文提出構(gòu)造分輸計劃新方案的兩階段法，并改進(jìn)了目標(biāo)函數(shù)。

以單點(diǎn)注入、多點(diǎn)分輸成品油管道為研究對象，已知管道沿線各站場里程、各管段內(nèi)徑、各分輸站分輸流量上/下限值、各管段運(yùn)行流量上/下限值、計劃期內(nèi)管道起點(diǎn)輸入各批次的時間與流量、計劃期開始時刻管道內(nèi)各批次的油品種類、體積與位置、計劃期內(nèi)各中間分輸站每種油品的分輸需求體積，在滿足管道運(yùn)行約束的前提下，優(yōu)化管道沿線各中間站分輸計劃，包括站場分輸各批次的開始時刻、結(jié)束時刻、流量與體積。

1 并行模擬退火算法

1.1 算法總體框架

成品油管道分輸計劃優(yōu)化的并行SA算法邏輯框架(圖1)主要包含8個要素：目標(biāo)函數(shù)、初始解構(gòu)造方法、新解構(gòu)造方法、新解接收準(zhǔn)則、初始溫度Tmax、終了溫度Tmin、降溫速率r、Markov鏈長度S。

(1)分輸計劃初始方案構(gòu)造方法采用空間遞推法(Space Recursive Method, SRM)[11]。按照SRM法，從上游站場至下游站場，依次基于經(jīng)驗規(guī)則編制各分輸站分輸計劃初始方案。

(2)分輸計劃新方案構(gòu)造過程分為兩個階段：第一階段利用變量鄰域解構(gòu)造方法局部調(diào)整分輸計劃原方案；第二階段基于分輸作業(yè)銜接性規(guī)則微調(diào)第一階段產(chǎn)生的分輸計劃新方案。

(3)在每次迭代過程中，同時分別構(gòu)造多個新的分輸計劃，并選擇一個最好的新解作為本次迭代的最終解，從而實(shí)現(xiàn)SA算法并行計算。若同時存在多個最好的新解，則從中隨機(jī)選擇一個。

(4)分輸計劃新方案接收準(zhǔn)則為Metropolis準(zhǔn)則：若新解優(yōu)于舊解，則接收新解；否則，按概率接收新解。概率計算公式為：

其中，prof為差解被接收的概率；fn為新解目標(biāo)函數(shù)值；fo為舊解目標(biāo)函數(shù)值；T為SA算法的當(dāng)前溫度。

1.2 目標(biāo)函數(shù)

以計劃期內(nèi)管道沿線所有管段運(yùn)行流量波動幅度之和最小作為目標(biāo)函數(shù)。

其中，JS為站場總數(shù)；KS為時間節(jié)點(diǎn)總數(shù)；g1為計劃期內(nèi)各管段運(yùn)行流量總體波動幅度之和，m3/h；g2為計劃期內(nèi)各管段運(yùn)行流量瞬時波動幅度之和，m3/h；cpj為管段j權(quán)重系數(shù)，用于避免運(yùn)行流量的波動集中在某些管徑較大、里程較長的管段上，建議取值為管段j體積容量與管道總體積容量的比值；csj,k為0、1變量，csj,k=1 表示時段(tk-1，tk)和(tk，tk+1)內(nèi)管段j運(yùn)行流量不相等，csj,k=0則表示運(yùn)行流量相等；QDj,k為時段(tk，tk+1)內(nèi)管段j運(yùn)行流量，m3/h；QAj表示計劃期內(nèi)管段j運(yùn)行流量的平均值，m3/h。計劃期劃分方法：選取計劃期內(nèi)管道沿線各站場執(zhí)行分輸/注入作業(yè)的開始時刻與結(jié)束時刻作為時間節(jié)點(diǎn)，將計劃期劃分為若干時段。

1.3 分輸計劃初始方案構(gòu)造方法

為降低成品油管道分輸計劃優(yōu)化問題的求解難度，假設(shè)計劃期內(nèi)任一中間站對同一批次至多執(zhí)行一次分輸作業(yè)。實(shí)際上，目前中石油成品油管道分輸計劃的編制基本上符合該假設(shè)。基于該假設(shè)，采用SRM法[11]構(gòu)造分輸計劃初始方案。按照SRM法，從上游站場至下游站場，依次基于經(jīng)驗規(guī)則構(gòu)造管道沿線各分輸站分輸計劃初始方案。確定某站場是否分輸某批次的原則：計劃期內(nèi)該批次可以到達(dá)該站場、該站場對該批次所含油品存在分輸需求。確定某站場分輸各批次體積的原則：按計劃期內(nèi)批次過站體積的比例將該站場分輸每種油品的需求體積分配至各個批次。確定某站場分輸某批次流量的原則：優(yōu)先選取分輸流量允許范圍內(nèi)的中間值。確定站場分輸某批次時間的原則：為保證初始方案有較大的調(diào)整空間，將每次分輸作業(yè)安排在批次過站時間窗的中間時段。

1.4 分輸計劃新方案構(gòu)造方法

分輸計劃新方案構(gòu)造過程(圖2)可分為兩個階段。第一階段為采用變量鄰域解構(gòu)造方法調(diào)整分輸計劃原方案；第二階段為基于分輸作業(yè)銜接性規(guī)則微調(diào)第一階段產(chǎn)生的分輸計劃新方案。

圖1 成品油管道分輸計劃優(yōu)化的并行SA算法示意圖Fig. 1 Schematic diagram for optimizing delivery schedules of a multiproduct pipeline by parallel SA algorithm

(1)基于變量鄰域解構(gòu)造方法調(diào)整分輸計劃原方案

A 隨機(jī)選擇一個中間分輸站sj。

B 隨機(jī)選擇一個批次si，其油品類型為sp。

C 若計劃期內(nèi)被選中站場未分輸被選中批次，則增加一個新的分輸作業(yè)；否則，從分輸時間、分輸流量、分輸體積3個變量中隨機(jī)選擇一個變量，利用變量鄰域解構(gòu)造公式調(diào)整被選中分輸作業(yè)的執(zhí)行計劃。變量鄰域解構(gòu)造公式[23]如下：

其中，γ為決策變量，Xa≤γ≤Xb；γ’表示γ的鄰域解；T表示SA算法的當(dāng)前溫度；σ為0～1隨機(jī)數(shù)；sgn為符號函數(shù)。

(2)基于分輸作業(yè)銜接性規(guī)則微調(diào)分輸計劃

為提高分輸計劃新方案的質(zhì)量，第一階段基于變量鄰域解構(gòu)造法調(diào)整被選中分輸作業(yè)后，嘗試將該分輸作業(yè)與其它中間站執(zhí)行的分輸作業(yè)進(jìn)行銜接。若第一階段調(diào)整的分輸作業(yè)或者新增分輸作業(yè)的執(zhí)行開始/結(jié)束時刻與相鄰某站場所執(zhí)行分輸作業(yè)的開始/結(jié)束時刻的時間間隔小于φ小時，則銜接這兩個分輸作業(yè)。φ一般可取5 h。

圖2 分輸計劃新方案構(gòu)造流程Fig. 2 Flowchart for generating a new delivery schedule

1.5 參數(shù)設(shè)置

根據(jù)SA算法理論可知，算法迭代總次數(shù)由外循環(huán)與內(nèi)循環(huán)共同控制。初始溫度、終止溫度與降溫函數(shù)可用于控制外循環(huán)執(zhí)行次數(shù)，Markov鏈長度可用于控制內(nèi)循環(huán)執(zhí)行次數(shù)。根據(jù)文獻(xiàn)[11]中SA算法參數(shù)設(shè)置方法進(jìn)行參數(shù)設(shè)置，具體如下：

初始溫度：Tmax= 1000 ℃。

終止溫度：Tmin= 1 ℃。

Markov鏈長度：S =1000。

降溫函數(shù)：T’=r×T。其中，r表示降溫速率，0＜r＜1。需根據(jù)具體測試情況，綜合比較優(yōu)化結(jié)果的質(zhì)量與計算耗時確定r最佳取值。

并行線程數(shù)PN：需根據(jù)具體測試情況，綜合比較優(yōu)化結(jié)果的質(zhì)量與計算耗時確定PN最佳取值。

2 算例分析

以某成品油管道PA為例，驗證并行SA算法優(yōu)化管道分輸計劃的有效性。管道PA總長2000 km、內(nèi)徑544.8 mm、沿線共建有8個站場(編號分別為E0～E7)；管道沿線各管段運(yùn)行流量允許范圍為400～1200 m3/h，沿線各中間分輸站分輸流量范圍見表1。管道PA測試的計劃周期為0～670 h，已知計劃期內(nèi)管道起點(diǎn)輸入計劃(表2)、計劃期內(nèi)管道沿線各中間站每種油品分輸需求體積(表3)和初始時刻(0 h)管道內(nèi)各批次的油品種類、體積及位置(表4)，利用并行SA算法優(yōu)化管道PA分輸計劃。

2.1 冷卻系數(shù)r的影響

影響SA算法穩(wěn)定性的最主要因素是冷卻系數(shù)r。一般來說，r越大，SA算法優(yōu)化結(jié)果質(zhì)量越高，優(yōu)化耗時越長，算法穩(wěn)定性也越高。根據(jù)具體管道測試情況，綜合優(yōu)化結(jié)果質(zhì)量與耗時可確定r推薦值。不考慮并行計算，設(shè)定不同的r值，應(yīng)用SA算法優(yōu)化管道PA分輸計劃(圖3)。當(dāng)0.05≤r≤0.85時，隨著r增加，優(yōu)化結(jié)果質(zhì)量呈明顯上升趨勢。當(dāng)0.85≤r≤0.95時，隨著r增加，優(yōu)化結(jié)果質(zhì)量提升幅度較小，但計算耗時顯著增加。因此，將r推薦值設(shè)定為0.85。

表1 管道PA沿線中間分輸站分輸流量允許范圍Table 1 Delivering flow rate limits of intermediate stations along pipeline PA

表2 計劃期內(nèi)管道PA起點(diǎn)輸入計劃Table 2 Inputting schedules for the initial station of pipeline PA during the scheduling horizon

表3 計劃期內(nèi)管道PA中間分輸站對各油品的分輸需求體積Table 3 Product demands of intermediate stations along pipeline PA during the scheduling horizon

表4 初始時刻管道PA內(nèi)各批次油品的種類、體積及位置Table 4 Product type, volume and location of every batch in pipeline PA at the initial moment

2.2 并行線程數(shù)PN的影響

令r=0.85，設(shè)定不同的并行線程數(shù)目PN，應(yīng)用并行SA算法優(yōu)化管道PA的分輸計劃(圖4)。

(1)當(dāng)1≤PN≤8時，隨著PN增加，并行SA算法計算耗時呈線性增長。并行SA算法采用的并行策略為：每次迭代時，同時產(chǎn)生若干個新解，并選擇質(zhì)量最高的解作為本次迭代的最終解。總迭代次數(shù)不變，但每次迭代時由于增加了各線程之間的等待耗時，使得并行SA算法的計算耗時隨PN增加而增大。當(dāng)PN=8時，計算時間約為6 min，滿足單條管道分輸計劃編制工作允許計算時間的要求。

(2)當(dāng)1≤PN≤6時，隨著PN增加，優(yōu)化結(jié)果質(zhì)量呈明顯上升趨勢，目標(biāo)函數(shù)值由5010 m3/h降低至4490 m3/h。當(dāng)6≤PN≤8時，隨著PN增加，優(yōu)化結(jié)果質(zhì)量提升幅度較小，目標(biāo)函數(shù)值由4490 m3/h降低至4480 m3/h，但計算耗時卻由325 s增加至386 s，增長約19%。

綜合考慮優(yōu)化結(jié)果質(zhì)量與計算耗時，將PN推薦值設(shè)定為6。與PN=1相比，當(dāng)PN=6時，優(yōu)化結(jié)果目標(biāo)函數(shù)值降低約10%。

圖3 不同冷卻系數(shù)下，SA算法優(yōu)化管道PA分輸計劃的結(jié)果Fig. 3 The optimization results for pipeline PA provided by SA algorithm with different cooling factors

圖4 當(dāng)r=0.85時，不同并行線程數(shù)下，并行SA算法優(yōu)化管道PA分輸計劃的結(jié)果Fig. 4 The optimization results for pipeline PA provided by parallel SA algorithm with different parallel threads when r=0.85

2.3 優(yōu)化結(jié)果展示

令r=0.85、PN=6，20次應(yīng)用并行SA算法優(yōu)化管道PA分輸計劃，平均耗時約5.4 min，優(yōu)化結(jié)果目標(biāo)函數(shù)平均值為4490 m3/h。圖5為某次利用并行SA算法產(chǎn)生的管道PA分輸計劃，優(yōu)化結(jié)果目標(biāo)函數(shù)值為4386 m3/h。圖5中P1、P2、P3分別表示0#柴油、92#汽油、95#汽油。計劃期內(nèi)，中間分輸站E1—E2、E3—E4、E5—E6之間所執(zhí)行的分輸計劃具有很好的銜接性，保證了管道沿線各管段流量運(yùn)行平穩(wěn)。

2.4 并行SA算法與CPLEX對比

文獻(xiàn)[10]采用離散時間方程建立了成品油管道分輸計劃優(yōu)化的MILP模型，并利用CPLEX求解模型。基于文獻(xiàn)[10]建立的MILP模型，利用CPLEX優(yōu)化管道PA的分輸計劃。實(shí)際生產(chǎn)中，單條管道分輸計劃編制工作允許的計算時間一般不超過10 min。將計算時間設(shè)定為10 min，CPLEX成功產(chǎn)生了管道PA的分輸計劃，優(yōu)化結(jié)果目標(biāo)函數(shù)值為6021 m3/h。CPLEX與并行SA算法的計算過程如圖6所示。與CPLEX相比，并行SA算法優(yōu)化結(jié)果的目標(biāo)函數(shù)值降低約27%，計算耗時降低約46%。從優(yōu)化結(jié)果來看，并行SA算法明顯優(yōu)于CPLEX。

圖5 當(dāng)r=0.85、PN=6時，應(yīng)用并行SA算法產(chǎn)生的管道PA分輸計劃Fig. 5 The delivery schedules of pipeline PA provided by parallel SA algorithm, when r=0.85 and PN=6

圖6 CPLEX、并行SA算法優(yōu)化管道PA分輸計劃的迭代過程Fig. 6 Computational processes about optimizing delivery schedules of pipeline PA using CPLEX and parallel SA algorithm

3 結(jié)束語

通過在SA算法中引入并行計算，建立了適用于成品油管道中間站分輸計劃優(yōu)化的并行SA算法，其綜合性能不僅優(yōu)于非并行SA算法，而且優(yōu)于CPLEX商業(yè)優(yōu)化求解器。不同于非并行SA算法，并行SA算法最大的特點(diǎn)是：從當(dāng)前迭代點(diǎn)出發(fā)同時構(gòu)造并檢驗多個鄰域解，選擇其中最好者作為本次迭代產(chǎn)生的新迭代點(diǎn)。對于某成品油管道PA，將并行計算線程數(shù)取為6個，并行SA算法可以保證在計劃編制工作可接受的計算時間內(nèi)(不超過10 min)產(chǎn)生較優(yōu)的成品油管道分輸計劃。基于并行SA算法產(chǎn)生的管道PA分輸計劃具有很好的銜接性，符合計劃編制人員編制批次計劃的要求，保證了管道沿線各管段流量運(yùn)行平穩(wěn)。與非并行SA算法相比，并行SA算法的優(yōu)化結(jié)果目標(biāo)函數(shù)值降低約10%。與CPLEX優(yōu)化求解器相比，在10 min的計算時間內(nèi)，并行SA算法的優(yōu)化結(jié)果的目標(biāo)函數(shù)值降低約27%。