邵 荃，梁斌斌，朱 燕，張海蛟，張金石

(1.南京航空航天大學(xué) 民航學(xué)院，江蘇 南京210016;2.中國民用航空新疆管理局，新疆 烏魯木齊820016)

民航應(yīng)急物資調(diào)度是指利用民航運輸機(jī)進(jìn)行應(yīng)急物資調(diào)度的過程，因其速度快、航程遠(yuǎn)的特點在應(yīng)急救援體系中起著重要作用。一般而言，民航應(yīng)急物資調(diào)度需要多機(jī)場、多物資協(xié)同參與，利用有限運力將有限的物資運抵災(zāi)區(qū)。國外學(xué)者針對地面應(yīng)急物資調(diào)度問題作了大量研究，普遍采用優(yōu)化算法進(jìn)行求解［1-4］;國內(nèi)學(xué)者針對多目標(biāo)調(diào)度問題進(jìn)行了研究［5-6］，設(shè)計了不同的路徑優(yōu)化數(shù)學(xué)模型［7-8］。OZDAMAR［9］基于航路管理程序(RMP)對直升機(jī)應(yīng)急物資調(diào)度問題進(jìn)行求解。祁明亮等［10］研究了車輛與直升機(jī)聯(lián)合調(diào)度問題?？梢钥闯觯F(xiàn)有研究大多是關(guān)于車輛或者直升機(jī)應(yīng)急調(diào)度，對民航應(yīng)急物資調(diào)度的研究十分有限。不同于前兩種調(diào)度方式，在民航應(yīng)急物資調(diào)度中，機(jī)場無法保障的機(jī)型無法飛往該機(jī)場，短航程機(jī)型無法進(jìn)行遠(yuǎn)距離運輸，飛行地面保障時間無法忽略，這些特點決定了有必要對民航應(yīng)急物資調(diào)度進(jìn)行單獨研究。遺傳算法是解決物資調(diào)度問題的常用優(yōu)化算法之一，但其穩(wěn)定性較差且容易早熟。由其衍生出的完整元胞遺傳算法有較好的穩(wěn)定性和抑制早熟的能力［11-12］。當(dāng)前，少有研究采用元胞遺傳算法求解應(yīng)急調(diào)度問題，常用的元胞遺傳算法極少考慮個體適應(yīng)度對元胞狀態(tài)演化的影響，在尋優(yōu)能力和收斂性方面均有較大缺陷，有待進(jìn)一步改進(jìn)。

根據(jù)民航應(yīng)急物資調(diào)度的特點，建立不同于地面車輛和直升機(jī)調(diào)度的多機(jī)場多物資聯(lián)合調(diào)度模型，并基于適者生存的生態(tài)法則提出一種與個體適應(yīng)度相關(guān)的元胞狀態(tài)演化規(guī)則，可有效提高算法的尋優(yōu)能力和收斂性能，便于更好地求解模型。實例表明，該模型合理可行，算法較標(biāo)準(zhǔn)遺傳算法和傳統(tǒng)元胞遺傳算法均具有優(yōu)勢。

1 民航應(yīng)急救災(zāi)物資調(diào)度原理

民航應(yīng)急救災(zāi)物資調(diào)度問題可描述為：存在若干個救災(zāi)機(jī)場(運出物資)和若干個受災(zāi)機(jī)場(接收物資)。當(dāng)災(zāi)情發(fā)生后，飛機(jī)在救災(zāi)機(jī)場裝載物資，運輸物資前往具有保障能力的受災(zāi)機(jī)場，卸載物資之后返回任一具有保障能力的救災(zāi)機(jī)場，如此往返直至救災(zāi)任務(wù)結(jié)束。為方便理論建模，假設(shè)物資為統(tǒng)一規(guī)格的包裝單元，飛機(jī)總是裝載整數(shù)個單元的物資;運輸過程中飛機(jī)總是以巡航速度飛行。該問題的數(shù)學(xué)表達(dá)式如下：

救災(zāi)機(jī)場集合D={d1，d2，…，dn}，表示共有n個救災(zāi)機(jī)場，di表示第i個救災(zāi)機(jī)場;

