侯 艷，黃康煥，張億仙，伍乃騏

(廣東工業(yè)大學 1.計算機學院；2.機電工程學院，廣東 廣州 510006)

原油一次加工過程具有生產規(guī)模大、約束多、過程連續(xù)、裝置多和工藝復雜等特點，是一個包含離散事件和連續(xù)過程的混合系統(tǒng)，其短期生產計劃和調度問題屬于NP-hard問題[1]。短期生產計劃是煉廠根據市場需求和國際油價變化等因素制定的7~10天生產計劃。目前仍然缺乏成熟的商用軟件和技術用于制定原油處理過程短期生產調度計劃，大部分工作需要決策者手工完成。而優(yōu)化的調度計劃能夠顯著提高煉廠的利潤。因此，如何制定優(yōu)化的調度計劃成為研究的熱點之一。

部分學者通過數學規(guī)劃方法對調度問題進行求解。Wu等[2]通過線性規(guī)劃方法對輸油管道的輸油過程進行建模，降低管道轉運原油過程中消耗的能量；Hamisu等[3]通過將煉油調度周期等分為若干個時間間隔，即離散時間表示，對原油處理過程建立混合整數線性規(guī)劃模型，降低油輪卸油成本、油輪在海上的等待成本、油罐庫存成本和切換成本；為了保證求解的精度，離散時間表示需要將調度周期等分為盡可能小的時間片，導致產生大量的二元決策變量?，F有軟件不能很好地求解大規(guī)模問題，而連續(xù)時間表示能夠有效減少決策變量的數量。Zhang等[4]考慮管道的長距離運輸，以原油卸載成本、庫存成本、原油混合成本等為優(yōu)化目標，采用連續(xù)時間表示，建立混合整數非線性規(guī)劃模型。陳璇[5]對沿海型和內陸型原油調度問題建立連續(xù)時間模型，降低油輪海上等待成本等操作成本。連續(xù)時間表示雖然能有效減少決策變量的數量，但是容易產生非線性約束，導致求解非常困難，在實際應用中不可行。另外，數學規(guī)劃方法用于求解調度問題時，要求短期生產計劃已知。然而在制定詳細調度前，短期生產計劃未知，因此，數學規(guī)劃方法不能直接應用于原油短期生產調度。

為了解決這個問題，Wu等[6]采用Petri網對原油處理過程進行建模和分析，將生產調度問題劃分為上下2個層次。上層采用線性規(guī)劃得到蒸餾塔煉油計劃以最大化生產速率；下層制定原油詳細調度計劃以實現上層的煉油計劃。文獻[7-8]分別采用啟發(fā)式算法和混合整數規(guī)劃模型對原油調度問題建模，求解優(yōu)化的調度計劃。文獻[9-12]根據兩層模型和控制理論，采用Petri網對系統(tǒng)的可調度性進行分析，得到各種情況下的可調度性條件。然而據此得到的詳細調度計劃不一定是一個最優(yōu)的調度計劃。

為了制定優(yōu)化的詳細調度計劃，Hou等[13]將原油短期調度問題轉換為供油罐和蒸餾塔的指派問題，以系統(tǒng)的可調度性條件為安全邊界，通過指派蒸餾塔和供油罐執(zhí)行調度決策，啟發(fā)式地尋找原油的詳細調度計劃以實現上層的煉油計劃。同時，以供油罐使用個數成本、罐底混合成本、管道混合成本以及蒸餾塔的供油罐切換成本為優(yōu)化目標，通過采用NSGA-II算法[14]尋找問題的Pareto最優(yōu)解。

原油處理過程短期生產調度問題是多目標組合優(yōu)化問題。進化算法常用于求解多目標優(yōu)化問題。陳祥等[15]采用文化基因算法求解開放車間調度問題；Jiang等[16]采用粒子群算法求解云制造供應鏈調度的多目標優(yōu)化問題。對于高維多目標進化算法，多樣性和收斂性是評價算法性能的2個重要指標。為了提高進化算法的多樣性和收斂性性能，本文在文獻[13]已經將原油短期調度問題轉換為供油罐到蒸餾塔的指派問題的基礎上，提出基于改進的骨干粒子群算法(improved bare bone particle swarm optimization, I-BBPSO)和NSGA-II算 法協(xié) 同進化 的雙種群進化算法(bare bone particle swarm optimization and NSGA-II, BPGA)對原油短期調度問題進行求解。

