盧卓桓, 仰燕蘭, 葉 樺

1(東南大學(xué) 自動化學(xué)院, 南京 210096)

2(復(fù)雜工程系統(tǒng)測量與控制教育部重點實驗室, 南京 210096)

郵件集散中心 (簡稱集散中心) 作為全國郵政網(wǎng)絡(luò)的重要樞紐, 主要負(fù)責(zé)對郵寄實物進(jìn)行轉(zhuǎn)發(fā)和運輸. 數(shù)據(jù)預(yù)測, 2021至2035年全社會貨運量年均增速約為2%. 穩(wěn)中有升的貨物運輸需求之下, 高價值、小批量、時效強的貨物運輸需求快速攀升. 在此背景下, 航空快件運輸憑借其時效優(yōu)勢成為集散中心的重要業(yè)務(wù)之一[1].當(dāng)前, 集散中心的航空運力主要為全貨機 (固定運力)和客機腹艙 (備選運力). 由于兩種運力方式在時間、成本、載運量等方面各異, 不同運力調(diào)度方案的運輸效率、綜合成本差異較大. 因此, 合理調(diào)度航空運力成為提升集散中心航空快件保障能力、處理效率, 實現(xiàn)降本增效的重要舉措.

目前, 國內(nèi)外關(guān)于航空運力調(diào)度問題 (air capacity scheduling problem, ACSP) 的研究主要集中在民航公司的飛機排班問題上, 包括機型指派, 機隊分配, 航班時刻優(yōu)化等. Hane等研究了一般形式下的機隊分配問題, 并建立整數(shù)規(guī)劃模型成功求解了某大型航空公司的機型分配問題[2]. Rexing等提出了帶時間窗的機隊分配模型, 并開發(fā)了速度優(yōu)先的直接算法和內(nèi)存優(yōu)先的迭代算法進(jìn)行求解[3]. Khanmirza等針對航空公司的時刻表設(shè)計和機隊分配提出了一種并行主從遺傳算法進(jìn)行求解, 該算法能夠在短時間內(nèi)解決大規(guī)模的調(diào)度問題[4]. 孫宏等在分析了國內(nèi)航空公司的航線結(jié)構(gòu)和飛機排班要求后引入了“航班節(jié)”的概念, 提出了用于描述單樞紐航線結(jié)構(gòu)下飛機排班問題的排序模型, 并設(shè)計一種分階段指派算法用于求解[5].

相比于目前的研究方向, 集散中心的航空運力調(diào)度問題 (air capacity scheduling problem of mail distribution center, ACSP-MDC) 具有一定的差異性. 一方面,集散中心航空運力調(diào)度的重點不再是客機的排班, 而是綜合調(diào)度兩種運力, 使既定目標(biāo)最優(yōu). 另一方面, 飛機排班問題是根據(jù)航班計劃和現(xiàn)有資源的情況為每個航班指定一架具體的飛機. 而ACSP-MDC不僅需要為航線選擇合適的飛機, 還需要為選中飛機分配載運量,使集散中心當(dāng)天所有的郵寄實物能順利運輸?shù)礁鱾€目的地.

ACSP-MDC作為ACSP的擴展問題, 其基本目標(biāo)是在一定航空運力資源的約束下得到合理的調(diào)度方案,使空運成本最低. 近年來, 元啟發(fā)式算法因其具有較好的適用性和魯棒性, 已經(jīng)成為解決此類問題的主流方法. 李耀華等將飛機排班計劃一體化優(yōu)化模型問題歸結(jié)為多目標(biāo)車輛路徑問題, 并提出一種自適應(yīng)單親遺傳算法進(jìn)行求解[6]. 李紅啟等構(gòu)建了航空快遞的地面集貨運力和航空運力協(xié)同調(diào)度問題的混合整數(shù)規(guī)劃模型,并設(shè)計一種兩階段啟發(fā)式算法進(jìn)行求解[7]. 本文在對比了各種元啟發(fā)式算法之后, 考慮到烏鴉搜索算法 (crow search algorithm, CSA) 具有結(jié)構(gòu)簡單, 控制參數(shù)少, 尋優(yōu)能力較強等優(yōu)點[8], 且算法的有效性和實用性已在作業(yè)車間調(diào)度、工件設(shè)計、圖像分割等不同工程領(lǐng)域中得到證明[9-11], 選擇其作為本文的基礎(chǔ)算法. 本文針對ACSP-MDC的特點, 引入懲罰函數(shù)法以消除部分約束,并在CSA的基礎(chǔ)上, 對初始化種群、烏鴉位置更新等環(huán)節(jié)進(jìn)行改進(jìn), 同時引入交叉變異機制, 提出混合交叉策略和自適應(yīng)變異策略以豐富種群的多樣性.