受災(zāi)機(jī)場集合E={e1，e2，…，em}，表示共有m個受災(zāi)機(jī)場，ej表示第j個受災(zāi)機(jī)場;

應(yīng)急物資種類集合C={c1，c2，…，cl}，表示共有l(wèi)種應(yīng)急物資，cI表示第I種應(yīng)急物資;

救災(zāi)機(jī)場物資儲備量矩陣，dmiI表示應(yīng)急物資I在救災(zāi)機(jī)場i的儲備數(shù)量：

受災(zāi)機(jī)場物資需求量矩陣，emjI表示應(yīng)急物資I在受災(zāi)機(jī)場j的需求數(shù)量：

受災(zāi)機(jī)場物資需求緊要度矩陣，ujI表示物資I在受災(zāi)機(jī)場j的緊要程度，其中…，m。

可用飛機(jī)集合A={A1，A2，…，Ap}，表示共有p架飛機(jī)執(zhí)行調(diào)度任務(wù)，Ah表示第h架飛機(jī);

飛機(jī)載運量集合L={L1，L2，…，Lp}，Lh表示第h架飛機(jī)的載運量;

飛機(jī)航程參數(shù)集合R={R1，R2，…，Rp}，Rh表示第h架飛機(jī)的航程;

飛機(jī)巡航速度集合v={v1，v2，…，vp}，vh表示飛機(jī)h的巡航速度;

飛機(jī)地面保障時間t={t1，t2，…，tp}，th表示飛機(jī)h所需的地面保障時間;

機(jī)場飛機(jī)可用度矩陣：

(Oij=1 或0)，其中Oij=1 表示機(jī)場j可以保障飛機(jī)i，Oij=0 表示機(jī)場j不能保障飛機(jī)i。表示救災(zāi)機(jī)場i處飛機(jī)h的可用度，Oehj表示受災(zāi)機(jī)場j處飛機(jī)h的可用度;

航線距離集合r={rij|i=1，2，…，n;j=1，2，…，m}，rij表示救災(zāi)機(jī)場i到受災(zāi)機(jī)場j的航線距離;

飛機(jī)h的一次飛行任務(wù)τs包括出發(fā)的救災(zāi)機(jī)場、運載的物資種類及數(shù)量和抵達(dá)的受災(zāi)機(jī)場。飛 行 任 務(wù)表示飛機(jī)h在飛行任務(wù)τs中從救災(zāi)機(jī)場i向受災(zāi)機(jī)場j運載第I種物資的數(shù)量;

應(yīng)急調(diào)度方案φ 由所有飛機(jī)的調(diào)度方案組成，即φ= {φ1，φ2，…，φp}，其中φh表示飛機(jī)h的調(diào)度任務(wù)序列，即φh= {τ1，τ2，…，τN}，表示任務(wù)序列φh由N次飛行任務(wù)組成;

應(yīng)急調(diào)度方案φ 的完成時間為Tφ，飛機(jī)h執(zhí)行其調(diào)度任務(wù)序列φh的時間為：

如果任務(wù)序列中存在無法保障飛機(jī)h的機(jī)場則對應(yīng)的任務(wù)完成時間為無窮大(Th=∞)。由此可得到整個調(diào)度方案的完成時間Tφ=max{Th}(h=1，2，…，p)。

救災(zāi)效果滿意度指數(shù)γ 表示運輸方案能夠滿足受災(zāi)需求的程度。統(tǒng)計受災(zāi)機(jī)場物資需求量矩陣EM中各個物資需求元素被滿足的程度。元素emjI被滿足的程度可以通過物資運抵總量與機(jī)場物資需求量的比值進(jìn)行衡量，即為其中i，j∈τs，I=0，1，…，l。

則有：

γ 反映不同物資需求緊要度約束下的物資需求滿足程度。γ =0 表示所有機(jī)場的物資需求均不被滿足，救災(zāi)效果滿意度最低;γ =1 表示所有機(jī)場的物資需求均被滿足，救災(zāi)效果滿意度最高。

