改進的螢火蟲算法求解阻塞流水線調度問題

郭麗萍，李向濤，谷文祥，2，殷明浩

(1.東北師范大學計算機科學與信息技術學院，吉林長春130117;2.長春建筑學院基礎教學部，吉林長春130607)

流水線調度問題是一類重要的組合優化問題，阻塞流水線調度問題作為其中的一個重要分支，由于具有很強的工程背景和實際意義，而得到了越來越多學者的關注，目前人們已經證明包含2臺機器以上的阻塞流水線調度問題的計算復雜性為NP-hard［1］.該問題的解決方法主要分為構造性啟發式方法［2-3］和元啟發式方法［4-7］.螢火蟲算法［8］是由劍橋大學的Yang在2008年提出來的一個模擬自然界中螢火蟲發光行為的仿生算法，目前該算法在函數優化方面表現出較優的性能［9］.為了使螢火蟲算法適用于求解阻塞流水線調度問題，本文對螢火蟲算法進行了一些改進，加強了算法的搜索性能，且避免了螢火蟲算法早熟的問題.

1 阻塞流水線調度問題

阻塞流水線調度問題可以描述為:n個作業J1，J2，…，Jn要在 m 臺機器 M1，M2，…，Mm上運行，每個作業包含m道工序，每個作業的工序必須按照統一順序在m臺機器上執行，每個作業在每臺機器上的運行時間是固定的，要求協調不同作業的執行順序，使得某種生產性能指標達到最優.阻塞流水線調度問題的約束條件包括如下:

1)每個作業在每臺機器上只能加工1次;

2)每臺機器每次最多只能加工1道工序;

3)每臺機器上的作業加工順序相同;

4)作業的某道工序在某臺機器上的加工一旦開始，便不能終止;

5)一個作業完成某道工序后，該作業將在當前機器上滯留直到下游機器變為可用為止.

在該問題中，用 π ={π1，π2，…，πn}表示一個調度序列，oj，k(j=1，2，…，n;k=1，2，…，m)表示作業j在機器k上執行一個無間斷的工序，其所需時間為 pj，k(j=1，2，…，n;k=1，2，…，m)，ej，k(j=1，2，…，n;k=0，1，…，m)表示第j個作業在第k臺機器上的離開時間.則有:

本文選取最大完工時間最小為優化指標，那么最大完工時間Cmax(π)=en，m，且它的計算復雜性為O(mn).

以3個作業J1、J2、J3為例，下面舉例說明如何求解阻塞流水線調度問題，3個作業在3臺機器上運行所需時間如表1所示.第1個作業在第1臺機器上運行所需要的時間為3，第1個作業在第2臺機器上運行所需要的時間為4，依此類推.

表1 時間矩陣Table 1 Time matrix of an example

由式(1)～(5)可以計算出該問題的最大完工時間的最小值 e3，3為21，對應的作業執行順序為{J3，J2，J1}，如圖 1 所示.

圖1 問題的甘特圖表示Fig.1 The Gantt chart of the example

2 螢火蟲算法

螢火蟲算法是一種通過模擬自然界中螢火蟲發光行為的仿生算法.目前，人們對自然界中螢火蟲的發光行為的意義還存在很多爭議，不過有2點是確定的:吸引配偶和吸引潛在的獵物.

一般情況下，可見光的強度隨著當前位置到光源距離的增大而減小，另外，空氣還要吸收一部分光源.因此，螢火蟲發出光只在一定的距離范圍內可以被其他螢火蟲看見.螢火蟲算法約定［9］:1)螢火蟲之間不存在性別之分，每只螢火蟲都可以吸引其他的螢火蟲，也可以被其他的螢火蟲吸引;2)螢火蟲的吸引力和它們發出的光的亮度成正比，因此，對于任意2只螢火蟲來說，亮度較弱的螢火蟲會朝著亮度較強的螢火蟲飛行;如果當前可見光距離范圍內，沒有比自己更亮的螢火蟲時，則該螢火蟲實行隨機飛行策略;3)光的亮度一般由求解問題的目標函數決定.