1 ACSP-MDC 問題模型

1.1 問題描述

集散中心經(jīng)過郵件分揀、郵件過檢、封發(fā)裝箱等前置環(huán)節(jié)之后, 需要調(diào)度航空運力, 將郵寄實物發(fā)往不同的城市. 可供集散中心調(diào)度的航空運力包括集散中心所擁有的全貨機和民航客機的腹艙. 全貨機通常在夜間飛行, 其優(yōu)勢是起飛時間靈活、運載量大, 且變動成本(單位: 元/(t·km)) 較低, 但是每啟動一架全貨機,需要額外支出維修費、起降停場費、地面服務(wù)費等固定成本, 因此, 全貨機的載運率越高, 啟動全貨機運輸會更加劃算. 客機航班的起飛時間由航空公司制定, 變動成本由航空公司或第三方貨運代理制定, 每架客機能提供的貨運量受到飛機型號和當(dāng)日旅客的行李量所限制. 因此, 集散中心在選擇備選運力時只能在航空公司制定的航班計劃中選擇合適的航班運載郵寄實物.

綜上所述, 郵件集散中心航空運力調(diào)度問題 (ACSPMDC) 可定義為: 集散中心從固定運力和備選運力中選擇合適的運力資源, 合理安排每種資源的飛行航線和載運量, 在滿足一定約束的前提下 (載重約束、航線約束等) 最小化運輸成本.

1.2 數(shù)學(xué)模型

為了便于描述ACSP-MDC的數(shù)學(xué)模型, 首先引入如下的符號定義:

K: 航線集合, K ={1,2,···,l};

I: 全貨機集合, I ={1,2,···,n};

Jk: 航線k 的客機航班集合, Jk={1,2,···,mk};

dk: 航線k的空距;

wk: 航線k的待運貨物量;

cargo_mli: 全貨機i 的最大載運量;

cargo_fci: 啟用全貨機i 的固定成本;

cargo_vci: 全貨機i 的變動成本;

civil_mlkj: 航線k 航班 j 的最大載運量;

civil_vckj: 航線k 航班 j 的變動成本;

M: 充分大的正數(shù);

xik: 全貨機i 飛 航線k 的載運量;

yik: 0-1變量, 是否使用全貨機i 飛 航線k;

ukj: 航線k航班 j 的載運量.

ACSP-MDC的數(shù)學(xué)模型為:

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

約束條件:

目標(biāo)函數(shù)式(1)由3個部分組成, 分別是全貨機的變動成本、全貨機的固定成本和客機航班的變動成本,其中, 變動成本的計算方法為: 飛機載運量×航距×飛行的單位成本; 約束(2)是航線載重約束, 表示每條航線全貨機和客機航班的載運量之和需要滿足該航線郵寄實物的待載運總量; 約束(3)-(4)分別是客機和全貨機的載重約束, 表示每種運力資源的載運量大于0且不能超過其最大載重; 約束(5)表示一架全貨機最多只能飛一條航線; 約束(6)約束了 xik和yik之間的關(guān)系, 表示只有選中全貨機 i 飛 航線k , 才會為其分配飛往航線k 的載運量, 否則全貨機i 到航線k 的載運量為0.

2 改進(jìn)烏鴉搜索算法

烏鴉搜索算法是基于烏鴉藏食和竊食的社會性行為進(jìn)行開發(fā)的, 其位置更新公式如式(7)所示:

樸素CSA的算法流程如圖1所示.

圖1 樸素CSA算法流程圖

為了使CSA能夠有效地解決ACSP-MDC, 本文根據(jù)具體場景的約束, 對CSA在求解過程中的各個步驟進(jìn)行一定程度的改進(jìn).