民航應(yīng)急物資調(diào)度受到時間、物資儲備數(shù)量和飛機(jī)性能的約束。一般認(rèn)為，災(zāi)害發(fā)生后3 天內(nèi)為黃金救援時間，因此要求調(diào)度方案的完成時間Tφ不能大于72 h。救災(zāi)點的物資儲備量是有限的，從每個救災(zāi)點運出的應(yīng)急物資總量不得超過救災(zāi)點的儲備總量。同時，受制于飛機(jī)的運載能力，飛機(jī)每次運載的物資總量不得大于飛機(jī)的額定載運量。考慮飛機(jī)的有限航程，飛機(jī)h的調(diào)度任務(wù)序列φh中單次飛行的最遠(yuǎn)航線距離定為max{rφh}，則應(yīng)有max{rφh}≤Rh。

應(yīng)急調(diào)度有如下兩個決策目標(biāo)：

(1)要求應(yīng)急調(diào)度方案的執(zhí)行總時間最短，即F1=min { max{Th}}(h=1，2，…，p);

(2)要求應(yīng)急調(diào)度方案的救災(zāi)效果最好，即F2=max{ γ} 。

為方便建立單目標(biāo)數(shù)學(xué)模型，將F2變形為F2=min {1 -γ }。由此建立如下數(shù)學(xué)模型：

式(4)中，ω1，ω2分別為決策因子z1和z2的權(quán)重系數(shù)，設(shè)計如下權(quán)重系數(shù)：

為確保決策因子受權(quán)重系數(shù)的約束，對決策因子進(jìn)行歸一化處理，得到單目標(biāo)函數(shù)：

其中，T*=72，表示運輸時間上限。

2 民航應(yīng)急調(diào)度優(yōu)化算法

2.1 算法流程

常用的元胞狀態(tài)演化只與鄰域內(nèi)元胞的生死狀態(tài)有關(guān)，而與個體適應(yīng)度無關(guān)，這極有可能導(dǎo)致優(yōu)秀個體的迅速消亡，從而降低算法的尋優(yōu)能力?；谶m者生存的生態(tài)法則，提出一種基于個體適應(yīng)度的元胞狀態(tài)演化規(guī)則，算法流程如圖1 所示。

圖1 算法流程圖

2.2 遺傳編碼

將所有飛機(jī)的調(diào)度任務(wù)序列組成一條染色體。每架飛機(jī)的飛行任務(wù)序列組成一個基因片段，每個基因片段中的單次飛行任務(wù)為一個基因單元，每個基因單元中的救災(zāi)機(jī)場、物資運輸量和受災(zāi)機(jī)場統(tǒng)稱為基因位點。例如，現(xiàn)有救災(zāi)機(jī)場d1，d2，d3，受災(zāi)機(jī)場e1，e2和應(yīng)急物資c1，c2，c3。則飛機(jī)h對應(yīng)的染色體基因片段編碼格式可以為表示飛機(jī)h從救災(zāi)機(jī)場d1出發(fā)前往受災(zāi)機(jī)場e2，運載的3 類應(yīng)急物資數(shù)量分別為之后返回救災(zāi)機(jī)場d2，再運載物資前往受災(zāi)機(jī)場e1，運載的物資數(shù)量分別為直至任務(wù)結(jié)束。由此，整個應(yīng)急調(diào)度方案可用如下編碼表示：

其中，“||”將每架飛機(jī)的基因片段分隔開。

2.3 元胞空間初始化

初始元胞群體應(yīng)隨機(jī)選取，元胞空間不能太大，否則影響算法收斂速度;但也不能太小，否則可能引起早熟現(xiàn)象。為提高算法效率，應(yīng)在種群生成時對染色體進(jìn)行可行性驗證，檢驗染色體是否滿足約束條件。種群初始化后，將個體依次放置在元胞網(wǎng)格內(nèi)，并隨機(jī)確定元胞的生死狀態(tài)。

2.4 元胞遺傳操作