1 問題描述和數學模型

1.1 問題描述

煉廠的運作過程有3個階段：原油一次加工過程、生產過程和成品油轉運。圖1為原油一次加工過程示意圖。首先，油輪將原油卸載到儲油罐中；接著，輸油管道將港區(qū)儲油罐中的原油轉運到廠區(qū)供油罐中；最后將供油罐中的原油供給蒸餾塔進行分餾。原油一次加工后得到的石油半成品經過后續(xù)的混合加工后，最終得到需要的石油產品，并將其運輸到市場。

圖 1 原油一次加工過程Figure 1 Crude oil operations

原油經過卸載、轉運和供油3個階段后，產生一系列的成本。

1) 油輪在港口等待卸油時，根據停泊時間長短需要向港口繳納相應的費用；

2) 油輪向不同的港區(qū)儲油罐卸載原油時，產生切換成本；

3) 輸油管道轉運原油時產生能耗成本；

4) 輸油管道轉運原油時，管道中不同種類的原油在交匯處產生混合成本；

5) 供油罐向外放完油后，不能擠出罐中的所有原油，罐底還剩余一部分原油；當該油罐進行下一次充油時，不同種類原油就會在罐底混合，混合后的原油品質下降，產生原油的混合成本；

6) 蒸餾塔在煉油期間，需要切換多個供油罐為其供油，產生切換成本；

7) 原油一次加工過程產生供油罐使用成本。

蒸餾塔煉油計劃為蒸餾塔DSk具體的煉油時間、分餾的原油種類，以及相應的原油體積。本文基于煉廠已有的蒸餾塔煉油計劃，假設港區(qū)儲油罐有足夠的原油可供調度，輸油管道以最大的輸油速率fpmax將原油從港區(qū)儲油罐轉運到廠區(qū)供油罐，制定優(yōu)化多個目標的原油詳細調度計劃，以實現蒸餾塔的煉油計劃。因此，本文目標是優(yōu)化上面提到的成本(4~7)。

1.2 數學模型

1.2.1 符號說明

Λ為待處理的原油類型總量；

K為蒸餾塔的數量；

?為供油罐個數；

OD為原油調度決策；

COT為原油類型；

S為原油輸入的來源；

D為原油輸出的目的地；

[a,b]為 OD的起始時刻和結束時刻；

V為 OD轉運的原油體積；

Θ={1, 2,···, Λ}，為待處理的原油類型集合；

DS={1, 2,···,K}，為蒸餾塔集合；

TK={1,2,···, ?}，為供油罐集合；

TKr為供油罐r；

DSk為蒸餾塔k；

Φ為原油在供油罐中的駐留時間；

τs、 τe分別為蒸餾塔煉油計劃的開始時刻和結束時刻；

dij為i號原油和j號原油在管道的混合成本系數；

αij為i號原油和j號原油在管道中的混合次數；

cij為i號原油和j號原油在供油罐中的罐底混合成本系數；

βij為i號原油和j號原油在供油罐中的罐底混合次數；

σ為一個供油罐的使用成本系數；

φ為原油一次加工過程使用的供油罐總個數；

ω為蒸餾塔切換一次供油罐的成本系數；

γk為DSk切換供油罐的次數；

Cr為TKr的容量；

vr(?)為? 時刻TKr中的原油體積；

Qr,COT(?)為? 時刻TKr中原油種類的數量；

Sk(?)∈{0,1}， ?時刻，1表示DSk正在煉油，0表示DSk沒有煉油；

QTKk(?)為? 時刻向DSk供油的供油罐數量；

Sr,in(?)∈{0,1}， ?時刻，1表示TKr正在充油，0表示TKr沒有充油；

Sr,out(?)∈{0,1}， ?時刻，1表示TKr正在放油，0表示TKr沒有放油。

1.2.2 問題數學模型