(1) 使用懲罰函數(shù)法消除部分約束, 拓寬烏鴉搜索的可行域, 避免在Step 5檢查烏鴉新位置可行性時, 過多烏鴉因為不滿足原有約束而被丟棄;

(2) 在Step 2初始化種群時引入基于冪函數(shù)載波的Logistic混沌映射, 提高初始種群的多樣性;

(3) 針對ACSP-MDC的特殊性, 提出了基于個體最優(yōu)追隨的位置更新策略取代Step 4的隨機跟隨策略,并嵌入正余弦算子, 提高CSA的搜索能力;

(4) 引入交叉變異機制, 根據(jù)ACSP-MDC的特點,提出兩種不同的交叉策略和一種自適應(yīng)變異策略, 維持搜索過程中種群的豐富性, 讓算法有更大概率跳出局部最優(yōu), 從而找到全局最優(yōu)解.

2.1 基于懲罰函數(shù)法的約束消除

在解決約束優(yōu)化問題時, 通過適當(dāng)?shù)匾爰s束處理技術(shù), 可以平衡目標(biāo)函數(shù)和約束條件的信息, 使優(yōu)化算法更加高效. 懲罰函數(shù)法是其中一種常見且有效的技術(shù), 其基本原理是: 將問題的約束條件轉(zhuǎn)換為一個懲罰函數(shù)并加在原目標(biāo)函數(shù)中, 構(gòu)造出一個增廣目標(biāo)函數(shù). 若當(dāng)前解不滿足約束條件, 增廣目標(biāo)函數(shù)就會加大懲罰, 淘汰可行域之外的解[12]. 通過引入懲罰函數(shù), 可以減少ACSP-MDC中的約束條件, 擴大解的可行域,避免部分烏鴉因為不滿足約束而被視為無效個體, 同時也可以省去算法中處理特定約束的步驟和時間.

本問題的約束條件為式(2)-式(6), 其中, 約束(3)和約束(4) 在搜索過程中通過邊界修復(fù)來實現(xiàn), 如式(8)和式(9)所示:

由于存在約束(5)-式(6), 使得xi=[xi0,xi1,···,xil]中只有一個維度的值大于0, 其余維度值均為0, 故可將式(9)轉(zhuǎn)化為式(10):

約束(2)和約束(5)則通過進(jìn)入懲罰項來消除, 結(jié)果如式(11)所示:

其中, σ1和 σ2為懲罰因子, 在本文中, 將其設(shè)置為足夠大的正數(shù)即可.

引入懲罰項后, 評價烏鴉位置優(yōu)劣的適應(yīng)度函數(shù)可定義為:

其中, Fc為 目標(biāo)函數(shù), Fp為懲罰項, 式(12)表明, 適應(yīng)度值越小, 烏鴉的位置越優(yōu).

2.2 初始化種群

初始化種群的好壞直接影響智能群體算法的尋優(yōu)精度和收斂速度, 樸素CSA采用隨機方法生成初始化種群, 這種生成方式有兩個弊端. 其一, 容易造成烏鴉的位置分布不均勻, 導(dǎo)致烏鴉種群難以搜索到整個空間, 易陷入早熟, 最終只能得到局部最優(yōu)解; 其二, 對于約束較為復(fù)雜的問題, 隨機生成的烏鴉容易逃逸出搜索空間的可行域. 盡管本文通過引入懲罰函數(shù), 邊界修復(fù)等方式消除了約束, 擴大了烏鴉搜索的可行域, 初始化高質(zhì)量的種群仍是必要的, 因為初始化烏鴉個體如果能盡量落在滿足原約束的可行域內(nèi), 它們將更靠近最優(yōu)解, 搜索效率會更高.

因此, 本文引入了基于冪函數(shù)載波技術(shù)的Logistic混沌映射[13]來提高烏鴉種群的多樣性, 同時借鑒文獻(xiàn)[14]對0-1變量的混沌初始化方法, 使初始的烏鴉個體更加靠近最優(yōu)位置, 為后續(xù)的搜索打下基礎(chǔ).

2.2.1 冪函數(shù)載波技術(shù)

Logistic映射方程為:

