求解廣義分配問題的拉格朗日蝙蝠算法

萬曉瓊，張惠珍，趙玉蘋

（1. 上海理工大學 管理學院，上海 200093；2. 國網上海市電力公司 物資公司，上海 200235）

廣義分配問題（generalized assignment problem，GAP）[1]是將給定的工作分配給有限定資源量的代理器來完成，以消耗最小的成本或取得最大利潤為目標，是經典的組合優化問題。GAP在實際生活中具有廣泛的應用，許多領域的容量約束問題都可以被抽象為GAP算例進行求解，如機器調度問題、有容量設施選址問題、供應鏈問題及車輛路徑問題等。

GAP在理論和應用價值上的意義引起了國內外學者的極大關注，多種算法被用來求解此問題，這些算法大致可以歸為3類：近似算法[2-3]、精確算法[4-6]和智能優化算法。由于GAP的復雜性，近似算法和精確算法對于求解大規模問題具有局限性，這促使了智能優化算法的廣泛應用，如Bai等[7]采用粒子群算法（PSO）求解兩階段模糊GAP；Tapkan等[8]在蜂群算法中引入多種模糊排序方式用于解決模糊多目標GAP；Liu等[9]提出一種可擴展并行遺傳算法（GA）求解GAP問題；Laguna等[10]在禁忌搜索算法（TS）中利用斥鏈定義鄰域，并用其求解多級GAP。它們在解決組合優化問題時取得了一定的成果，同時也顯現出單一智能優化算法的一系列缺點，如容易陷入局部最優解、穩定性差等。為了克服這些缺點，以及獲得更精確的解，運用混合智能優化算法解決GAP問題得到長足的發展，如吳良杰等[11]將禁忌搜索算法作為蟻群算法的局部搜索策略，提出了混合算法，有效地解決GAP；Chen等[12]結合小生境遺傳算法與蟻群算法，并建立螞蟻搜索的禁忌規則，提出了小生境遺傳螞蟻算法來求解GAP；Osman[13]將禁忌搜索算法與模擬退火算法相結合，用于求解GAP。

本文運用蝙蝠算法（BA）[14]融合拉格朗日松弛算法（LR）[15-16]來解決基本GAP。該算法利用蝙蝠算法產生GAP最優解的上界，采用拉格朗日松弛算法求出GAP最優解的下界，兩者結合，逐漸逼近GAP的最優解。

1 廣義分配問題

GAP可以描述為：將n個相互獨立的工作分配給m個有限定資源量的代理器，一個工作只能由一個代理器服務，一個代理器可以服務多個工作，且代理器完成工作所需的總資源量不能超過自身的限制。GAP的數學模型可表述為

2 蝙蝠算法

蝙蝠算法最早由劍橋大學的Yang于2010年提出，是一種基于微型蝙蝠群體回聲定位行為的搜索全局近似解的智能優化算法。蝙蝠個體是蝙蝠算法的基本單元，它們通過口腔或者鼻腔將從喉部產生的超聲波發射出去，然后依據發出超聲波和接收到回音之間的時差、到達兩耳之間的不同強度、回聲定位波形的變化來判斷獵物或障礙物的距離、方位以及性質，從而制定捕捉獵物和躲避障礙物的計劃，同時也使得蝙蝠即使在黑暗中也能找到位于裂縫中的棲息地。蝙蝠算法最早被用來解決連續問題，近幾年有學者利用它來解決離散問題，例如，旅行商問題[17]、車輛路徑問題[18]及背包問題[19]等。

本文嘗試運用蝙蝠算法來解決基本GAP問題。在求解過程中，假設蝙蝠群體隨機散布于d維搜索空間下，一個蝙蝠對應GAP問題的一個解，GAP問題的目標函數決定蝙蝠的適應度值，蝙蝠群體的其他個體通過調整超聲波的響度、脈沖發射率來追隨當前群體中的最優蝙蝠在解空間中進行搜索，從而使蝙蝠群體在求解空間中的運動由無序演變為有序，直至達到滿意解。

2.1 求解GAP問題的蝙蝠算法

