韋金銀，郭琪，史彬，鄢烈祥

（武漢理工大學化學化工與生命科學學院，湖北武漢430070）

隨著經(jīng)濟的快速發(fā)展，我國成品油需求量逐年增長。成品油在配送過程中往往需要經(jīng)過醫(yī)院、學校、小區(qū)等人口密集地區(qū)，一旦發(fā)生燃燒爆炸安全事故，后果不堪設想[1]，因此，在進行成品油二次配送時要綜合考慮配送風險和配送成本。事實上，成品油二次配送已經(jīng)成為影響石化企業(yè)利潤的重要因素[2]，研究成品油二次配送路徑優(yōu)化問題，選擇合適配送路徑，降低成品油二次配送成本和風險、提高配送效率對石油石化企業(yè)具有重要的現(xiàn)實意義。

目前，國內(nèi)外學者對成品油二次配送路徑優(yōu)化問題已開展理論與應用研究，Cornillier 等[3]建立了帶時間窗的多車場成品油配送路徑優(yōu)化模型，運用啟發(fā)式算法求解，并運用數(shù)組實驗對算法的性能進行評估；Avella 等[4]研究了車輛有限、司機有限的成品油二次配送問題，結(jié)合分支定價算法，設計了一種快速求解的啟發(fā)式算法，并用實例驗證了算法的有效性；全自強等[5]針對成品油主動配送模式的新形勢，利用百度地圖API獲取基于實時路況下的行駛時間，建立數(shù)學模型運用禁忌搜索算法進行求解，優(yōu)化出來的配送路徑可以提高車輛配送周轉(zhuǎn)率；孫麗君等[6]建立了考慮司機工作量均衡的成品配送路徑多目標優(yōu)化模型，運用了改進的非支配排序遺傳算法進行求解，基于實例求解分析，驗證了算法的有效性，優(yōu)化出來的配送方案，可以為企業(yè)提供決策支持，對于實際有類似問題有一定啟發(fā)意義。

上述研究工作主要考慮了多車場、車輛與司機有限、實時路況、工作量均衡等因素進行研究，是屬于傳統(tǒng)的車輛路徑問題（vehicle route problem，VRP）一類的問題[7-8]。實際上，成品油屬于危險品，配送風險和成本都是影響成品油二次配送路徑規(guī)劃的重要因素，現(xiàn)有的研究鮮有考慮風險的因素，同時成品油二次配送風險和成本兩者相互沖突，進行多目標優(yōu)化研究具有重要意義。鑒于此，本文結(jié)合成品油配送過程的真實情況，考慮到配送風險和車速的動態(tài)變化，建立了考慮時變風險的成品油二次配送路徑優(yōu)化模型，提出改進的多目標列隊競爭算法進行求解，以期為石化企業(yè)的成品油安全高效配送提供借鑒。

1 道路時變風險度量模型

圖1 風險影響區(qū)域變化

成品油配送過程中，油罐車到站后卸油，載重量會減少，產(chǎn)生事故的影響隨之減少，故成品油在配送過程中的風險是動態(tài)變化的，其風險影響程度可以定量地用如圖1 所示的風險影響區(qū)域面積表示。

配送過程中，城市相同區(qū)域不同時間段的人群集聚程度不同，油罐車在不同時間段經(jīng)過相同區(qū)域時風險量的大小不同，考慮到載重和人群聚集程度的動態(tài)變化，可將風險度量模型表示為式(1)[9-10]。

式中，RTij表示路徑（i，j）上的Tij時段的風險量，次·人/km；pij表示路徑（i，j）上事故的發(fā)生概率，次/km；wij表示油罐車運輸過程中的載重量，t；wo表示油罐車最初的載重量，t；r表示路徑（i，j）上風險影響半徑，km；qTij表示路徑（i，j）上Tij時段的人口密度，人/km2；Dij表示路徑（i，j）上的距離，km。

2 時變?nèi)丝诿芏群凸?jié)點距離