五元組OD=(COT,V,S,D,[a,b])表示原油調度決策。根據原油輸入的來源地和輸出的目的地，OD分為原油卸載決策ODU(S為油輪，D為港區(qū)儲油罐)、原油轉運決策ODT(S為港區(qū)儲油罐，D為廠區(qū)供油罐)以及供油決策ODF(S為廠區(qū)供油罐，D為蒸餾塔)。原油的詳細調度計劃需要保證煉油計劃在時間域[τs,τe]能夠完成。

目標函數為

其中，f1為管道中原油的混合成本；f2為供油罐罐底原油的混合成本；f3為使用供油罐的總成本；f4為蒸餾塔切換供油罐產生的成本。

約束條件為

式(1)表示相鄰2次ODT的原油種類相同時，管道中原油混合成本系數di j為0；式(2)表示TKr前后2次裝載原油的種類相同時，原油混合成本系數cij為0；式(3)表示任意時刻 ?，TKr中的原油體積滿足供油罐容量約束；式(4)表示任意時刻，TKr中只能裝有一種原油；式(5)表示煉油任務完成之前所有蒸餾塔持續(xù)煉油；式(6)表示任意時刻，至少有一個供油罐為DSk供油；本文中，在任意時刻，設為只有一個供油罐向DSk供油；式(7)保證了任意時刻TKr不能同時充油和放油；另外，供油罐充油完畢后，必須等待一段時間才能向外放油，以分離出原油中的海水，這段時間稱為駐留時間約束。

2 可調度性條件

當系統(tǒng)執(zhí)行某一 OD后，系統(tǒng)從狀態(tài)Zl轉換到狀態(tài)Zl+1。為了保證原油的詳細調度計劃能夠完成蒸餾塔的煉油計劃，要求系統(tǒng)在時間域[0,∞)滿足文獻[9]所描述的可調度性條件。

Hou等[13]通過指派一系列的供油罐和蒸餾塔執(zhí)行ODT和ODF，實現蒸餾塔的煉油計劃。假設a時刻系統(tǒng)指派TKr和DSk，根據DSk的煉油計劃、對應的原油庫存以及以往的調度記錄，可以確定DSk需要ODF和ODT轉運的原油類型COT；為了減少蒸餾塔切換供油罐的次數，ODF和ODT取可以轉運的最大油包體積V。

式中， ζki為DSk分餾i號原油的體積，Vki為a時刻已轉運至供油罐的i號原油的體積以及初始時i號原油體積之和；totalk為a時刻可用于DSk分餾的總原油體積。ODT的起始時刻為a，結束時刻b=V/fpmax+a。ODF的起始時刻為totalk/fdsk，結束時刻為(totalk+V)/fdsk。

另外，當原油調度決策不選擇供油罐或蒸餾塔，或者沒有空閑的供油罐時，考慮停運狀態(tài)，直至系統(tǒng)的可調度性邊界再執(zhí)行下一個調度決策。

3 BPGA算法

3.1 編碼

圖2為編碼示意圖， 決策變量X=(x1,x2,···,xm)表示1個粒子或者1條染色體，對應著1個原油詳細調度計劃。編碼xn的取值范圍在區(qū)間[0,1]內，n=1,2,···,m/2。其中，x2n?1、x2n為一組，分別為供油罐和蒸餾塔的編碼?？紤]到輸油管道存在停運的情況，并且每次調度決策取盡可能大的油包V。因此，取決策變量X的編碼長度m為

圖 2 粒子/染色體編碼Figure 2 Particle/chromosome coding

其中， ???為向上取整函數；Cmin為 TK中容積最小的供油罐的容量；ξk,COT為DSk需要管道轉運的COT號原油的體積。

3.2 解碼

為了使種群具有更好的局部搜索能力，本文引入了混沌理論?；煦缡谴嬖谟诜蔷€性系統(tǒng)中的一種普遍現象，混沌運動具有有界性、遍歷性、隨機性和對初值敏感等特點。混沌映射的數學表達公式為