由于GAP問題的特點，基礎的蝙蝠算法已不適用。本文在不改變蝙蝠算法原有概念和結構的基礎上，引入單點基因換位算子建立全新的速度、位置以及局部搜索更新規則，稱之為離散蝙蝠算法（discrete bat algorithm，DBA）。現從蝙蝠位置初始編碼、蝙蝠位置和速度更新、局部更新以及音量和脈沖的更新這4個方面來介紹離散蝙蝠算法。

2.1.1 蝙蝠位置初始編碼蝙蝠算法中蝙蝠的位置編碼代表GAP的可行解。設蝙蝠種群大小為s，第i只蝙蝠的速度為vi，初始位置為xi，迭代次數N。現介紹第i只蝙蝠個體初始位置編碼方式。

c. 從θi1開始進行分配，分配方式：將給定的代理器服務該工作時消耗的資源量按升序排列，從中隨機選取前兩個代理器中的任一個為其服務，若選中的代理器為w，則記，依次進行，直至全部工作分配完畢。

d. 若此時有代理器消耗的總資源量超出自身資源限制，則通過以下方式進行調整：搜索超出資源限制的代理器服務的工作，并計算代理器完成該工作時消耗單位資源花費的成本，將最大P值對應的工作重新分配給剩余資源量足夠且花費成本最小的代理器。

2.1.2 蝙蝠速度和位置更新

式中：x*表示當前最優解，它是通過比較(當前整)個蝙蝠種群對應的適應度值而得到；指和x*之間的漢明距離，具體為第i只蝙蝠在第t-1時刻位置編碼較最優蝙蝠位置編碼對應位不同的個數。例如，可以看出，在，x*中工作1，4，5均被分配給1，2，1代理器，而工作2，3分配給不同的代理器，則此時=2。

蝙蝠的位置更新采用單點基因換位算子，即在蝙蝠原位置編碼的基礎上隨機選取2個編碼位置進行交換，換位后的位置編碼對應的適應度值變優且符合約束條件時更新編碼，直至共進行vit次單點基因換位，位置更新結束。以m=3，n=7為例，隨機選取的2個編碼位置為2和5，則單點基因換位的方式如圖1所示。

圖1 單點基因換位方式Fig.1 Single point gene translocation

2.1.3 局部更新

類似于現存的其他智能優化算法，蝙蝠算法也可以進行局部搜索行為，且相對于其他算法，蝙蝠算法可以通過比較脈沖發生率ri和隨機數的大小，實現對全局搜索和局部搜索之間轉換的動態控制。

當隨機函數rand＞ri時，在當前最優解集中選擇一個最優蝙蝠，進行隨機游走，產生局部新解。游走方式同樣采用單點基因換位，具體是在當前最優解集中選擇一個最優蝙蝠，在該蝙蝠原有的位置編碼的基礎上進行一次單點基因換位，更新后的位置編碼對應的適應度值變優且在符合約束條件時更新編碼。

2.1.4 蝙蝠音量和脈沖發生率的更新

在蝙蝠接近獵物的過程中，聲波音量逐漸降低，同時脈沖發生率逐漸提高。音量Ait和脈沖發生率rit的更新公式為

式中：α為音量衰減系數；γ為脈沖發生率增加系數。對任意 0＜α＜1，γ＞0，可以看出逐漸接近最大脈沖發生率。

2.2 求解GAP問題的蝙蝠算法設計

通過上文方式進行基礎蝙蝠算法的改造，設計出的全新求解GAP問題的蝙蝠算法的具體步驟如下：

步驟 1 初始化蝙蝠種群。蝙蝠種群個數為s，第i只蝙蝠初始位置為xi，速度為vi，音量為Ai0，脈沖發生率為ri0，迭代次數N。

步驟 2 運用式（5）以及單基因換位算子進行蝙蝠速度和位置的更新。

步驟 3 若rand＞ri，進行局部搜索。

步驟 4 若rand＜Ai且，接受這個新解。

步驟 5 根據式（6）和式（7）進行蝙蝠脈沖發生率和音量更新。

步驟 6 排列蝙蝠并找到當前最優蝙蝠x*。