2.1 時變?nèi)丝诿芏全@取

在上述風險度量模型中，人口密度是重要的參數(shù)之一，一天之內(nèi)不同區(qū)域人口密度會發(fā)生動態(tài)變化，要準確獲取實時人口密度難度較大。現(xiàn)有研究中是以隨機的方式生成人口密度[10]，這種方式產(chǎn)生的人口密度數(shù)值并不準確，相應地降低了時變風險計算的有效性。為了更準確地評估時變風險，本文提出了使用百度熱力圖定量計算人口密度，主要流程如下。

①收集一天百度熱力圖的數(shù)據(jù)進行分析，在7:00—20:00 每時間間隔為30min 截取一張熱力圖，總截取27張熱力圖。

②用熱力值（1～7）代表人群聚集程度，熱力值越大表示人群程度聚集程度越大，根據(jù)區(qū)域內(nèi)的熱力值，劃分為高熱區(qū)（6～7）和次熱區(qū)（4～5）[11]。

③可根據(jù)人群集聚程度對在7:00—20:00 這個時間段進行劃分，劃分的時間段為7:00—9:30、9:30—12:30、12:30—18:00 和18:00—8:00 這4 個時間段。

④把城市分成N 個區(qū)域，以城市平均人口密度為基礎，不同區(qū)域在t時間段的人口密度計算公式如式(2)所示。

2.2 節(jié)點距離獲取

百度地圖開放平臺為開發(fā)者提供程序接口，開發(fā)者可以根據(jù)自己的需求調(diào)用百度地圖程序接口獲取相關服務。本文運用百度地圖WEB 服務的批量算路服務，根據(jù)起點和終點，批量計算兩點之間的實際道路距離和耗時，Matlab環(huán)境中求兩個點之間的實際距離相關代碼如式(3)所示。

url='http://api.map.baidu.com/routematrix/v2/driving?output=json&origins=x1&destinations=x2&tactics=13&ak= 您的AK';

代碼中x1和x2分別代表起點和終點的坐標。兩點之間的來回距離并不相等，因為在實際道路中考慮了單雙行道、是否禁止左右拐等因素，這樣的距離更具有實際意義。

3 BP神經(jīng)網(wǎng)絡預測車速

城市交通錯綜復雜，車輛在相同路段不同時間段有不同的車速，車速的變化直接影響到配送風險和配送成本的度量，因此需要進行車速的預測。本文以歷史車速數(shù)據(jù)為輸入，建立BP 神經(jīng)網(wǎng)絡，預測各個時段車速，此網(wǎng)絡隱含層神經(jīng)元的輸出如式(4)所示。

輸出層的輸出如式(5)所示。

式中，i ∈n，n表示輸入層單元個數(shù)；j ∈o，o表示隱含層的單元個數(shù)；k ∈h，h表示輸出層的個數(shù)；βij和βkj為連接權值；θj和θk為神經(jīng)元閥值；xi為輸入的歷史車速，km/h；bj為隱含層輸出值；yh為預測的車速，km/h；f(x)為傳遞函數(shù)。

4 配送模型和求解

4.1 問題描述

成品油二次配送主要是采用公路運輸方式將成品油從油庫運輸?shù)郊佑驼镜倪^程，其成品油配送路徑優(yōu)化問題描述如下。①每個加油站的需求量、油庫到各個加油站的距離和各個加油站之間的距離已知。②油庫的庫存足夠大，可以滿足各個客戶的需求，每輛車的載重量一定，都是同一種車型。③車的載重量大于任何一個加油站的需求量，一個加油站只能由一輛車服務，一輛車可以服務多個加油站。④車從油庫出發(fā)，經(jīng)過幾個加油站完成配送后回到油庫，而且每輛車有自己的服務時間窗，車最晚回到油庫的時間不能超過時間窗。⑤每個客戶只能被服務一次。

在滿足上述條件下，同步考慮配送成本和配送風險進行成品油二次配送路徑優(yōu)化。

4.2 雙目標優(yōu)化模型