在螢火蟲算法中，每只螢火蟲稱為一個“個體”，每個個體包含亮度、吸引力、位置等屬性.在該算法中，個體在某個固定位置的亮度I是固定的，該變量的設定取決于具體問題的目標函數.其他個體看到該個體的亮度則隨著它們之間距離的增大而減小，此外，傳播介質也會吸收一部分光線，因此，個體在距離當前位置r處的亮度還和介質的吸收率有關［9］.一般把一個個體在距離該個體r處的亮度表示為

式中:I0為個體在當前所處位置的亮度，γ為介質的吸收率.有時，為了減小個體亮度變弱的速率，式(6)也可用如下公式代替:

每個個體對其他個體的吸引力β也是相對的，它隨著2個個體之間距離的增大而減小;此外，吸引力還和介質的吸收率有關［9］.因此把個體在距離該個體為r的位置對其他個體的吸引力定義為

式中:β0表示2個個體距離為0時，當前個體對另外1個個體的吸引力，γ為介質的吸收率.

在該算法中，2個個體之間的距離r采用歐式距離.個體i向著比它更亮的個體j進行更新的公式定義為

式中:xi和xj分別表示個體i和個體j在整個種群更新之前的位置，x'i表示個體i朝著比自己更亮的個體j更新之后的位置，rij代表個體i和個體j之間的距離，β0e－γr2ij表示個體 j對個體 i的吸引力，α 是一個［0，1］內的隨機數，R為一個隨機數生成函數，生成［0，1］內的隨機數.

3 改進的螢火蟲算法求解阻塞流水線調度問題

雖然螢火蟲算法在函數優化方面表現出了較好的計算性能，但是，它在求解阻塞流水線調度問題中仍然存在早熟現象，為了提高算法的尋優性能，本文對螢火蟲算法進行了一些改進.

螢火蟲算法應用于求解阻塞流水線調度問題時，首先遇到的問題就是如何將實數編碼離散化.本文提出一種最小排序方法，這種方法可以將個體的實數編碼形式轉化成離散的調度序列.該方法對初始化后的實數序列從小到大進行排序，它們的排序序號構成一個離散序列，將這個序列作為該實數序列對應的離散調度序列.如表2所示，0.43是實數編碼序列中最小的實數，從小到大排序后，它的序號應該為1;同理，0.91是實數序列中最大的，因此其序號為5.新生產的序列{4，2，5，1，3}作為實數序列{0.85，0.56，0.91，0.43，0.76}對應的離散調度序列.

表2 個體表示形式Table 2 The presentation of an object

其次，在種群初始化時，本文設計了一種雙重初始化方法:分別生成2個種群，計算2個種群中每個個體的最小完工時間.讓2個種群中的個體按序號兩兩進行比較:種群1中的第1個個體和種群2中的第1個個體進行比較，選擇其中最小完工時間較小的個體，作為新種群的第1個個體;種群1中的第2個個體和種群2中的第2個個體進行比較，選擇其中最小完工時間較小的一個，作為新種群的第2個個體;依此類推，新種群作為初始化種群.此外，由于 NEH(Nawaz-Enscore-Ham)［10］啟發式算法是目前求解流水線調度最有效的啟發式方法之一，因此，在初始化種群時還引入了NEH的思想.種群的第1個個體由NEH啟發式方法生成，其他個體由雙重初始化方式生成，從而使得算法有一個較好的初始化環境.

目前，人們已經發現很多昆蟲存在萊維飛行(Lévy flights)［11-12］，且萊維飛行已經應用于優化領域，并取得了預期的良好效果.為了增強算法的全局搜索性能，避免種群陷入局部最優，在螢火蟲算法的搜索過程中，如果當前個體找不到比自己更優的個體時，選擇萊維飛行代替原算法中的隨機飛行.此外，對于那些種群中非最優的個體，對它們的飛行公式進行改進:當它們發現比自己更亮的個體時，首先由系統產生一個隨機數q，如果q小于0.5，則按照式(8)進行更新;否則，仍采用式(7)對個體的位置進行更新.

式中:xj仍然表示個體j在整個種群更新之前的位置，x'i表示的是第i個個體朝著前j－1個體中比自己亮的個體更新之后的新位置，x″i表示x'i朝著比自己亮的個體j更新之后的位置，rij代表個體i和個體j之間的距離，β0e－γr2ij表示個體 j對個體 i的吸引力.可以看出，式(8)是實時更新的，式(7)僅僅依賴于整個種群移動之前的個體位置.這樣的飛行更新方式能夠增加飛行的隨機性，有利于保持種群的多樣性，增大種群搜索空間.