(1)選擇。依次選擇生存元胞作為中心元胞，并在鄰域內(nèi)選取適應(yīng)度最大的生存元胞進(jìn)行交叉和變異操作，得到新元胞的適應(yīng)度，若其優(yōu)于中心元胞適應(yīng)度，則替代中心元胞，否則被舍棄。

(2)交叉。在雙親確定后，按照交叉概率pcross給每個基因片段分別選取一個交叉隨機(jī)位，對換父代染色體和母代染色體中的每個基因片段上對應(yīng)隨機(jī)位點之后的基因單元，形成兩個子代。例如，父代染色體為：

母代染色體為：

設(shè)定第一個基因片段的交叉隨機(jī)位為2，第二個基因片段的交叉隨機(jī)位為1，則交叉生成的兩條子代染色體分別為：

子代染色體1：

子代染色體2：

對生成的子代染色體進(jìn)行驗證，如果滿足驗證要求，則進(jìn)行后續(xù)操作;否則，將其舍棄，重新執(zhí)行交叉操作。

(3)變異。生成子代染色體后，按照變異概率pmutate給每個基因片段分別選取一個變異隨機(jī)位，將每個基因片段上變異隨機(jī)位對應(yīng)的基因位點數(shù)值突變?yōu)榱硪唤M數(shù)值，生成新的子代染色體。如果子代滿足驗證要求，則加入下一代種群;否則，將其舍棄，重新執(zhí)行變異操作。

2.5 元胞狀態(tài)演化

針對完整元胞遺傳算法忽視個體適應(yīng)度，降低了算法尋優(yōu)能力的缺陷，筆者基于適者生存法則對元胞狀態(tài)演化規(guī)則進(jìn)行了改進(jìn)。自然界中，個體存活取決于所處的生態(tài)環(huán)境和自身生存能力。生態(tài)環(huán)境越好，個體越容易存活;生存能力越強，個體也越容易存活。在元胞遺傳算法中，生態(tài)環(huán)境可以通過生存元胞和死亡元胞的數(shù)量關(guān)系進(jìn)行衡量。生存元胞越多，表明生態(tài)環(huán)境越好;反之，表明生態(tài)環(huán)境越差。同理，適應(yīng)度越優(yōu)，生存能力越強;適應(yīng)度越差，生存能力越弱。將鄰域空間(包含中心元胞)視為中心元胞的生存環(huán)境，將環(huán)境中的元胞分為生存(S=1 )和死亡(S=0 )兩個群體，并設(shè)定如下狀態(tài)演化規(guī)則：

當(dāng)元胞(i，j)在t時刻狀態(tài)為生存，即1 時，其在t+1 時刻的狀態(tài)及概率為：

當(dāng)元胞(i，j)在t時刻狀態(tài)為死亡，即0 時，其在t+1 時刻的狀態(tài)及概率為：

式中：S=1 或S=0 表示元胞處于生存或死亡狀態(tài);N為元胞個體生存環(huán)境中的元胞數(shù)量;Nt(S=1 )與Nt(S=0 )分別表示個體生存環(huán)境中狀態(tài)為生存與死亡的元胞數(shù)量;Rank1(i，j) 為元胞(i，j) 的適應(yīng)度在生存群體中的排名;Rank0(i，j) 為其在死亡群體中的排名，適應(yīng)度從差到優(yōu)進(jìn)行排列。

2.6 控制參數(shù)及終止條件

元胞遺傳算法的控制參數(shù)主要有鄰域類型、種群規(guī)模ps、交叉概率pcross、變異概率pmutate和終止代數(shù)等。選定的鄰域類型為Moore 類型，其模型包括中心元胞及周圍8 個鄰居元胞;群體規(guī)模ps=100，元胞空間大小為10×10;交叉概率pcross=0.9，變異概率pmutate=0.05。算法終止條件設(shè)為滿足兩個中的任意一個即可：①若迭代次數(shù)達(dá)到設(shè)定的最大值4 000;②若某代染色體的平均適應(yīng)度達(dá)到這一代最佳染色體適應(yīng)度的0.96。

3 實驗結(jié)果及分析