考慮時變風險成品油二次配送優(yōu)化模型見式(6)和式(7)。

約束條件見式(8)～(15)。

式(6)和式(7)表示最小化配送成本和風險；式(8)表示每輛車的配送任務不能超過車的最大載重量；式(9)表示車輛從油庫出發(fā)最后又回到油庫；式(10)和式(11)表示每個加油站只能被一輛車服務；式(12)表示車輛最晚回到油庫的時間不能超過T；式(13)表示車輛從加油站i 到加油站j 的時間關系；式(14)表示車輛早到和晚到的懲罰關系式；式(15)表示車輛在運輸過程中載重量的變化關系式。

4.3 車速時變性分析

根據(jù)路網(wǎng)把路段分成高速路段和主干路段，按照車輛一天速度的變化分成n個時段，各個時段表示為[T1,T2,…,Tn]，車輛在路段(i,j)車速有兩種情況：①車輛通過路段(i,j)時，恰好在第Tm個時間段內(nèi)，那么車輛在路段(i,j)的車速只有一種車速為②車輛通過路段(i,j)跨過了2 個及以上時間段時，則車輛的行車速度具有兩個以上的車速，車輛在第m 個時間段內(nèi)先以車速行駛，依此類推，到第Tm+n時段后以車速vTm+nij行駛，相應行駛距離為。綜上，車輛從加油站i 到加油站j 的時間可表示為式(16)。

根據(jù)一天內(nèi)人群的群聚集程度變化情況，劃分為m 個時段，各個時段表示為[Tt1,Tt2,…,Ttm]，車速的變化引起tij的變化，會導致車輛在路段(i,j)經(jīng)過不同的時段,同上，車輛在路段(i,j)上Tij時段的人口密度可表示為式(17)。

4.4 算法設計

列隊競爭算法（LCA）是一種并行搜索、多層競爭的全局優(yōu)化搜索算法，已經(jīng)成功應用于求解流程工業(yè)間歇生產(chǎn)調(diào)度和化工過程模擬等多種優(yōu)化問題[12-14]?？紤]時變風險成品油二次配送路徑優(yōu)化問題是雙目標優(yōu)化問題，針對這類問題使用多目標列隊競爭算法（MOLCA）求解，可以快速得到分布均勻的Pareto 最優(yōu)解[15-16]。本文在MOLCA 基礎上，引入非支配排序和擁擠度排序來分配家族搜索空間，在選擇父代家族無性繁殖后代時，只選擇非支配解，如果非支配解低于初始種群數(shù)，再隨機生成新的家族進行補充，這種方式可以很好地避免陷入局部最優(yōu)解，并且可以得到更多的Pareto 最優(yōu)解。算法基本步驟如下：①隨機生成m組實數(shù)序列代表m 個初始家族；②計算m 個家族的配送成本和風險；③進行非支配排序，處于同一非支配排序等級時，進行擁擠度計算，擁擠度距離較大的家族排在隊列的前面；④根據(jù)各個家族在列隊中的位置，處于前面的家族有較大的概率進行小規(guī)模的變異，處于末位的家族有較大概率進行無性繁殖；⑤計算子代的配送成本和風險；⑥進行非支配排序，選擇非支配解作為父代進行下一代繁殖，如果非支配解沒有達到初始家族數(shù)，再隨機生成新的家族進行補充；⑦判斷當前的迭代次數(shù)是否大于最大迭代次數(shù)，如果是輸出最優(yōu)解，否則轉(zhuǎn)②繼續(xù)迭代。

4.5 算法的變異和繁殖策略

本文采用的是實數(shù)型編碼，變異策略是在數(shù)列中隨機產(chǎn)生兩個點，并交換兩個點的位置，繁殖策略采取的是無性繁殖，在數(shù)列中隨機產(chǎn)生兩個點，并在這兩個點這間進行隨機排列，具體的過程如圖2所示。

圖2 變異繁殖操作示意圖

5 實例計算與結(jié)果分析