為了加速種群的收斂，本文還提出了一種α的更新方式，使得α隨著迭代次數逐漸減小，更新公式為

式中:α0選用0.9，t為迭代次數.

對于一個調度序列而言，目前主要有3種產生新序列的方法:插入、交換和逆序［6］.其中，插入方法被證明為最適合產生一個新序列的方法.為了加強螢火蟲算法的局部搜索能力，本文引入了一個局部搜索算法.局部搜索算法可以描述如下［7］:

2)i=i+1，如果 i＞n，則 i=i－n.

3)把πr賦給中間序列U，把序列U中的第i個作業從U中去除，將插入U中所有可以插入的位置，選取其中最好的調度序列πbest.

4)分別求解序列πbest和序列U的最小完工時間 f(πbest)和 f(U)，如果 f(πbest)＜f(U)，則 U=πbest，j=0;否則 j=j+1.

5)如果 j≤n，轉步驟2);否則，算法結束，πr=U.

在這個算法中，i和j用來控制循環次數，沒有實際意義，n是作業的數目，也是解的維數.

這種局部搜索算法的好處在于:首先，局部搜索策略并沒有直接對所選序列進行修改，而是應用于一個中間序列U，避免算法陷入局部最優;其次，在每次搜索過程中，只有找到更好的解時，才對中間序列進行更新，有利于算法逐步向更優的解靠近.因此，局部搜索算法在一定程度上加強了螢火蟲算法的局部搜索能力.

改進后的螢火蟲算法和原算法相比，初始化部分使得算法有一個較好的搜索域，局部搜索算法又進一步提高了算法的搜索性能.改進后的螢火蟲算法求解阻塞流水線調度問題的步驟如下:

1)初始化種群，包括迭代次數t=0、最大迭代次數tmax、隨機參數 α0、個體吸引力 β0、介質吸收率γ、局部搜索概率p.

2)按照式(1)～(5)計算每個個體的最小完工時間.

3)螢火蟲算法.

①更新α;

②對當前個體，如果有比自己更亮的個體，則按照改進后的更新方式，利用式(7)或(8)對當前個體的位置進行更新;否則，采用萊維飛行對個體位置進行更新.

4)按照式(1)～(5)計算每個個體的最小完工時間.

5)系統產生一個隨機數，如果這個隨機數小于局部搜索概率p，則對種群個體進行局部搜索，并更新種群.

6)t=t+1，如果t≥tmax或算法達到其他標準，則算法終止;否則，轉步驟3).

4 仿真實驗

實驗所用PC機的處理器為AMD Athlon64 X2 3600+1.91 GHz，內存為 960 MB，編程工具采用Matlab 7.0.

4.1 數據集

本文采用著名的Taillard數據集作為測試實例，該數據集包含12個子集，每個子集包含10個測試實例.為了檢驗算法的有效性，從這些子集中選取了部分具有代表性的實例作為測試數據，所選實例涵蓋了從20個作業5臺機器到200個作業10臺機器不同難度的測試實例，具有很好的代表性.所選實例的規模如表3所示.例如在Ta04實例中，20×5表示20個作業在5臺機器上運行.

表3 實例規模Table 3 The size of instances

4.2 實驗結果及分析

實驗中，種群大小為10，最大迭代次數tmax為100，α0、β0、γ 的設置參考了螢火蟲算法相關文獻［9，13-14］中的參數設置，α0設置為 0.9，β0設置為1.0，γ 設置為 0.9，局部搜索概率 p 設置為 0.2.對于前6個實例，每個實例進行10次獨立實驗;對于后4個實例，由于數據量較大，每個個體進行5次獨立實驗.為了測試改進后的螢火蟲算法在阻塞流水線調度問題中的性能，還對同等條件下的未改進的螢火蟲算法和標準的粒子群算法進行了測試.表4為未改進的螢火蟲算法的測試結果，表5為標準粒子群算法的測試結果，改進的螢火蟲算法的測試結果列于表6中.其中，在標準的粒子群算法中，2個學習因子c1和c2取值為2，慣性權重 ω設置為 0.7.