N′=N′+1

步驟 7 ，若 ，則結束搜索并輸出結果；否則，轉步驟2。 為此刻的迭代次數N′＞N N′

3 拉格朗日蝙蝠算法

3.1 求解GAP的拉格朗日松弛算法

拉格朗日松弛算法（Lagrangian relaxation，LR）是一種分解協調算法，由Herd和Karp于1970年提出，最早用于解決旅行商問題。隨后被用于各大領域，現已成為解決組合優化問題時產生解的下界的一種不可或缺的方法。拉格朗日松弛算法的基本原理是通過引入拉格朗日乘子λj，將原問題中的某些約束作為懲罰項吸收到目標函數中，使原問題成為較易解決的拉格朗日松弛問題。通過拉格朗日乘子松弛約束式（3），松弛后的GAP可表示為

s.t. 式（2）和式（4）

再通過標準次梯度優化算法修正拉格朗日乘子，求得拉格朗日對偶問題（Lagrangian dual，LD）的解。其中，拉格朗日對偶問題表示為

s.t. 式（2）和式（4）

利用標準次梯度優化算法，第t時刻對拉格朗日乘子的修正為

式（11）和（12）為拉格朗日乘子步長的更新；

式（13）表示拉格朗日乘子的更新。

通過以上方式求解拉格朗日對偶問題，解的目標函數值為原問題最優解的下界。但由于約束條件被松弛，擴大了GAP問題可行解的范圍，導致松弛后的解不一定是原問題的可行解。

為了保證解的可行性，本文對拉格朗日對偶問題的解進行可行化，具體方式為：搜索沒有被分配的工作，選擇完成該工作消耗量最小的代理器為其服務。同理，搜索有多個代理器服務的工作，選擇完成該工作時消耗的資源量最小的為其服務，執行此操作直至所有工作都只有一個代理器為其服務。檢驗每個代理器消耗的總資源量，若超出自身占有，按照d的調整方式進行修正。可行化后的解為原問題的可行解，對應的目標函數值為原問題最優解的上界。

3.2 求解GAP問題的拉格朗日蝙蝠混合算法設計

將離散蝙蝠算法和拉格朗日松弛算法相結合，提出全新的解決GAP問題的拉格朗日蝙蝠算法（Lagrangian bat algorithm，LR-DBA）。兩種算法取長補短，離散蝙蝠算法求得GAP問題最優解的上界，此解為標準次梯度提供優秀上界；拉格朗日松弛算法求得GAP問題最優解的下界，當可行化后的解優于蝙蝠算法的最差解時，替代最差解；求得的上界、下界相互逼近，逐漸接近GAP問題的最優解。該算法的具體步驟如下：

步驟 1 初始化蝙蝠種群。蝙蝠種群個數s，初始音量Ai0，初始脈沖發生率ri0，最大迭代次數N，初始化蝙蝠個體位置xi。

步驟 2 蝙蝠速度和位置更新。計算得出最優蝙蝠x*。通過式（5）更新速度，利用單點遺傳基因換位算子更新位置xnew。

步驟 3 局部搜索。若rand＞ri，從最優解集中選一個解，利用單點遺傳基因換位算子產生一個局部解xnew。

步驟 4 若rand＜Ai且，則更新當前滿意解及其對應的目標函數值，并通過式（6）增大 ri，式（7）減小 Ai。

步驟 5 排列蝙蝠并找到當前最優x*，計算GAP問題最優解的上界，記為UB。

步驟 6 運用標準次梯度算法求解拉格朗日對偶問題，求得的解為GAP問題最優解的下界，記為LB。

步驟 7 對拉格朗日對偶問題的解進行可行化，求得GAP問題最優解的上界，記為 U B′。若可行化后的解優于蝙蝠算法的最差解，則更新蝙蝠算法最差解。

步驟 8 N′=N′+1，若 N′＜N，則轉步驟2；否則，結束搜索并輸出結果。

4 仿真實驗與分析