其中, zn為生成的[ 0,1]區(qū) 間的混沌值, 迭代初始時z0為隨機生成的混沌初值, 但不可取0、0.25、0.75、1;μ∈[0,4], 為控制參數(shù), 當(dāng) 取4時, 系統(tǒng)處于完全混沌狀態(tài). 雖然Logistic映射產(chǎn)生的混沌變量具有遍歷性, 但由于其軌道點分布不均勻, 導(dǎo)致落在區(qū)間兩端的點要比區(qū)間內(nèi)部的點更多. 為了糾正這一問題, 充分發(fā)揮Logistic混沌變量的遍歷性, 文獻(xiàn)[14]提出了式(14)所示的冪函數(shù)載波技術(shù):

其中, 0 ＜a＜b＜1, 0 ＜p＜1, q ＞1, zn為式(13)產(chǎn)生的混沌變量, zn′為載波后的混沌變量. 通過冪函數(shù)載波,靠近區(qū)間左邊的點會右移, 靠近區(qū)間右邊的點會左移,從而使得載波后的點分布更加均勻, 遍歷性更好, 圖2展示了這一點. 本文取a =0.15, b =0.96, p=0.46, q =13.6.

圖2 冪函數(shù)載波前后Logistic映射的遍歷性

2.2.2 基于冪函數(shù)載波的初始化方法

本文對變量 xik和ukj進(jìn)行實數(shù)編碼, 組成一個烏鴉個體. 對于備選運力ukj, 直接使用式(15)進(jìn)行初始化.

對于固定運力 xik, 由于其受到了一架全貨機只能飛行一條航線的限制, 因此, 可以先對0-1變量 yik進(jìn)行初始化, 使其滿足約束(5), 然后再根據(jù) yik的值對 xik進(jìn)行賦值, 如式(16)所示:

假設(shè)全貨機數(shù)量為n, 航線數(shù)量為l, 則對于每只烏鴉, 初始化yik(i=1,2,···,n,k=1,2,···,l)的步驟如下:

Step 1. 賦給式(13)隨機初值, 生成n 維的混沌序列seq.

Step 2. 對混沌序列應(yīng)用式(14), 將其載波成新的序列 s eq′.

Step 3. 將區(qū)間[ 0,1] l +1等分, 生成區(qū)間編號為0-l.

Step 4. 判斷序列 s eq′[i]所在的區(qū)間范圍, 并置:

由上面4個步驟實現(xiàn)的yik的初始化可以滿足約束(5), 然后應(yīng)用式(16)初始化 xik的值, 應(yīng)用式(15)初始化ukj的值, 將 xik和ukj組成一個完整的烏鴉個體, 由此生成的烏鴉比隨機生成的烏鴉有更大的概率在原約束下的可行域內(nèi)或更加靠近可行域, 且分布更加均勻. 重復(fù)迭代即可生成整個種群.

2.3 改進(jìn)的位置更新策略

2.3.1 個體最優(yōu)追隨機制

初始生成的烏鴉已經(jīng)滿足了約束(5), 此時得到的個體是較好的, 如果仍然應(yīng)用式(7)所示的樸素CSA的隨機跟隨策略更新烏鴉位置, 會導(dǎo)致大量烏鴉在探索新位置時, 遠(yuǎn)離原約束下的可行域.

假設(shè)如下情況: 集散中心共有2架全貨機, 當(dāng)天有2條航線需要飛行. 在初始化種群后, 進(jìn)入搜索階段. 僅考慮全貨機的調(diào)度, 設(shè)某兩個烏鴉個體初始化情況如式(18)所示, 烏鴉i 跟隨烏鴉 j的可能結(jié)果如圖3所示.可以看到, 跟隨后烏鴉i 已經(jīng)不滿足約束(5), 即同一架全貨機飛行兩條航線, 此時已偏離最優(yōu)解, 可認(rèn)為是一次無效跟隨.

i j圖3 烏鴉 跟隨烏鴉 的個體編碼

因此, 為了保證烏鴉每次跟蹤都能更加靠近最優(yōu)解, 本文提出了個體最優(yōu)追隨機制, 即在跟蹤階段, 每只烏鴉不再隨機跟蹤種群中的另一只烏鴉, 而是追隨上一代的個體最優(yōu)位置, 這種調(diào)整能夠保證烏鴉不斷靠近最優(yōu)解, 并提高算法的收斂速度. 改進(jìn)后的位置更新策略如式(19)所示.