由表4、5可以看出，無論是最小值、最大值，還是平均值，螢火蟲算法的求解結果普遍優于標準粒子群算法.對比表4和表6，可以發現，改進的螢火蟲算法在求解效果上要遠遠優于基本的螢火蟲算法，而且求解結果更加穩定.綜合表4～6，改進后的螢火蟲算法的求解效果明顯優于基本的螢火蟲算法和標準的粒子群算法.并且隨著種群的增大，這種差距更加明顯，充分體現了算法的有效性，而且新算法的求解效果更加穩定.

表4 螢火蟲算法測試結果Table 4 The experiment results of firefly algorithm

表5 粒子群算法求解結果Table 5 The experiment results of particle swarm algorithm

表6 改進的螢火蟲算法求解結果Table 6 The experiment results of improved firefly algorithm

5 結束語

本文提出了一種改進的螢火蟲算法，設計了雙重初始化方法，引入了NEH啟發式方法，重新設計了算法中個體位置的更新方式，并引入了萊維飛行來增大種群的搜索域.在求解阻塞流水線調度問題時，設計了一種將實數編碼離散化的機制，實現了實數編碼算法對離散化問題的求解.此外，本文還引入了一種局部搜索算法來加強算法的局部搜索能力.實驗結果表明，改進后的螢火蟲算法在求解阻塞流水線調度問題時，性能有了明顯提高，而且隨著問題規模的增大，這種提高更加明顯，體現了算法的有效性和魯棒性.

由于阻塞流水線調度問題具有很重要的實際應用價值，因此將來的工作會考慮把其他優秀的算法引入到這個問題中進行求解.另外，由于元啟發式算法普遍存在計算量大的問題，將來的工作還應該考慮如何提高求解效率.

［1］ALLAHYERDI A，NG C T，CHENG T C E，et al.A survey of scheduling problems with setup times or costs［J］.European Journal of Operational Research，2008，187(3):985-1032.

［2］MCCORMICK S T，PINEDO M L，SHENKER S，et al.Sequencing in an assembly line with blocking to minimize cycle time［J］.Operations Research，1989，37(6):925-935.

［3］RONCONI D P.A note on constructive heuristics for the flowshop problem with blocking［J］.International Journal of Production Economics，2004，87(1):39-48.

［4］CARAFFA V，IANES S，BAGCHI T P，et al.Minimizing makespan in a blocking flowshop using genetic algorithms［J］.International Journal of Production Economics，2001，70(2):101-115.

［5］RONCONI D P.A branch-and-bound algorithm to minimize the makespan in a flowshop with blocking［J］.Annals of Operations Research，2005，138(1):53-65.

［6］GRABOWSKI J，PEMPERA J.The permutation flow shop problem with blocking.A tabu search approach［J］.Omega，2007，35(3):302-311.

［7］WANG Ling，PAN Quanke，SUGANTHAN P N，et al.A novel hybrid discrete differential evolution algorithm for blocking flow shop scheduling problems［J］.Computers ＆ Operations Research，2010，37(3):509-520.

［8］YANG Xinshe.Nature-inspired metaheuristic algorithms［M］.Frome，UK:Luniver Press，2008:81-96.

［9］YANG Xinshe.Firefly algorithms for multimodal optimization［C］//Proceedings of the 5th International Conference on Stochastic Algorithms:Foundations and Applications.Berlin/Heidelberg，Germany:Springer-Verlag，2009:169-178.

［10］NAWAZ M，ENSCORE E E J，HAM I.A heuristic algorithm for the m-machine，n-job flow shop sequencing problem［J］.Omega，1983，11(1):91-95.

［11］BROWN C T，LIEBOVITCH L S，GLENDON R.Lévy flights in Dobe Ju/＇hoansi foraging patterns［J］.Human Ecology，2007，35(1):129-138.

［12］PAVLYUKEVICH I.Lévy flights，non-local search and simulated annealing［J］.Journal of Computational Physics，2007，226(2):1830-1844.

［13］YANG Xinshe.Firefly algorithm，Lévy flights and global optimization［C］//Research and Development in Intelligent Systems XXVI.London，UK:Springer，2010:209-218.

［14］YANG Xinshe.Firefly algorithm，stochastic test functions and design optimization［J］.International Journal of Bio-Inspired Computation，2010，2(2):78-84.