式中，t為進化的代數；yt為第t代的出生數；μ為控制參數；當μ∈(3.569 94,4]時，映射處于混沌狀態(tài)，T為迭代的代數。隨機生成第1個出生數y1后，剩下的(T?1)個出生數根據混沌映射依次生成，得到一個長度為T的混沌序列。

決策變量X在解碼前，需要根據混沌映射函數進行映射?；煦缬成浜瘮礖i(x)為

式中，Y(?)=y?；cei l( )為向上取整函數。X的編碼xi通過混沌映射函數，找到混沌序列上相應的值Y(?)。將X的編碼依次映射后，得到映射后的X′為

某一時刻，根據以往的調度記錄，可以得到需要管道轉運原油的蒸餾塔集合 DS和空閑的供油罐集合 TK；根據 TK中供油罐的個數和 DS中蒸餾塔的個數，將連續(xù)的[0,1]分成若干個區(qū)間長度相等的實數槽，若X′的編碼取值在某一槽中，就取槽對應編號的供油罐或蒸餾塔。例如某一時刻 TK有2個空閑的供油罐TK1和TK3，那么將[0,1]等分為2個實數槽，TK1的實數槽為[0,0.5]，TK3的實數槽為(0.5,1]，如果映射后的編碼取值為0.6，那么，指派TK3執(zhí)行ODT和ODF。類似地，根據X′的編碼指派相應的DSk。

蒸餾塔在完成煉油任務前，至少需要有一個供油罐為其供油。那么，在調度周期[τs,τe]內至少有K個供油罐處于被占用狀態(tài)。假設煉廠有H個供油罐可供使用，H≥2K+1，則在任何時刻，系統(tǒng)至多有H?K個空閑的供油罐。假設給定的系統(tǒng)初始狀態(tài)為Z0，由于管道存在停運的情況，因此Z0至多有K(H?K)+1個種子狀態(tài)。

當系統(tǒng)某一時刻根據X′的編碼指派TKr和DSk執(zhí)行ODT和ODF后，系統(tǒng)處于不可調度狀態(tài)，系統(tǒng)需要通過回溯遍歷，將不可調度狀態(tài)糾正為可調度狀態(tài)，其執(zhí)行步驟如下。

Step1初始化管道輸油速率fpmax、供油罐集合TK狀態(tài)信息、原油駐留時間 Φ、蒸餾塔煉油速率fdsk以及煉油計劃信息。

Step2判斷系統(tǒng)的初始狀態(tài)是否滿足可調度性條件。是，執(zhí)行Setp3；否，結束。

Step3根據X′的編碼指派TKr和DSk執(zhí)行調度策略。

Step4判斷調度是否可行。否，標記當前策略已使用，修改編碼，執(zhí)行Step5；是，執(zhí)行Step6。

Step5判斷當前狀態(tài)所有策略是否均已使用。是，彈出棧頂狀態(tài)，并標記當前策略已使用，修改編碼，執(zhí)行Step5；否，執(zhí)行Step3。

Step6當前狀態(tài)進棧。

Step7判斷是否完成煉油任務。是，結束；否，取下一組編碼，執(zhí)行Step3。

3.3 Pareto最優(yōu)解集

對于一個多目標問題，其解集為 Set。任意2個解X1,X2∈Set且X1X2，M為優(yōu)化的目標數量，以目標最小化為例，如果存在以下關系

則稱X1占優(yōu)于X2。對于X?∈Set，若解集 Set中不存在占優(yōu)于X?的解，則稱X?為問題的Pareto最優(yōu)解。由所有Pareto最優(yōu)解構成的集合，稱為Pareto最優(yōu)解集。Pareto最優(yōu)解集在目標空間形成的曲面，稱為Pareto前沿面(Pareto optimal front，PF)。

為了避免I-BBPSO的儲備集中解的個數過少，I-BBPSO與NSGA-II都采用快速非支配排序算法(fast non-dominated sort，FNS)構造或更新Pareto最優(yōu)儲備集[14]。

3.4 種群更新

為了提高算法的全局搜索能力和局部搜索能力，本文提出基于I-BBPSO和NSGA-II的BPGA算法，協(xié)同搜索原油詳細調度問題的Pareto最優(yōu)解。