2.3.2 正余弦優(yōu)化策略

考慮到固定的參數(shù)設(shè)置使得樸素CSA的搜索能力受限, 本文引入正余弦算法(sine cosine algorithm,SCA)作為局部優(yōu)化算子來提高CSA的全局探索和局部開發(fā)的能力[15,16]. 位置更新公式如式(20)所示:

其中, R1為控制參數(shù), 決定位置更新的方向, 其計算方法為.R1的大小控制著烏鴉的搜索范圍, 隨著迭代次數(shù)增加, R1逐漸減小, 烏鴉的搜索范圍隨之減小, 最終收斂在一個最優(yōu)位置, a為常數(shù), 本文取a=2; R2, R3, R4是 均勻分布的隨機數(shù), R2∈[0,2π], 控制位置更新的步長, R3∈[0,2], 控制個體最優(yōu)位置的影響程度, R4∈[0,1], 控制算法是使用正弦還是余弦操作.

正余弦算法本身具有參數(shù)自適應(yīng)的機制, 能夠很好地彌補CSA參數(shù)固定的缺點. 在迭代前期, R1較大, 算法探索全局的能力較強, 在迭代后期, R1較小, 算法的局部開發(fā)能力較強; 同時R2的隨機性也控制著正余弦函數(shù)的振幅, 使得同一個方向上的不同區(qū)段范圍都有可能被搜索到, 以此更好地平衡搜索過程中的全局探索能力和局部開發(fā)能力, 并最終收斂于最優(yōu)解或較滿意的解.

2.4 交叉變異機制

改進(jìn)后的位置更新策略能夠大大加快算法的收斂速度, 但卻不利于維持搜索過程中種群的多樣性, 導(dǎo)致算法易陷入局部最優(yōu)解. 為了讓算法有更大的機會跳出局部最優(yōu), 本文受遺傳算法啟發(fā), 引入了交叉變異機制來提高算法全局尋優(yōu)的性能.

2.4.1 混合交叉策略

在ACSP-MDC中, 全貨機和客機的約束不同, 故本文提出一種混合交叉策略, 將烏鴉個體按全貨機和客機分為兩段, 分別記為X段和U段, 各自采用不同的交叉策略.

對于X段, 為了保證交叉后的個體滿足約束(5),采用如下交叉策略: 將兩個父個體的X段按航線數(shù)量進(jìn)行切分, 切分點即為交叉點, 對于每個交叉點, 隨機生成 [ 0,1] 之 間均勻分布的隨機數(shù)r andom , 若random＜Pc(給定概率), 則進(jìn)行交叉操作.

對于U段, 如果生成的[ 0,1]區(qū) 間的隨機數(shù)random＜Pc, 則采用式(21)進(jìn)行實數(shù)交叉, 更新個體的U段.

Step 1. 在種群中隨機選取另一只烏鴉j ;

Step 2. 按照X段的交叉策略進(jìn)行交叉, 交叉后可得到烏鴉i的X段和兩個子X段;

Step 3. 按照U段的交叉策略生成新的U段;

Step 4. 將Step 2和Step 3生成的X段和U段組合起來, 生成3只烏鴉, 應(yīng)用適應(yīng)度函數(shù)并保留適應(yīng)度值最小的烏鴉.

混合交叉示意圖如圖4所示.

圖4 混合交叉示例

2.4.2 自適應(yīng)變異策略

常見的變異操作包括兩點交換變異和單點特殊變異, 受本文實際場景的約束, 一架全貨機只能飛行一條航線, 交換變異有很大概率會造成不可行解, 故本節(jié)采用單點特殊變異, 具體操作如式(22)所示:

為了保持種群的優(yōu)良性, 對較好的個體減少變異的概率, 對較差的個體增加變異的概率, 本文提出了如下的自適應(yīng)變異策略.