3.1 初始條件設(shè)定

現(xiàn)有P1，P2，P3這3 種機(jī)型共7 架飛機(jī)從救災(zāi)機(jī)場d1，d2，d3運輸藥品、食品、生活用品3 種應(yīng)急物資前往受災(zāi)機(jī)場e1，e2。已知機(jī)場d1無法保障P2機(jī)型，機(jī)場e2無法保障P1機(jī)型;救災(zāi)機(jī)場與受災(zāi)機(jī)場之間的航線距離如表1 所示。

表1 救災(zāi)機(jī)場與受災(zāi)機(jī)場航線距離 km

救災(zāi)機(jī)場的物資儲備和受災(zāi)機(jī)場的物資需求如表2 所示;機(jī)型參數(shù)如表3 所示。

表2 物資儲備和物資需求表 單元

表3 機(jī)型參數(shù)表

3.2 調(diào)度方案分析

結(jié)合應(yīng)急調(diào)度原理和數(shù)學(xué)模型的決策目標(biāo)，將調(diào)度方案執(zhí)行總時間和救災(zāi)效果滿意度作為調(diào)度方案可行性評價參數(shù)。執(zhí)行改進(jìn)元胞遺傳算法，迭代4 000 次得到最滿意調(diào)度方案，如表4 所示。

表4 改進(jìn)元胞遺傳算法求解的最滿意調(diào)度方案

該方案要求各架飛機(jī)執(zhí)行任務(wù)的次數(shù)分別為9，9，8，11，11，6，7，運輸效率較高的大機(jī)型P1和P2的任務(wù)次數(shù)較多，運輸效率最低的小機(jī)型P3的任務(wù)次數(shù)最少。一般而言，機(jī)型越大，載重越大，速度越快，運輸效率越高。因此，大機(jī)型被調(diào)用的優(yōu)先級別高于小機(jī)型。受機(jī)場保障能力的限制，3 架P1型飛機(jī)#1，#2 和#3 只能前往受災(zāi)機(jī)場e1;兩架P2型飛機(jī)#4，#5 只能在救災(zāi)機(jī)場d2和d3裝載物資。受機(jī)型航程的限制，兩架P3型飛機(jī)#6，#7 只能從救災(zāi)機(jī)場d2裝載物資前往受災(zāi)機(jī)場e1，e2。執(zhí)行該方案后各機(jī)場的物資量變動情況如表5 所示。

對比表2 和表5 可知，該方案的救災(zāi)效果滿意度指數(shù)為1，各個受災(zāi)機(jī)場的物資需求均被滿足，救災(zāi)效果達(dá)到最佳。最滿意方案的執(zhí)行總時間為52.551 h，明顯低于72 h。實驗結(jié)果表明，該數(shù)學(xué)模型合理可行。

表5 最滿意方案下各機(jī)場的物資變動表 單元

3.3 不同算法下的方案求解及對比分析

將標(biāo)準(zhǔn)遺傳算法(SGA)、完整元胞遺傳算法(CEGA)和改進(jìn)元胞遺傳算法(MCGA)分別用于求解該算例。每個算法執(zhí)行100 次，每次迭代4 000次并生成一組調(diào)度方案，每個算法均生成100 組調(diào)度方案。引入最滿意函數(shù)值、整體平均值和函數(shù)值均方差作為執(zhí)行一次算法的性能指標(biāo)。對執(zhí)行100 次算法的性能指標(biāo)結(jié)果取平均值，得到如表6 所示的性能指標(biāo)平均值對比表。

表6 3 種算法的性能指標(biāo)平均值對比表