骨干粒子群算法(bare bone particle swarm optimization，BBPSO)是由Kennedy[17]提出的，粒子的位置通過高斯采樣進行更新。為了保證算法后期的收斂性，本文引入模擬退火算法p(t)，采用以下改進的位置更新公式。

粒子群中的每個粒子記錄著其目前所找到的最優(yōu)解。該高斯采樣公式保證粒子始終在其個體引導者和全局引導者之間運動，使得粒子朝著一定的方向探索，可以保證粒子的全局搜索能力。

NSGA-II儲備集中每個個體都有相應的適應度，較優(yōu)的個體適應度較大。算法的思想是：首先隨機選擇種群中的2個個體，適應度越大的個體被選擇的概率越大，然后按照一定的概率使染色體交叉和變異，從而產生子代種群。

NSGA-II算法能夠保存種群中的較好個體，并使種群能夠朝向一定的方向進化。

3.5 移民策略

為了使BPGA具有較好的多樣性，防止算法過早陷入局部最優(yōu)，BPGA采用移民策略實現2個種群的協(xié)同計算，并且引入種群差熵，動態(tài)控制種群進行移民操作。

3.5.1 種群差熵

對于雙種群進化算法，以往的思路是通過手動調節(jié)參數來控制種群移民操作，將源種群中的最優(yōu)個體代替目標種群中的最差個體，實現多種群協(xié)同進化。本文通過引入種群差熵[18]的方法實時監(jiān)測當前種群的多樣性。種群差熵的計算公式為

式中，E(t)為種群在第t代時的Pareto熵；D(t)為第t代種群與上一代種群的Pareto差熵；O為種群中不同個體的數量;Pi表示目標值相同的個體在種群中的概率。

3.5.2 移民操作

Cmax是一個固定的常數。Pop1、Pop2分別為IBBPSO、NSGA-II的儲備集。當I-BBPSO或者NSGAII新生種群的差熵不大于0的次數達到Cmax時，取出2個種群的儲備集中一定比例的最優(yōu)個體，與另一個種群的儲備集合并，然后經過快速非支配排序FNS后，得到更新的儲備集Pop1、Pop2，實現種群間的交流。

3.6 BPGA的算法流程

Step1初始化種群。首先隨機生成N個決策變量X1-XN，解碼后得到N個個體；然后將個體1?N/2、(N/2+1)?N分別進行快速非支配排序FNS后，放入Pop1、Pop2中。2個儲備集的大小均為N/2。cnt1、cnt2均為0，分別用于統(tǒng)計D1(t)和D2(t)小于0的次數。

Step2種群更新。I-BBPSO、NSGA-II更新儲備集得到新種群R1和R2，分別與Pop1和Pop2合并后，經過快速非支配排序FNS，得到更新后的Pop1和Pop2。

Step3計算Pop1、Pop2的種群差熵D1(t)、D2(t)和cnt1、cnt2。

Step4cnt1或 cnt2是否等于Cmax。 是，進行移民操作，等于Cmax的計數器置0；否，執(zhí)行Step5。

Step5是否滿足結束條件。是，將Pop1和Pop2合并，經過快速非支配排序FNS后得到大小為N的Pareto最優(yōu)解集E；否，執(zhí)行Step2。

4 算例分析

圖3為中國南部某煉油廠的一個10 d煉油計劃，該煉油計劃共有3個蒸餾塔DS1～DS3，以及10個供油罐TK1~TK10。

圖 3 蒸餾塔10天煉油計劃Figure 3 The distiller feeding schedule for 10 days

表 1 供油罐初始狀態(tài)信息Table 1 The initial status information of charging tanks

為了說明BPGA的性能，將BPGA與多目標骨干粒 子 群(bare-bones multi-objective particle swarm optimization, BB-MOPSO)算 法[19]、多 目 標 粒 子 群(multiple objective particle swarm optimization,MOPSO) 算法[20]、NSGA-II和NSGA-III[21]等4種算法進行實驗對比。實驗均基于PlatEMO[22]平臺，采用Matlab編程實現。實驗環(huán)境為Intel(R) Core(TM)CPU 1.8 GHz 4.00 RAM，Windows 7。對于每組實驗，獨立重復執(zhí)行30次，取平均值。