其中, i為烏鴉按適應(yīng)度值的大小升序排列后的序號,N 為種群大小, pmin為設(shè)定的最小變異概率, pmax為設(shè)定的最大變異概率, Pm的 取值范圍為( pmin,pmin+pmax).式(23)所示的自適應(yīng)變異概率是線性遞減的, 當(dāng)個體越優(yōu)時, Pm越小, 變異的機會小, 有利于優(yōu)良個體的保留; 當(dāng)個體較差時, Pm變大, 變異的機會隨之增大, 有利于較差個體的進(jìn)化. 采用自適應(yīng)變異策略可以在保持種群優(yōu)良性的同時, 豐富種群的多樣性.

2.5 算法總流程

改進(jìn)烏鴉搜索算法的詳細(xì)步驟如圖5所示.

圖5 改進(jìn)CSA算法流程圖

3 問題求解

以某物流公司所擁有的全貨機數(shù)據(jù)和該公司某段時間的物流情況為參考, 設(shè)計如表1-表4的測試用例,以驗證本文算法的有效性.

表1給出了航線的信息, 表2給出了可供調(diào)度的全貨機信息, 表3給出對應(yīng)航線的客機航班信息, 表4給出了7組運輸量組合, 每組都給出了不同航線需要運輸?shù)呢浳锪康那闆r.

表1 航線信息

表2 全貨機信息

表3 客機航班信息

表4 航線運輸量(t)

設(shè)置種群數(shù)為60, 迭代次數(shù)為100, 烏鴉跟蹤步長fl 為2, 感知概率 A P 為0.1, 正余弦算子的輔助常數(shù)α 為2, 交叉概率 Pc為 0.5, 最大變異概率 pmax為0.5, 最小變異概率 pmin為0.

使用本文改進(jìn)的CSA算法求解結(jié)果如表5所示,其中, 平均值和最優(yōu)值是運行50次所得到的結(jié)果, 最優(yōu)的調(diào)度方案記錄為: 航線:飛機編號(載重). 從表中可以看到, 每組算例的調(diào)度方案都能滿足各自航線的載重需求, 同時, 每組結(jié)果的平均值和最優(yōu)值相差很小,最大的一組W7僅相差了318.96, 表明本文改進(jìn)的CSA算法具有較好的魯棒性.

表5 改進(jìn)CSA算法求解結(jié)果

圖6展示了樸素CSA, SCA-CSA和改進(jìn)的CSA三種算法在求解W1-W7時的表現(xiàn). 從圖中可以看出改進(jìn)的CSA相比于樸素CSA和SCA-CSA具有更好的求解精度, 在面對不同問題時改進(jìn)的CSA也表現(xiàn)出更好的魯棒性.

圖6 不同算法在不同貨物量組合下的最小成本對比

針對W7這組測試用例, 圖7展示了3種算法的收斂曲線, 可以看出改進(jìn)CSA算法的求解精度和收斂速度都優(yōu)于其他兩種算法.

圖7 W7測試用例下3種算法的收斂曲線

進(jìn)一步擴展問題的規(guī)模, 表6所示的8個問題分別是對航線, 全貨機和客機航班進(jìn)行一定程度的擴充,并為每條航線隨機分配該航線的貨物量. 對每個問題都進(jìn)行50次的測試, 測試結(jié)果如表7所示. 從表中可以看出, 本文改進(jìn)的CSA在求解ACSP-MDC問題時的表現(xiàn)優(yōu)于樸素的CSA和SCA-CSA.

表6 擴充仿真實例

表7 3種算法結(jié)果(適應(yīng)度值)對比

4 結(jié)束語

本文針對郵件集散中心航空運力調(diào)度問題, 在航空運力資源充足的情況下, 建立了以最小化運輸成本為目標(biāo)的優(yōu)化模型. 針對ACSP-MDC的特點, 本文在樸素CSA的基礎(chǔ)上, 做了一系列的改進(jìn), 包括引入懲罰函數(shù)消除部分約束, 擴大烏鴉搜索的可行域, 根據(jù)問題的實際場景提出新的位置更新策略, 并引入交叉變異機制, 豐富種群的多樣性, 提高全局尋優(yōu)的性能. 改進(jìn)后的CSA相較于樸素CSA和SCA-CSA的收斂速度更快, 求解結(jié)果更優(yōu), 魯棒性更好, 且能有效地解決郵件集散中心的航空運力調(diào)度問題, 為集散中心空運業(yè)務(wù)降本增效提供合理的解決方案.