本文以武漢某石化企業(yè)2018年12月11日的成品油配送路徑優(yōu)化問題為例，以石化企業(yè)的油庫為配送中心，取16 個加油站為該油庫的配送客戶，所選取的加油站如圖3所示?？蛻艚邮艿姆諘r間s=30min，客戶時間窗在[0,660]內(nèi)取值，其中0代表時刻8:00，660代表時刻19:00，時間窗轉(zhuǎn)換的方式是以時刻8:00為基準，某一時刻相對于基準過去的時間，比如時刻9:10 相對時刻8:00 過去了70min，那么70 代表的是時刻9:10，油罐車的最大載重量為8t。實例數(shù)據(jù)如表1所示。

5.1 BP神經(jīng)網(wǎng)絡車速預測結(jié)果分析

圖3 加油站分布圖

表1 實例數(shù)據(jù)

在成品油二次運輸過程中，大多時間都是在高速路/快速路和主干路上行駛，所以本文主要考慮油罐車在運輸過程中高速路/快速路和主干路上的行車速度，假設次干路速度和主干路行車速度相等。BP 神經(jīng)網(wǎng)絡傳遞函數(shù)選擇雙極S 型函數(shù)，對于訓練函數(shù)選取學習率可變的最速下降BP 算法，網(wǎng)絡的迭代次數(shù)設置為1000，目標誤差設置為0.01，不同數(shù)據(jù)的最佳網(wǎng)絡結(jié)構(gòu)不同，所以對于主干路車速的預測，輸入層節(jié)點數(shù)為5，隱含層節(jié)點數(shù)為7，輸出層節(jié)點數(shù)為1 的網(wǎng)絡結(jié)構(gòu)，對于高速路/快速路的預測，輸入層節(jié)點數(shù)為5，隱含層節(jié)點數(shù)為14，輸出層節(jié)點數(shù)為1的網(wǎng)絡結(jié)構(gòu)。數(shù)據(jù)收集選擇從2018 年10 月29 日到12 月14 日工作日的35天數(shù)據(jù)，每天從7:00 到20:00 每30min 收集一次高速路/快速度和主干路的平均車速，樣本數(shù)據(jù)為10月29日到12月9日工作日的數(shù)據(jù)，12月10日到12月14日為預測數(shù)據(jù)。本文以12月11日周二的預測為例，對預測結(jié)果進行分析，采用MATLAB 軟件對數(shù)據(jù)進行處理以及分析，抽取周二組的數(shù)據(jù)對網(wǎng)絡進行訓練，高速路/快速路和主干路網(wǎng)絡的R 值分別為0.95596、0.97337，可見預測精度較高。為了進一步評估模型，引入平均相對誤差（MRE），見式(18)。

式中，yi為實際車速，km/h；y'為預測的車速數(shù)據(jù)，km/h。通過計算快速路的MRE 值為7.99%，主干路的MRE 值為5.68%，預測程度具有較高的可信度，兩者的速度對比圖如圖4所示。

圖4 預測車速和實際速度對比圖

城市的交通雖然每時每刻都在變化，但是有的時候在一個時間范圍內(nèi)車輛速度的變化并不大，根據(jù)武漢交通委員會交通周報，早高峰時間段為7:00—9:00 和晚高峰時間段為17:00—19:00，以及收集到的數(shù)據(jù)進行分析，把一天的時間分為7:30—9:00、9:00—10:30、10:30—12:30、12:30—14:30、14:30—17:00、17:00—19:00 和19:00—20:00 這7 個時間段，劃分后時間段的車速預測結(jié)果如表2所示。

表2 車速預測結(jié)果

5.2 實例計算結(jié)果分析