在求解實際問題時，通常問題的PF前沿面是未知的。在實際的實驗中，通常通過收集所有運行的數據構造問題的PF前沿面，然后通過最大?最小標準化的方式將數據歸一化。Zitzler等[23]采用(hypervolume，HV)指標評價多目標進化算法的多樣性和收斂性的綜合性能。假設解集 Set為非支配集，本文選定參考點為(1,1,1,1)，解集 Set的HV指標值為由解集 Set中所有點與參考點在目標空間中所圍成的超立方體的體積。因此，算法求解問題得到的HV值越大，算法的綜合性能越優(yōu)。本實驗中種群大小Size為250，由圖4可以看出，當評價次數超過10 000次后，5種算法的HV值趨于穩(wěn)定，BPGA算法的綜合性能較優(yōu)。

Zhang等[24]采用C指標評價多目標進化算法得到解的質量。其公式為

圖 4 5種算法的HV指標Figure 4 HV indicator of five algorithms

式中，u為算法B迭代后儲備集中的個體；v為算法A迭代后儲備集中的個體；分子為算法B迭代結束后儲備集中被算法A中儲備集的解占優(yōu)的解的個數；分母為算法B儲備集中解的個數。因此，當C(A,B)>C(B,A)時，說明算法A儲備集的解的質量更高。 表2中，A為BPGA，B為對比算法。從表2看出，在列出的5種算法中，BPGA求得的Pareto最優(yōu)解的質量更優(yōu)。

表 2 5種算法的C指標Table 2 C indicators of five algorithms

另外，本文算法與文獻[13]實驗結果相比，不僅能得到更多的非支配解，同時部分非支配解占優(yōu)于其中的部分支配解。表3和表4列出詳細的各個非支配解的各個目標函數值。表3中方案1、6與表4中方案1、5所對應的Pareto最優(yōu)解相等；表3中方案2、3、5所對應的Pareto最優(yōu)解要優(yōu)于表4中方案2、3、4、6所對應的Pareto最優(yōu)解。從以上數據可以看出，BPGA在求解原油一次加工過程多目標優(yōu)化問題時，具有較優(yōu)的性能。值得注意的是，算法得到最終的Pareto最優(yōu)解集后，決策者需要根據實際情況選擇最適合的生產調度計劃。表4中，當設備維護致使供油罐不足時，必須犧牲其他成本以保證生產連續(xù)穩(wěn)定進行，此時應選擇方案6；當供油罐充足時，可以選擇方案1以最大限度降低其他成本。

表 3 BPGA多次迭代獲得的非支配解的各個目標函數值Table 3 The objective function values of the non-dominated solutions obtained by BPGA algorithm

表 4 文獻[13]中得到的Pareto最優(yōu)解集Table 4 Pareto optimal solution set obtained in the literature [13]

5 結束語

為了降低原油處理過程短期生產計劃問題的復雜性，本文把該問題分解為上下2層：上層求解一個目標煉油計劃以獲得最大生產率；下層求解一個詳細煉油計劃以實現上層目標煉油計劃。上層目標煉油計劃相對來說比較容易得到，但其下層詳細生產調度的獲得十分困難。本文在已有文獻將下層詳細調度問題轉化為供油罐到蒸餾塔的指派問題的基礎上，探索求解下層詳細調度過程中多個優(yōu)化目標的Pareto最優(yōu)解的新方法。為了提高多目標進化算法的多樣性和收斂性，本文提出一種基于I-BBPSO和NSGA-II協(xié)同計算的雙種群進化算法BPGA。與傳統(tǒng)的單種群進化算法相比，BPGA具有更好的收斂性和多樣性，同時能獲得質量較優(yōu)的解，且具有調節(jié)參數少、能根據種群狀態(tài)自動維護種群多樣性等優(yōu)點。本文算法在求解原油一次加工過程的詳細調度時，能夠向決策者提供多種Pareto最優(yōu)的詳細調度計劃，在實際的煉油生產中，對提高企業(yè)的利潤具有重要的現實意義。