為了驗證算法的有效性，測試了Beasely的OR-library中的16組標準測試算例。其中，包括12組小規模算例，每組小規模算例中包含5個小算例；4組較大規模算例，每組較大規模算例中包含 3個小算例。在 Intel（R） Core（TM）i5-3320 CPU@2.60 GHz 2.60 GHz內存，64 位Windows 8操作系統的環境下利用Matlab R2015a進行編程實現和數據處理。該混合算法的初始參數設置為：蝙蝠種群數量s=20，迭代次數N=500，初始音量ri0=0.25，初始脈沖發生率 ri0=0.5， α =β=0.9， π =2，拉格朗日乘子λj的初始值參照文獻[20]的取值方式，即 λj=max cij,j∈ J。

為了測試LR-DBA混合算法的求解性能，將其與DBA算法進行比較。表1和表2分別列出了小規模和較大規模算例的相關數據。表中，運行時間為500次迭代的總時間；算例的目標函數值為隨機運行5次，取最好一次為滿意解；G為本文求得的滿意解和文獻[21]的最優目標函數值之間的差值；平均運行時間和平均G為一組算例中5個小算例的平均值；達到最優算例個數為一組算例中5個小算例中達到最優解的個數。C為目標函數值，D為最優目標函數值。G的計算式為

由表1可以看出，DBA算法求解的60個算例中只有2個達到最優解，而用LR-DBA混合算法求解時有10個達到最優，且12組算例的G均小于1%，幾乎達到最優解，表明LR-DBA混合算法在求解小規模算例中結果較好。表2中，對比DBA算法和LR-DBA混合算法的G，很明顯看出LR-DBA混合算法在求解較大規模GAP算例中也有很大的優勢，LR-DBA混合算法求解的12個算例中除gapc3以外，其他算例的G均在5%以內，而DBA算法求解所得G最大達到15%，且12個算例中有7個G都大于5%。分析全部算例，可以看出，LR-DBA混合算法在求解GAP問題算例時具有更好的全局搜索能力，且收斂精度較好。

表1 小規模GAP算例計算結果Tab.1 Calculating results of a small-scale GAP

表2 較大規模GAP算例計算結果Tab.2 Calculating results of a large-scale GAP

綜合分析表1和表2的實驗數據，LR-DBA混合算法求解GAP問題性能優異，尤其針對小規模算例。但由于LR-DBA混合算法結合了拉格朗日松弛算法，在求解小規模GAP問題時，LR-DBA混合算法的平均運行時間幾乎是DBA算法運行時間的8倍，在較大規模的求解中，LR-DBA的平均運行時間也是DBA算法的6倍，說明混合算法的收斂速度慢。

為了更直觀地展現LR-DBA混合算法的求解效果，本文以gapa3為例，給出求解GAP問題的迭代過程圖，如圖2所示。

圖2 算例gapa3的迭代過程圖Fig.2 Iterative process of the numerical example gapa3

由圖2可以看出，混合算法在求解過程中迭代的前期收斂速度快，能在較小的迭代次數得到一個較好解，但后期收斂速度較慢。綜合比較迭代次數、目標函數值和運行時間，可以看出，拉格朗日松弛算法一方面擴展了蝙蝠算法的搜索范圍，但另一方面也使得混合算法的運行時間變長。

5 結束語

應用BA算法來解決基本GAP問題，并針對GAP問題的特點，在基本BA算法框架基礎上引入遺傳算法單基因換位算子，提出DBA算法；在此基礎上，將DBA算法與LR算法結合組成LR-DBA混合算法，具有較強的全局搜索能力和魯棒性。通過大量的仿真實驗結果，表明LR-DBA混合算法求解GAP問題的性能優于DBA算法。

本文不僅為廣義分配問題提供了一種優異的求解方法，同時也拓寬了蝙蝠算法在離散問題上的進一步應用。但是，較現有的求解GAP問題的其他智能優化算法（如禁忌搜索和蟻群混合算法[11]、混合遺傳算法[22]、禁忌搜索算法[23]），LR-DBA混合算法的性能較差，并且如仿真實驗部分所示，LR-DBA混合算法在求解效率上也有待提高，這將是本文今后進一步研究的內容。