模型參數(shù)設計，單位運輸成本C 取值為1 元，啟動單個車輛成本B取值為200元，早到懲罰費用E取值為10元，晚到懲罰費用L取值為40元，事故發(fā)生概率pij取值為每106km 發(fā)生7.75 次[17]，風險半徑r取值為1.6km[1]，考慮到大多加油站分布在武漢市主城區(qū)，故-Q的取值為武漢市主城區(qū)平均人口密度6362 人/km2，wt本人選擇在1～3 范圍內(nèi)進行取值，在matlab2009 環(huán)境中分別用MOLCA 和改進的多目標列隊競爭算法（I-MOLCA）求解。求得的Pareto最優(yōu)解如圖5所示。

由圖5 可以看出，改進的MOLCA 的Pareto 最優(yōu)解完全支配了MOLCA 的Pareto 最優(yōu)解，且得到更多的Pareto 最優(yōu)解。Pareto 最優(yōu)解集只是為決策者提供了多種決策方案，為了選擇合適的解，決策者可以從最優(yōu)解集中選取一個妥協(xié)解，妥協(xié)解可以根據(jù)經(jīng)驗、實際情況或某種規(guī)則進行選取，本文妥協(xié)解取Pareto 最優(yōu)解中與理想解的相對距離最小解[18]，相應的函數(shù)表達式如式(19)所示。

圖5 Pareto最優(yōu)解

式中，PCS為妥協(xié)指數(shù)，(W,E)為Pareto 最優(yōu)解的任一最優(yōu)解，(W?,E?)為理想解，其中W*和E?分別表示以成本為單目標和以風險為單目標的最優(yōu)解。求得的妥協(xié)解配送成本1280.2，配送風險為37.3，相應路線如表3所示。

表3 妥協(xié)解的相應路線

為了驗證算法的有效性，本文用I-MOLCA 求解文獻[19]的案例，得到結(jié)果對比如表4 所示，IMOLCA求得的非支配解，完全支配文獻[19]的兩個非支配解。

從表3 看出一共6 條路徑，但是只用了2 輛油罐車。由于受到載重的限制，油罐車在完成一條線路配送任務后返回配送中心，還在油罐車的配送時間窗內(nèi)，如果停止配送會導致部分油罐車閑置，且現(xiàn)實情境中，多啟動一輛油罐車的費用是比較高的，考慮到上訴因素本文允許在油罐車的配送時間窗內(nèi)可以進行多次配送。

表4 I-MOLCA求解文獻結(jié)果對比

本文通過求解不考慮時變風險的模型，求得Pareto 最優(yōu)解和考慮時變風險的Pareto 最優(yōu)解進行比較，如圖6。由圖6 可知，圖中不考慮時變風險的曲線的坡度相對于考慮時變風險的曲線小了很多；求得的妥協(xié)解的路徑比較由圖7 和圖8 可知，不考慮時變風險求得的路線是7條，考慮時變風險求得的路線是6條，其中不同顏色代表不同的行車路徑，兩邊優(yōu)化的路線差異是較大的；綜上兩種比較結(jié)果表明動態(tài)變化的人群聚齊程度對路線優(yōu)化的結(jié)果影響是比較大的。

圖6 Pareto最優(yōu)解對比圖

圖7 考慮時變風險的路徑優(yōu)化路線圖

6 結(jié)論

圖8 不考慮時變風險的路徑優(yōu)化路線圖

考慮到成品油二次配送的特性，采用BP 神經(jīng)網(wǎng)絡預測車速，以配送過程的風險和成本為目標，提出了改進的多目標列隊競爭算法，用于求解成品油二次配送路徑優(yōu)化雙目標優(yōu)化問題，并通過計算某武漢石化企業(yè)和文獻的實例，獲得以下結(jié)論：①預測的高速路段和主干路段的車速與實際車速的平均相對誤差分別為7.99%和5.68%，預測結(jié)果具有較高的預測精度；②與文獻差分進化算法求得的結(jié)果進行比較，本文提出的多目標列隊競爭算法求得的結(jié)果優(yōu)于文獻求得的結(jié)果；③求解得到了同步考慮配送過程的時變風險和成本的成品油二次配送路徑優(yōu)化的解集，企業(yè)可以根據(jù)實際情況進行選擇。

符號說明