張 建
(中國科學(xué)技術(shù)大學(xué) 管理學(xué)院,合肥 230026)
在工業(yè)生產(chǎn)中,集成電路芯片的制造包括生產(chǎn)、探測、裝配和測試環(huán)節(jié)[1-3].在測試過程中,需要將半導(dǎo)體芯片置于高溫烤箱中進(jìn)行耐溫測試以檢測出潛在的風(fēng)險.為了提高測試效率,芯片被分成不同的批次集中放在烤箱中處理.這里,半導(dǎo)體芯片就類似于工件,烤箱類似于機(jī)器.這是一個典型的批處理機(jī)調(diào)度的例子.
批調(diào)度問題不同于傳統(tǒng)的機(jī)器加工問題,一個批處理機(jī)能同時加工處理多個工件,從而有效提高資源的利用率.研究批調(diào)度問題,通過優(yōu)化調(diào)度有效提高生產(chǎn)效率,降低生產(chǎn)成本,對于提高生產(chǎn)管理水平、獲取更高的經(jīng)濟(jì)效益有著重要的現(xiàn)實意義[4].
人們對批調(diào)度問題的研究由來已久.Gimple 和Ikura[5]最先研究了工件尺寸相同、到達(dá)時間和交貨期相一致的單機(jī)批調(diào)度問題.Lee 和Uzsoy[6]提出了動態(tài)規(guī)劃的方法來解決此類調(diào)度問題,優(yōu)化目標(biāo)是最小化最大延遲時間和最小化延遲工件數(shù),局限性在于這種基于動態(tài)規(guī)劃的精確求解方式需要較長的計算時間.Uzsoy[7]首次將工件尺寸相同的性質(zhì)擴(kuò)展到差異工件.在此基礎(chǔ)上,Zhou[8]融入了工件動態(tài)到達(dá)的限制,并結(jié)合距離矩陣提出了構(gòu)造性啟發(fā)式算法,這些一次性啟發(fā)式算法在處理大規(guī)模工件問題時雖然在運行時間上占優(yōu),但解的質(zhì)量并不太好.
在實際生活中,工廠往往采用多臺機(jī)器并行處理的方式來提高生產(chǎn)效率.Chang[9]、Damodaran 和Chang[10]等考察了并行批處理機(jī)調(diào)度問題,優(yōu)化目標(biāo)是最小化制造跨度.其中Chang 采用模擬退火算法,通過交換初始解中不同批中的工件生成鄰域解,這導(dǎo)致初始解的質(zhì)量對最終解的質(zhì)量有很大影響,從而影響算法的性能.Chen[11]在并行機(jī)基礎(chǔ)上納入了動態(tài)到達(dá)的約束,并結(jié)合蟻群算法和遺傳算法以及ERT-LPT 規(guī)則來優(yōu)化完工時間,并且取得了不錯的效果.
元啟發(fā)式算法也被越來越多的學(xué)者應(yīng)用到調(diào)度問題中.例如基于構(gòu)建性的蟻群算法[12],基于鄰域搜索的模擬退火算法[3,9]和進(jìn)化算法,其中進(jìn)化算法包含但不限于遺傳算法[13]、粒子群算法[14]、差分進(jìn)化算法[15]等.分布式估計算法作為一種較新的基于種群的進(jìn)化算法,在求解組合優(yōu)化問題方面具有良好的表現(xiàn),比如置換流水車間問題[16],護(hù)士排班問題[17],但是在并行機(jī)調(diào)度方面的應(yīng)用研究卻很有限.因此,本論文采用分布式估計算法求解考慮差異工件的并行批處理機(jī)調(diào)度問題.
論文組織如下:第1 章節(jié)為問題描述;第2 章節(jié)詳細(xì)介紹分布式估計算法在考慮差異工件的并行批處理機(jī)調(diào)度中的應(yīng)用;第3 章節(jié)為仿真實驗;第4 章節(jié)為結(jié)論與展望.
論文考慮包含差異工件的并行批處理機(jī)調(diào)度問題,相關(guān)假設(shè)如下:
(1)將n個工件分配到m個批處理機(jī)加工;
(2)工件具有不同的尺寸和加工時間;
(3)各臺機(jī)器容量一致,工件j只需在其中一臺機(jī)器上加工,且它在所有機(jī)器上加工時間相同;
(4)工件集中分批,按照批次順序安排在機(jī)器上加工;
(5)批的尺寸等于批中所有工件尺寸之和,且不能超過機(jī)器容量限制,批的加工時間等于批中工件最大加工時間;
(6)一旦某個批次開始加工,加工過程不能中斷且不能有新的工件加入到當(dāng)前批中;
(7)優(yōu)化目標(biāo)是makespan,也就是使最后離開加工系統(tǒng)的批次的完工時間最小;
該問題采用調(diào)度理論中的三元組法[18]可表示為Pm|batch,B,sj|Cmax,以下給出相關(guān)參數(shù)定義及混合整數(shù)規(guī)劃模型:
① 集合
工件集合:{j|j=1,2,3,···,n}
批集合:{b|b=1,2,3,···;b≤n}
機(jī)器集合:{k|k=1,2,3,···,m}
② 參數(shù)
sj:工件j尺寸;
pj:工件j加工時間;
B:機(jī)器容量;
③ 決策變量
PTbk:在機(jī)器k上加工的批b的加工時間;
CTbk:批b在機(jī)器k上的完工時間;
Cmax:最終完工時間;
xjbk:工件j分在批b中安排在機(jī)器k上加工時為1,否則為0;
數(shù)學(xué)模型如下:
式(1)為優(yōu)化目標(biāo);式(2)保證一個工件只能分配在一個批中只能在一個機(jī)器上加工;式(3)保證批的容量不超過機(jī)器的容量限制;式(4)說明批的加工時間為批中最長加工時間;式(5)表明前后批完工時間之間的關(guān)系;式(6)為Cmax的計算方法;式(7)是對變量xjbk的取值約束.
為了能更有效評估后續(xù)算法中解的質(zhì)量,我們考慮這樣一個下界:假設(shè)存在某種單位工件ej具有單位尺寸單位工時,一個尺寸為sj、加工時間為pj的普通工件可表示為sj*pj個ej,此時問題轉(zhuǎn)化成若干個性質(zhì)相同的單位工件的并行機(jī)調(diào)度問題.理想狀態(tài)下所有單位工件能完美分配到批中,即批中所有工件尺寸之和恰好等于批的容量,然后均勻分配到各臺機(jī)器上加工,可得:
顯然這樣計算得出的結(jié)果是一個有效的下界,其中批的利用率達(dá)到最高.
分布式估計算法是一種基于概率模型的進(jìn)化算法,最早于1996年由Mühlenbein,Bendisch[19]等人提出,是當(dāng)前國際進(jìn)化計算領(lǐng)域的研究熱點.根據(jù)變量之間的相關(guān)性,算法可分為三類:變量無關(guān)、雙變量相關(guān)和多變量相關(guān);按照變量的類型可分為離散變量的分布式估計算法和連續(xù)域上的分布式估計算法.
分布式估計算法不同于經(jīng)典的遺傳算法,它摒棄了遺傳算法中采用選擇、交叉、變異的方式生成子代種群的策略,轉(zhuǎn)而采用概率模型來更新下一代種群,通過收集統(tǒng)計全局信息,學(xué)習(xí)優(yōu)勢個體特征來更新概率模型進(jìn)行采樣,然后生成新的個體.該算法基本流程如圖1所示.
論文采用基于工件序列的編碼方式,每個工件有相應(yīng)的序列編號,根據(jù)不同的排列可得到不同的調(diào)度方案.
圖1 EDA 基本流程圖
給定一個工件序列,首先我們采用FF (First-Fit)規(guī)則進(jìn)行分批,然后根據(jù)LPT (Longerst Processing Time)方法選擇當(dāng)前處于空閑狀態(tài)的機(jī)器進(jìn)行調(diào)度.FF 規(guī)則如下:
步驟1:選擇工件序列中第一個工件加入到當(dāng)前批中,若違反機(jī)器容量約束,考慮下一個工件;
步驟2:每當(dāng)工件完成分批,則從序列中去除該工件;若當(dāng)前批容量已滿或不能容納其它任何工件,則新建一個批;
步驟3:重復(fù)步驟1~2,直到所有工件完成分批;
分批完成后,按批次的加工時間由長到短重新排列,并根據(jù)此排列順序選擇非加工狀態(tài)的機(jī)器進(jìn)行加工,若存在多臺非加工狀態(tài)的機(jī)器,則任選其中之一進(jìn)行加工.
考慮10 個工件在兩臺機(jī)器上加工,假設(shè)各臺機(jī)器容量均為15,工件序列為[4,5,1,3,6,2,9,10,7,8],工件屬性及分批調(diào)度結(jié)果如圖2所示.
算法采用概率矩陣P來表征分布估計模型,Pij(t)表示第t代中工件i出現(xiàn)在位置j周圍的概率,這里用“周圍”是因為我們考慮了幾種不同的更新策略來計算工件i在位置j左右的概率,而不是僅僅在固定位置j上.
圖2 10 工件調(diào)度示意圖
算法初始階段所有概率值均設(shè)為1/n,以保證工件能被均勻抽樣.對于每一代種群,設(shè)種群規(guī)模為Q,我們根據(jù)個體的Cmax值計算適應(yīng)度值(fitness),然后挑選最優(yōu)SP_Size個個體估計其概率分布,通過對其學(xué)習(xí)來更新概率矩陣.這里的fitness及SP_Size計算如下:
其中,α為優(yōu)勢個體選擇比例,0<α<1.
論文中采用FF-LPT 規(guī)則解碼,因此工件的位置對于批的形成有很大影響,相鄰位置的工件就越有可能分配到同一個批中.我們考慮4 種更新機(jī)制來完善概率矩陣,第一種更新機(jī)制學(xué)習(xí)優(yōu)勢個體中當(dāng)前位置工件的概率分布;第二種更新機(jī)制學(xué)習(xí)當(dāng)前位置及其之前位置工件的概率分布;第三種更新機(jī)制與第二種相反,通過學(xué)習(xí)優(yōu)勢個體中當(dāng)前位置及其之后位置工件的概率分布;第四種更新機(jī)制考慮當(dāng)前位置及其鄰域內(nèi)工件的概率分布.具體更新方式如下:
其中,β為學(xué)習(xí)率,0<β<1;vj為單邊鄰域大小,例如當(dāng)vj= 2 時,位置5 的鄰域為位置3,4,5,6,7,|v5|為5;位置1 的鄰域為位置1,2,3,|v1|為3.
為了闡明概率矩陣更新過程,我們用一個包含5 個工件的算例示范.假設(shè)挑選4 個優(yōu)勢個體,如圖3所示.設(shè)vj為1,根據(jù)4 種更新機(jī)制統(tǒng)計各個位置的工件信息,構(gòu)建優(yōu)勢個體的概率分布矩陣,結(jié)果如圖4所示.通過學(xué)習(xí)優(yōu)勢個體的概率分布特征來更新概率矩陣P,然后基于比例選擇采用輪盤賭的方式依次決定各個位置上工件編號,即可獲得新的個體.
圖3 挑選的4 個優(yōu)勢個體
綜上,面向差異工件的并行批處理機(jī)調(diào)度問題的分布式估計算法如算法1 所示.
算法1.面向差異工件的并行批處理機(jī)調(diào)度問題的分布式估計算法初始化參數(shù)Q,α,β,gen(迭代次數(shù)),隨機(jī)生成初始種群;t= 0 ;Whilet<gen采用FF-LPT 規(guī)則對種群中個體進(jìn)行解碼,計算適應(yīng)度值;根據(jù)適應(yīng)度值選取最優(yōu)的SP_Size個個體,并估計其概率分布;選取2.3 中的更新策略對概率矩陣進(jìn)行更新;根據(jù)新的概率矩陣采樣生成新的種群;令t=t+1;end輸出最優(yōu)調(diào)度方案及makespan.
本實驗根據(jù)文獻(xiàn)[9]提出的方法生成隨機(jī)算例,算例規(guī)模按3 個標(biāo)準(zhǔn)劃分:工件數(shù)量n,工件尺寸s,加工時間p;其中s、p均服從離散均勻分布;機(jī)器數(shù)量m為2 或4,容量均為20.具體參數(shù)見表1.
每類算例提供一個編號.例如包含50 個工件,工件尺寸服從[1,10]分布,加工時間滿足[1,20],在兩臺機(jī)器上加工的一類算例編號為J2S3P2M1.每類算例隨機(jī)生成10 個子算例,每個算例計算10 次.
圖4 不同更新機(jī)制下優(yōu)勢個體概率矩陣計算結(jié)果
表1 算例規(guī)模
影響EDA 性能的參數(shù)主要有4 個:Q,α,β,gen.其中,優(yōu)勢個體的選擇比例α對概率矩陣的更新有很大影響,若比例選擇過小,子代種群與當(dāng)代個體相似度較高,則容易陷入局部最優(yōu);若比例選擇過大,個體的優(yōu)勢特征弱化,導(dǎo)致概率矩陣的學(xué)習(xí)效果差.為此,論文考慮3 種不同的α值進(jìn)行預(yù)實驗.為了選取合適的參數(shù)值展開后續(xù)實驗,我們采用田口方法對不同水平的參數(shù)進(jìn)行評估,將參數(shù)作為因子,并對每個因子設(shè)置3 水平的參考值,具體的參數(shù)設(shè)置見表2.
表2 參數(shù)設(shè)置
算法提供了4 種更新策略,為方便表述,依次記為EDA1,EDA2,EDA3,EDA4,其中EDA3 表示采取第三種更新策略.由于EDA4 中需要考慮單邊鄰域v,因此為EDA4 單獨設(shè)立5 因素3 水平的實驗,參數(shù)順序為Q,α,β,v,gen.v具有3 水平,分別是[2,3,5],其余因素同表2.
察爾汗鹽湖是我國最大可溶性鉀鎂鹽礦床。1958年,第一批鹽湖人深入戈壁荒漠,正式開啟我國鉀肥工業(yè)生產(chǎn)序幕。60年來,我國鉀肥工業(yè)從無到有、從小到大、從弱到強(qiáng),形成了科研、設(shè)計、設(shè)備制造、施工安裝、生產(chǎn)、銷售、農(nóng)化服務(wù)等一套完整工業(yè)體系,產(chǎn)品市場競爭力不斷增強(qiáng),鉀肥生產(chǎn)技術(shù)和單套生產(chǎn)能力均達(dá)國際領(lǐng)先水平。
根據(jù)實驗參數(shù)和因子水平,選取正交陣列L9(34)和L27(35),其中L9(34)表示針對3 因素4 水平的問題設(shè)計9 種實驗.我們針對50 規(guī)模的工件,生成12(3×2×2)個算例,每個算例采取9 次試驗,算5 次,一共計算12×9×5 次;對于EDA4,需要計算12×27×5 次.最后的平均響應(yīng)變量值(ARV)計算方式如下:
其中,n為算例個數(shù).顯然,ARV 越小越好.
圖5中各子圖橫坐標(biāo)為因子水平值,縱坐標(biāo)為ARV,分析可知,EDA1 最優(yōu)參數(shù)組合為:Q=60,α=0.2,β=0.1,gen=500;EDA2 最優(yōu)組合為:Q=60,α=0.1,β=0.1,gen=500;EDA3 最優(yōu)組合為:Q=50,α=0.1,β=0.3,gen=500;EDA4 最優(yōu)組合為:Q=60,α=0.1,β=0.3,ν=2,gen=500.
為了評估算法的性能,我們將和文獻(xiàn)[9]提出的模擬退火算法(SA)以及經(jīng)典的遺傳算法(GA)作對比.SA 中各項參數(shù)與原文獻(xiàn)一致,GA 中種群規(guī)模和迭代次數(shù)同EDA 一樣,變異概率設(shè)為0.1.所有的算法均由Matlab 實現(xiàn),實驗環(huán)境為4 核3.4 GHz、8 GB RAM 的Windows10 系統(tǒng).
表3、表4分別給出了2 臺和4 臺機(jī)器環(huán)境下所有算例的結(jié)果.實驗中分別統(tǒng)計了各個算法運行得到的最優(yōu)值、平均值、最差值及運行時間等,這里只給出了比值ratio,計算公式為:
ratio越小,說明結(jié)果越接近下界,所得到的解越接近最優(yōu)解.
圖5 不同更新機(jī)制下EDA 的因子水平趨勢
從表3、表4中可以看出,對于絕大部分算例,我們提出的4 種EDA 得到的解均是要優(yōu)于SA 算法的.對于2 臺機(jī)器,EDA1 的ratio均值為1.24,GA 為1.27,SA 為1.44.隨著工件規(guī)模的增加,算法的解的質(zhì)量也越來越好,尤其是在較多并行機(jī)的情況下.
值得注意的是,4 種EDA 算法中,EDA1 的性能是最好的,其次是EDA2,EDA3 的性能最差.EDA2 與EDA3 的更新機(jī)制是截然相反的,EDA3 考慮當(dāng)前位置及其之后位置工件的概率,但是由于采樣時我們是按照順位比例選擇的方式依次決定各個位置的工件,這導(dǎo)致EDA3 效率較低.而EDA1 只考慮當(dāng)前位置的工件的概率分布,使得概率分布極為“精純”,所以性能較好.
另外,GA 的性能比想象中好很多,對于小規(guī)模算例,GA 的解甚至比EDA1 還要占優(yōu);但對于大規(guī)模工件,EDA1 得到的解的質(zhì)量要優(yōu)于GA.
表3 M1 類算例結(jié)果
表4 M2 類算例結(jié)果
運行時間方面,4 種EDA 均是要長于GA、SA 的.這是由于EDA 生成下一代種群需要反復(fù)對概率矩陣采樣生成新的個體,這一過程相對于SA 中鄰域搜索或GA 中交叉、變異要繁雜很多,因此造成時間上的冗余.
圖6給出了J2S2P1M1 類典型算例的收斂趨勢圖.從中可以看出EDA1 的收斂性更好.其他類算例的收斂趨勢與此類似.
圖6 J2S2P1M1 類典型算例收斂趨勢圖
本文研究了分布式估計算法在考慮差異工件的并行批處理機(jī)調(diào)度中的應(yīng)用問題.首先提出了一個混合整數(shù)規(guī)劃模型和一個下界,然后采用FF-LPT 規(guī)則實現(xiàn)對工件的分批和排序;對于基于概率矩陣的分布式估算法,提出了4 種不同的更新機(jī)制,并通過仿真實驗和SA、GA 對比,驗證了算法的有效性.經(jīng)比較知,對于采取不同更新機(jī)制的EDA,EDA1 的性能最好,EDA2 次之,EDA3 性能最差.
未來的研究可考慮通過局部優(yōu)化來進(jìn)一步優(yōu)化算法的性能,以及考慮多階段并行機(jī)調(diào)度問題或不同的優(yōu)化目標(biāo)等.