對100 組調(diào)度方案的最滿意函數(shù)值取平均值，得到3 種算法的最滿意函數(shù)值的平均值分別為0.787，0.746 和0.675，MCGA 算法搜尋出的最滿意函數(shù)值明顯優(yōu)于SGA 和CEGA 算法。對100組調(diào)度方案中的所有迭代結(jié)果取平均值，得到3種算法的整體平均值分別為0. 815，0. 772 和0.698，說明MCGA 算法的尋優(yōu)效率明顯增強。前兩組數(shù)據(jù)表明CEGA 和MCGA 在尋優(yōu)能力和整體求解精度上較SGA 均有較大提高，而MCGA比CEGA 更勝一籌。對100 組調(diào)度方案得到的函數(shù)值均方差取平均值，得到均方差平均值分別為0.033，0.042 和0.037，說明MCGA 在進(jìn)化收斂性上也比CEGA 略勝一籌。

4 結(jié)論

結(jié)合民航應(yīng)急物資調(diào)度的特點，針對多機(jī)場民航應(yīng)急物資聯(lián)合調(diào)度問題，建立了以飛機(jī)性能和物資需求為約束，以運輸時間最短、救災(zāi)效果最好為決策目標(biāo)的數(shù)學(xué)模型。求解過程中提出一種改進(jìn)元胞遺傳算法，針對傳統(tǒng)元胞遺傳算法尋優(yōu)能力不足的缺點，基于生態(tài)法則對元胞狀態(tài)演化規(guī)則進(jìn)行改進(jìn);最后設(shè)置算例進(jìn)行求解。算例的最滿意調(diào)度方案用時較短且取得了最佳的救災(zāi)效果，結(jié)果表明運輸效率較高的大機(jī)型應(yīng)優(yōu)先調(diào)用，算例結(jié)果與實際情況相符，證明數(shù)學(xué)模型合理可行。對比標(biāo)準(zhǔn)遺傳算法和完整元胞遺傳算法，改進(jìn)元胞遺傳算法在尋優(yōu)能力和收斂性兩方面均有明顯優(yōu)勢。

［1］ BERKOUNE D，RENAUD J，REKIK M，et al. Transportation in disaster response operations［J］. Socio -Economic Planning Sciences，2012，46(1)：23 -32.

［2］ LIU N，YE Y.Humanitarian logistics planning for natural disaster response with Bayesian information updates［J］. Journal of Industrial and Management Optimization，2013，10(3)：665 -689.

［3］ BOZORGI -AMIRI A，JABALAMELI M S，AL -E -HASHEM S M J. A multi -objective robust stochastic programming model for disaster relief logistics under uncertainty［J］.OR Spectrum，2013，35(4)：905-933.

［4］ SHEU J B.Dynamic relief-demand management for emergency logistics operations under large-scale disasters［J］. Transportation Research Part E：Logistics and Transportation Review，2010，46(1)：1 -17.

［5］ 文仁強，鐘少波，袁宏永，等.應(yīng)急資源多目標(biāo)優(yōu)化調(diào)度模型與多蟻群優(yōu)化算法研究［J］. 計算機(jī)研究與發(fā)展，2013，50(7)：1464 -1472.

［6］ 王劍，王紅衛(wèi).多目標(biāo)資源受限的運輸調(diào)度問題研究［J］.武漢理工大學(xué)學(xué)報，2008，30(5)：155- 158.

［7］ 梁勤歐.基于改進(jìn)免疫算法的有能力約束車輛路徑問題［J］. 武漢理工大學(xué)學(xué)報(信息與管理工程版)，2011，33(5)：763 -767.

［8］ 高鴻鶴，唐辰. 基于配送時間最短的應(yīng)急物流路徑規(guī)劃［J］.物流工程與管理，2014，36(2)：75 -78.

［9］ OZDAMAR L. Planning helicopter logistics in disaster relief［J］.OR Spectrum，2011，33(3)：655 -672.

［10］ 祁明亮，秦凱杰，趙琰.雪災(zāi)救援物資車輛-直升機(jī)聯(lián)合運送的調(diào)度問題研究［J］.中國管理科學(xué)，2014，22(3)：59 -67.

［11］ 黎明，王瑩，陳昊，等.基于捕食機(jī)制的元胞遺傳算法［J］.應(yīng)用科學(xué)學(xué)報，2012，30(6)：669 -676.

［12］ KORNYAK V V. Symmetric cellular automata［J］.Programming and Computer Software，2007，33(2)：87 -93.