章 旸，王鋮愷，劉 冉

(上海交通大學(xué) 機械與動力工程學(xué)院，上海 200240)

0 引言

柔性流水車間(Flexible Flow Shop, FFS)常見于半導(dǎo)體制造、鋼鐵制造、汽車制造等大量工業(yè)環(huán)境。作業(yè)調(diào)度是流水車間生產(chǎn)過程中的重要工作和經(jīng)典問題，也是典型的NP-hard問題。合理的作業(yè)調(diào)度安排對于生產(chǎn)效率與資源利用均有很大的意義。對于可分批加工的生產(chǎn)任務(wù)，批量流方法可有效縮短最大加工完成時間。將每種工件的總加工批量劃分為若干個子批，其中每個子批包含若干個單工件[1]，通過對子批的單獨調(diào)度，實現(xiàn)同種工件在多個工序中的并行加工，可縮短整個生產(chǎn)任務(wù)的完成時間。這類調(diào)度問題一般稱為批量流調(diào)度問題。

批量流調(diào)度能夠有效縮短生產(chǎn)任務(wù)的最大完成時間，但這種模式下的生產(chǎn)調(diào)度也更加困難。一般工件的分批越細(xì)，最大完成時間的縮短情況越明顯，但過多的子批數(shù)會增加問題的規(guī)模和調(diào)度難度。同時，在實際生產(chǎn)中，每臺機器在子批切換過程中需要進行換模，換模時間一般與加工順序密切相關(guān)，這進一步加大了問題的復(fù)雜性。而且，隨著生產(chǎn)環(huán)境的復(fù)雜化與不確定性的提高，生產(chǎn)過程的隨機性也越發(fā)重要且必須被調(diào)度所考慮，例如工件的隨機來料，工件在機器上的加工時間隨機以及機器故障等?？紤]隨機性的車間調(diào)度問題已經(jīng)非常復(fù)雜，考慮隨機因素的批量流調(diào)度顯然更加困難。

柔性流水車間內(nèi)批量流調(diào)度問題(Flexible Job-shop Scheduling Problen with Lot Streaming, FJSP-LS)已經(jīng)引起學(xué)界的重視和研究。LOW等[2]論述了車間調(diào)度中應(yīng)用批量流方法的優(yōu)勢。對于FJSP-LS的求解，根據(jù)批量劃分方式的不同，可分為均等子批、一致子批和可變子批[3]。目前，大多數(shù)柔性車間分批調(diào)度研究中，多采用一致子批的劃分策略，白俊杰等[4]采用通過以最小化加工完成時間為目標(biāo)對子批加工順序進行分批與排序，指出在子批數(shù)量相同的情況下,柔性分批調(diào)度方法比均等分批調(diào)度方法可以更有效地縮短交貨期。由于可變子批的復(fù)雜性，目前在該方面的研究較少，WANG等[5]提出一種三階段方法來解決有換模時間的FJSP-LS，通過最小批量進行劃分后以加工完成時間最小為目標(biāo)確定調(diào)度方案，再優(yōu)化批量，加工過程中子批批量可變。但為了縮短加工時間，往往需要更多的子批數(shù)據(jù)從而增加傳輸次數(shù)，生產(chǎn)成本變大。根據(jù)求解方法可將FJSP-LS分為精確算法與啟發(fā)式算法兩個研究方向。其中精確算法主要用來求解規(guī)模較小的批量流調(diào)度問題，MORTEZAEI等[6]通過建立混合整數(shù)規(guī)劃模型實現(xiàn)了FJSP-LS的求解，并指出隨著加工階段的增多，批量流方法的效果更為顯著，但是該模型對于所求解的問題規(guī)模要求較高，最多為兩種工件，每種工件最多可分為3個子批;ALFIERI等[7]針對考慮換模時間的流水車間批量流調(diào)度問題提出一種動態(tài)規(guī)劃算法，可在短時間內(nèi)實現(xiàn)精確求解，但是該算法僅針對單種工件的兩階段加工生產(chǎn)線的分批調(diào)度，并不能應(yīng)對多種類、多生產(chǎn)階段的生產(chǎn)調(diào)度問題。近年來，為了應(yīng)對越來越復(fù)雜的制造環(huán)境，以及對于大規(guī)模問題求解的要求，各種啟發(fā)式算法被不斷引入到批量流生產(chǎn)調(diào)度問題的研究中，VENTURA等[8]、TSENG等[9]分別采用遺傳算法與粒子群算法進行求解，并且有大量學(xué)者通過對各啟發(fā)式算法的改進實現(xiàn)了批量流問題的高效求解，如ALFIERI等[10]提出一種基于動態(tài)規(guī)劃的精確與啟發(fā)式相結(jié)合的算法來求解單種工件、兩臺機器的批量流調(diào)度問題;劉雪紅等[11]針對柔性車間可變子批的特點建立了以最小化最大完成時間和最小化批次數(shù)目為優(yōu)化目標(biāo)的多目標(biāo)柔性作業(yè)車間調(diào)度模型，并設(shè)計改進的候鳥算法進行求解，其中禁忌搜索(Tabu-Search, TS)是求解FJSP的經(jīng)典算法[12-13]，且與各啟發(fā)式算法相比，TS能夠在合理的時間內(nèi)給出較好的調(diào)度方案。近年來，越來越多的學(xué)者將TS運用到FJSP-LS上，并獲得了高質(zhì)量的解。HAMED等[14]采用TS針對多種類工件的批量流調(diào)度問題進行求解，并設(shè)計了3種鄰域搜索來加速算法的收斂，取得了良好的效果;SHAHVARI等[15]將TS與混合整數(shù)規(guī)劃相結(jié)合，通過TS探索解空間，在加工路線確定后將混合整數(shù)規(guī)劃問題轉(zhuǎn)化為線性規(guī)劃對目標(biāo)解進行評估。

對于隨機條件下的車間調(diào)度問題，現(xiàn)有的相關(guān)研究并不多，主要的研究方向為考慮工件加工時間隨機、工件隨機來料或設(shè)備故障的條件下最小化總加工完成時間。其中：WANG等[16]提出一種基于兩階段仿真的分布評估算法，用于在工件加工時間隨機的情況下最小化流水車間的總加工時間；ZANDIEH等[17]針對工件隨機來料，考慮機器故障的車間調(diào)度問題引入了變領(lǐng)域搜索算法(Variable Neighborhood Search, VNS);FU等[19]研究了工件加工時間服從正態(tài)分布，最小化最大完成時間與拖期的多目標(biāo)優(yōu)化問題，建立了混合整數(shù)規(guī)劃模型，并開發(fā)了煙花算法來加速模型的求解。以上方法均僅考慮了工件的單獨調(diào)度，并未考慮工件的分批問題，且不能應(yīng)對復(fù)雜的隨機環(huán)境。HAN等[18]對工件加工時間隨機的批量流調(diào)度問題采用蒙特卡洛采樣方法，將問題轉(zhuǎn)化為確定的多目標(biāo)優(yōu)化問題，并提出一種多目標(biāo)遷徙鳥優(yōu)化算法進行求解，但是該問題并未考慮最小批量問題，即工件的子批批量可以為符合條件的任意整數(shù)，且不考慮批次切換過程中的換模時間。

Benders分解方法由BENDERS[20]首先提出，用于解決混合整數(shù)規(guī)劃問題，近年來也越來越多被用于大規(guī)?；蜉^為復(fù)雜的車間調(diào)度問題的求解。如HARJUNKOSKI等[21]提出一種求解鋼鐵行業(yè)大規(guī)模的復(fù)雜調(diào)度問題的分解算法;WALDHERR等[22]通過分解算法求解工件在加工過程中需要額外資源、有并行機器且考慮機器換模時間的復(fù)雜流水車間調(diào)度問題。Benders分解還可與SAA(sample average appropriation)方法結(jié)合，實現(xiàn)對隨機問題的求解[23-24]。

雖然目前對批量流問題和隨機車間調(diào)度問題均分別已有相關(guān)研究，但將兩者結(jié)合起來的研究較少，且在考慮隨機性時大多為針對某特定分布的研究，不能適應(yīng)更為復(fù)雜的隨機環(huán)境，并且目前還沒有用Benders分解來求解隨機環(huán)境下的FJSP-LS的相關(guān)文獻(xiàn)。本文研究了復(fù)雜隨機來料的可分批柔性生產(chǎn)線作業(yè)調(diào)度問題，基于SAA方法建立了混合整數(shù)規(guī)劃模型，該模型適用于復(fù)雜隨機環(huán)境，并考慮模型的難求解性以及Benders分解在求解難解的混合整數(shù)規(guī)劃問題上的優(yōu)越性，本文采用該方法來實現(xiàn)對復(fù)雜隨機環(huán)境下柔性車間批量流調(diào)度問題的求解，同時設(shè)計禁忌搜索算法實現(xiàn)近似求解，可應(yīng)對復(fù)雜隨機環(huán)境以及大規(guī)模的批量流調(diào)度問題的求解，通過數(shù)值實驗驗證了方法的有效性。

1 問題描述與數(shù)學(xué)模型

1.1 問題描述

本文研究對象是由K臺機器組成的生產(chǎn)線，其中機器位置固定且緩沖區(qū)容量無限，每臺機器可進行一種或多種工藝的加工，但同一時間僅可加工一個工件，所有機器在系統(tǒng)時刻0處于空閑狀態(tài)，不考慮機器故障。生產(chǎn)線需加工I種不同種類的工件，對于一個工件i∈I，其總加工任務(wù)量為di，di可被分為若干個子批，所有工件的來料時間不確定，服從某隨機分布。工件來料后方可進行加工，在進入生產(chǎn)前即分為若干個批量不等的子批進入機器進行加工，子批在每臺機器上加工時批量保持不變，每一子批加工前有與加工順序相關(guān)的機器換模時間。子批批量規(guī)定為最小批量minS的整數(shù)倍，同時規(guī)定di為minS的整數(shù)倍，工件i最多可被分為ui個子批，ui=di/minS。所有子批按照加工順序依次經(jīng)過所有機器，只有在前一個子批全部加工完成并離開機器后，下一子批才可進入該機器并進行換模。工件i在機器k上的加工時間確定，為PTik，機器k上工件i到工件j上的換模時間STijk與i，j相關(guān)，不考慮工件在機器間的運輸時間。

問題需確定每類工件應(yīng)分為多少個子批、每個子批包含多少個最小批量以及所有子批的加工順序，以最小化生產(chǎn)線在隨機環(huán)境下的最大完成時間期望。

如圖1所示：Cmax表示該調(diào)度方案的加工完成時間，J11表示工件1的第一個子批，J12表示工件1的第二子批，以此類推。其中：圖1a～圖1c為工件來料時間均為0時3種調(diào)度方式的加工完成時間比較，圖1a為不分批，圖1b與圖1c分別為兩種分批調(diào)度方案。由圖1a與圖1b、圖1c的對比可見，對工件進行分批操作有利于縮短工件的加工完成時間;而對比圖1b與圖1c可見，不同的分批調(diào)度策略縮短加工時間的效果明顯不同，分批1優(yōu)于分批2。當(dāng)工件來料時間發(fā)生變化時，如圖1d與圖1e所示，其中分批2的加工完成時間小于分批1。綜合對比圖1b、圖1c與圖1d、圖1e，由不同的對比結(jié)果可知工件的來料時間對調(diào)度方案的效果有較大影響。在不同的來料時間背景下，最優(yōu)的加工順序與分批方案不同，最優(yōu)的加工完成時間也不同，因此當(dāng)來料時間不確定時，難以確定最優(yōu)調(diào)度方案。

1.2 基于SAA的隨機規(guī)劃模型

在隨機環(huán)境中，工件來料時間、子批在各機器上的加工完成時間與生產(chǎn)線最大完成時間都具有隨機性質(zhì)，很難直接建立數(shù)學(xué)模型求解。SAA方法根據(jù)工件來料時間的隨機特性生成S個場景，目標(biāo)是最小化這S個場景的最大完成時間均值。

建立本問題隨機混合整數(shù)規(guī)劃模型MIPSAA如下：

(1)優(yōu)化變量

nij為整數(shù)變量，工件i的第j個子批的最小批量個數(shù)，可以為0;

Sij為整數(shù)變量，工件i的第j個子批的工件總數(shù)，可以為0;

zijgh為0-1變量，zijgh=1時，工件i的第j個子批和工件g的第h個子批均不為0，且工件i的第j個子批為工件g的第h個子批的前一個子批，否則zijgh=0;

Cs為連續(xù)變量，場景s的最大完成時間;

tijks為連續(xù)變量，場景s中工件i的第j個子批在機器k上的加工完成時間;

yij為0-1變量，nij>0時，yij=1;nij=0時，yij=0。

(2)隨機變量

ATis為場景s中工件i的來料時間。

(3)其他變量

uij為輔助變量，用于順序記錄。

MIPSAA為:

(1)

s.t.

tijks≥tghks+SijPTik+STgik-(1-zghij)M,

0≤i≤I,1≤j≤ui,0≤g≤I,

1≤h≤ug，1≤k≤K,1≤s≤S,

(i,j)≠(g,h);

(2)

tijks≥tij(k-1)s+SijPTik+STgik-(1-zghij)M,

0≤i≤I,1≤j≤ui,0≤g≤I,1≤h≤ug，

1≤k≤K,1≤s≤S,(i,j)≠(g,h);

(3)

yij≤nij≤yijM,0≤i≤I,1≤j≤ui;

(4)

Sij=nij·minS,0≤i≤I,1≤j≤ui;

(5)

(6)

Cs≥tijKs,1≤i≤I,1≤j≤ui,1≤s≤S;

(7)

tij0s≥ATis,1≤i≤I,1≤j≤ui,1≤s≤S;

(8)

(9)

(10)

(11)

uj-ugh+M·zghij≤M-1,

0≤i≤I,1≤j≤ui,0≤g≤I,1≤h≤ug;

(12)

y01=1;

(13)

yij,zijgh∈{0,1},tijks≥0,Cs≥0,

intSij≥0,intnij≥0,uij≥0，

0≤i≤I,1≤j≤ui,0≤g≤I,

1≤h≤ug,0≤k≤K,1≤s≤S。

(14)

其中：目標(biāo)函數(shù)(1)為最小化S個場景的最大完成時間均值;約束(2)表示同一機器上前后兩個子批之間的約束;約束(3)表示同一子批在相鄰兩臺機器上的加工約束，M為一個很大的正數(shù)，只有zijgh等于1時約束(2)、(3)才生效;約束(4)表示輔助變量yij與nij之間的關(guān)系;約束(5)表示每種工件每個子批包含Sij的工件數(shù)量為最小批量的整數(shù)倍;約束(6)表示同種工件的所有子批大小之和等于該種工件總批量;約束(7)表示最大完成時間不小于所有批次在最后一臺機器上完工時間;約束(8)表示所有子批在該種工件來料后才可加工;約束(9)與約束(10)分別約束了每一個批次的緊前與緊后工件最多只有一批;約束(11)表示虛擬工件0為第一個子批;約束(12)保證了由變量zijgh決定的加工順序是一條無閉環(huán)的完整的生產(chǎn)線;約束(13)表示虛擬工件0的第一個批次批量大于0;約束(14)表示各變量的取值范圍。

2 求解算法

MIPSAA是一個非常難以求解的混合整數(shù)規(guī)劃模型，車間調(diào)度問題(Job-shop Scheduling Problem, JSP)為經(jīng)典的NP-hard問題[25]，而本文研究的柔性車間調(diào)度還需考慮工件的分批以及工件的隨機來料，進一步增加了問題的求解難度，顯然也具有NP-hard的性質(zhì)。即使對很小的問題算例，MIPSAA包含大量的約束和決策變量，例如對于3種工件，每種工件最多可分3批，場景數(shù)為75的小規(guī)模問題包含的約束多達(dá)14萬條，變量6 000多個，直接通過CPLEX等求解工具無法在有限時間內(nèi)求解。本文引入Benders分解方法并設(shè)計了兩套算法來對該調(diào)度問題進行優(yōu)化。首先簡單介紹Benders分解算法，其核心思想是通過分離容易處理的變量與復(fù)雜變量，將原問題分解為子問題與主問題，對兩類變量分別求解。模型MIP為：

mincTx+dTy。

(15)

s.t.

Ax+By≥b;

(16)

x≥0,y∈Y。

(17)

(18)

(19)

x≥0。

(20)

(21)

s.t.

ATu≤c;

(22)

u≥0。

(23)

遍歷所有的y不現(xiàn)實，雖然MIP的可行域復(fù)雜，但在求解最優(yōu)解的過程中，僅需要有限的cut(割平面)去逼近原問題。cut構(gòu)造方法如下：

dTy+(b-By)Tup≤w。

(24)

式(24)稱為optimalitycut。

(2)若DSP無界，記ur為DSP的一條極線，根據(jù)對偶理論，SP不可行，故增加限制減小使DSP無界的區(qū)域，等價于縮小了使SP不可行的區(qū)域，增加的cut為：

(b-By)Tup≤0。

(25)

式(25)稱為feasibilitycut，與optimalitycut共同完成主問題的構(gòu)造。MIP可重寫為：

(26)

構(gòu)造僅含整數(shù)變量的主問題(Master Problem, MP)如下：

allthefeasibilitycuts}。

(27)

主問題為由cut組成的整數(shù)規(guī)劃，cut由DSP問題的極點與極線給出。Benders分解通過MP與DSP求解與迭代直至問題上下界相等，找到最優(yōu)解。

2.1 基于Benders分解的精確算法

2.1.1 主問題與子問題構(gòu)建

minCs。

(28)

s.t.

約束(7)～(8);

(29)

0≤i≤I,1≤j≤ui,0≤g≤I,

1≤h≤ug，1≤k≤K,(i,j)≠(g,h);

(30)

0≤i≤I,1≤j≤ui,0≤g≤I,1≤h≤ug，

1≤k≤K,(i,j)≠(g,h);

(31)

tijks≥0,Cs≥0，intSij≥0,0≤i≤I,

1≤j≤ui,0≤g≤I,1≤h≤ug,0≤k≤K。

(32)

SPSAA-s的對偶問題DSPSAA-s如下，其中：aijghks、bijghks、cijs和dijs為SPSAA-s的對偶變量;DSPSAA由S個DSPSAA-s構(gòu)成。

(33)

s.t.

(34)

(35)

aij01ks≤0,1≤i≤I,1≤j≤ui,

1≤k≤K-1;

(36)

aij01Ks≤0,1≤i≤I,1≤j≤ui;

(37)

(38)

aijghks≥0,bijghks≥0,cijs≥0,dijs≥0,1≤i≤I,

1≤j≤ui,1≤g≤I,1≤h≤ug,1≤k≤K。

(39)

其中

(40)

主問題MPSAA如下，T為當(dāng)前迭代次數(shù)。

minw。

(41)

s.t.

約束(4)～(6)，(9)～(14);

(42)

其中

b0=nij·PTik+STgik-(1-zghij)·M。

(43)

Benders分解求解MIPSAA的計算流程總結(jié)如下：

2.2.2 算法加速策略

盡管以上算法可對MIPSAA進行精確求解，但實際計算中算法在每一次迭代中僅添加一個cut效率過低，導(dǎo)致迭代次數(shù)過多，且伴隨著算法的迭代，MPSAA越來越復(fù)雜，每次迭代時間隨之增長。為了加速求解，設(shè)計了3種加速策略。

(1)預(yù)設(shè)optimalitycut策略

Benders分解的目的在于用有限的cut，去逼近原問題，為了加速求解，在算法開始前隨機生成一定數(shù)量cut，并在步驟1首先將預(yù)設(shè)的cut添加進MPSAA中，是對整個Benders分解的一個良好預(yù)熱。

(2)鄰域搜索添加多個cut

1)交換：從搜索起點的第一位開始交換Sc0相鄰數(shù)對。

2)增批：在Sc0的每一位依次插入分批未達(dá)被分至最細(xì)的工件的一個子批，并隨機選擇該工件的另一子批將其批量減去一個最小子批。

3)減批：依次刪除Sc0的每個子批，并將被刪除的子批的批量隨機加至該種工件的一個剩余子批上，當(dāng)Sc0中僅包含一個子批的工件時，不可作刪除處理。

每一次鄰域搜索可向主問題一次性添加若干cut，大大縮小可行域，但是頻繁的搜索會使得約束過多，從而加重主問題求解的負(fù)擔(dān)，故設(shè)計當(dāng)上界在LS次迭代內(nèi)未得到提升，則進行鄰域搜索，引入計數(shù)參數(shù)count_ub記錄上界重復(fù)次數(shù)，每次迭代結(jié)束后更新。

(3)有效不等式(Valid Inequalities, VIs)

(44)

(45)

(46)

(47)

以上均為MPSAA的可行解，卻對應(yīng)同一個Sc，為了避免重復(fù)解，可采用以下兩組VIs。

(1)對于同種工件，優(yōu)先對編號靠前的子批安排批量，約束如下：

yij≥yih,1≤i≤I,1≤i

(48)

(2)同一種工件的不同批量在加工順序的排列中依據(jù)該種工件的批次編號升序排列。

為了表達(dá)出每種工件在整個加工順序中的排序，引入ind變量作為輔助變量。ind無范圍約束，VIs如下：

indij-indgh≥(zghij-1)·M+1,1≤i≤I,

1≤j≤di,1≤g≤I,1≤h≤dg;

(49)

indi,j+1≥indij+1,1≤i≤I,1≤j≤di。

(50)

添加約束(49)和約束(50)至MPSAA中。改進后的Benders分解流程圖如圖3所示。

2.3 啟發(fā)式求解算法

上述精確算法難以在有限時間內(nèi)求解大規(guī)模問題，本文在此基礎(chǔ)上提出基于禁忌搜索的啟發(fā)式算法，將問題分解為以下兩個子問題:

(1)確定每種工件分批數(shù)與分批后所有子批的加工順序。

(2)加工順序確定后每個子批的批量。

子問題(1)通過禁忌搜索算法對加工順序進行優(yōu)化，子問題(2)采用SAA與Benders分解法結(jié)合求解。

2.3.1 禁忌搜索算法

算法對加工順序進行實數(shù)編碼，數(shù)字表示工件種類，數(shù)字在序列中的位置表示順序，數(shù)字在序列中出現(xiàn)的次數(shù)表示工件的子批個數(shù)，例如序列Sd={1 122}表示需加工兩種工件，工件1與工件2均被分為2個子批，工件1的第一個子批首先加工，其次是工件1的第2個子批、工件2的第1個子批、工件2的第2個子批。如圖4所示，每一次迭代從搜索起點Sd0出發(fā)對其進行鄰域搜索并將下一次迭代的搜索起點更新為本次的鄰域最優(yōu)解，鄰域搜索為增批與減批。

增批為選擇編碼中子批個數(shù)小于其子批上限的工件，在Sd0的所有相鄰位中插入該工件的一個子批，增批可使工件的總子批數(shù)加1。減批為依次刪除Sd0中工件編號在編碼中出現(xiàn)次數(shù)大于1的工件的每一個子批，減批可使工件的總子批數(shù)減1。

為避免搜索重復(fù)解，建立禁忌列表(Tabu List, TL)，將Sd0中被破壞的子序列存入TL內(nèi)，舍棄被TL禁忌的鄰域解，若一定迭代次數(shù)(禁忌步長)內(nèi)鄰域最優(yōu)均未出現(xiàn)該子序列，則對其解禁，從TL中刪除，如圖5所示。本文的禁忌搜索算法通過鄰域搜索對加工順序進行優(yōu)化，終止條件為計算時間。其中對加工順序的評價，即Sd的解碼為問題(2)的求解，并通過求解的加工完成時間對該加工順序進行評價，算法終止時輸出最優(yōu)解。

禁忌搜索算法的基本過程如下：

步驟1初始化。迭代次數(shù)計數(shù)器iteration=0。隨機生成Sd0并計算其最大完成時間Ma0。禁忌列表初始化為空，上界UB=Ma0，當(dāng)前最優(yōu)解Sdb=Sd0。

步驟2iteration=iteration+1，對Sd0進行鄰域搜索并舍去被禁忌的解，得到鄰域最優(yōu)Sdlb及其最大完成時間Malb。

步驟3Sd0=Sdlb，UB=min{UB,Malb}。

步驟4更新禁忌列表TL。

步驟5檢查計算時間是否滿足終止條件，若滿足，則迭代中止，輸出Sdb;否則返回步驟2。

2.3.2 解碼：基于SAA與Benders分解的子批量確定

對于每次禁忌搜索確定的工件全部子批的加工順序，其對應(yīng)子批批量無法直接得到，故基于SAA建立隨機規(guī)劃模型MIPTS，并采用Benders分解算法求解。

(1)輸入變量

L為編碼長度，子批編號l∈L;

π為加工順序，π(l)表示第l批工件的工件種類;

ali為輔助變量，ali=1表示第l個子批工件種類為i，否則ali=0。

(2)待優(yōu)化變量

nl為整數(shù)變量，第l個子批包含的最小子批數(shù)，大于0;

tlks為連續(xù)變量，第s個場景下第l個子批在第k臺機器上的加工完成時間。

(51)

s.t.tlks-t(l-1)ks≥nl·minS·PTπ(l)k+STπ(l-1)π(l)k,

2≤l≤L,1≤k≤K,1

(52)

tlks-tl(k-1)s≥nl·minS·PTπ(l)k+STπ(l-1)π(l)k,

1≤l≤L,1≤k≤K,1

(53)

tl0s≥ATπ(l)s,1≤l≤L,1

(54)

(55)

tlks≥0,intnl>0,ali∈(0,1),1≤l≤L,

0≤i≤I,0≤k≤K,1≤s≤S。

(56)

目標(biāo)函數(shù)(51)為最小化S個場景下，最后一個子批在最后一臺機器上的完成時間均值。約束(52)表示同一臺機器上相鄰子批在該機器上的加工完成時間之間的關(guān)系;約束(53)表示同一子批在相鄰兩臺機器上的加工完成時間之間的關(guān)系;約束(54)表示所有工件均在來料后才可開始加工;約束(55)是對各種工件各子批總量的約束;約束(56)是對所有變量取值范圍的約束。

mintLKs。

(57)

s.t.

約束(54);

STπ(l-1)π(l)k,2≤l≤L,1≤k≤K;

(58)

STπ(l-1)π(l)k,1≤l≤L,1≤k≤K;

(59)

tl0s≥ATπ(l)s,1≤l≤L。

(60)

SPTS-s的對偶問題DSPTS-s如下，其中λlks，γlks和δls分別為MIPTS-s的對偶變量。

(61)

s.t.

-γl1s+δls≤0,1≤l≤L;

(62)

λlk-λ(l+1)ks+γlks-γl(k+1)s≤0,

2≤l

(63)

λlKs-λ(l+1)Ks+γlKs≤0,2≤l

(64)

-λ2ks+γ1ks-γ1(k+1)s≤0,1≤k

(65)

-λ2Ks+γ1Ks≤0;

(66)

λLks+γLks-γL(k+1)s≤0,1≤k

(67)

λLKs+γLKs≤1;

(68)

λlks≥0,γlks≥0,δls≥0,1≤l≤L,1≤k≤K;

(69)

(70)

minw。

(71)

s.t.

約束(55);

(72)

nl>0,1≤l≤L。

(73)

其中:約束(72)由所有的optimalitycut組成;約束(73)為變量的取值約束。

計算流程可總結(jié)如下：

步驟2求解DSPTS，通過DSPSAA的目標(biāo)值均值更新UB并構(gòu)造cut更新MPTS。

3 數(shù)值實驗

本文算例規(guī)模用I-N-K表示，即需要在K臺機器上加工I種工件，每種工件最多可分為N個子批。同一算例中不同種工件的來料時間服從不同的隨機分布，分布表示如下：正態(tài)分布-N,指數(shù)分布-E,三角分布-T。最小批量minS=1。本章中，PTik采用來自TAILLARD[26]與MARIMUTHU等[27]的兩組benchmark數(shù)據(jù)，分別記為數(shù)據(jù)集Tail和數(shù)據(jù)集Mari，同時根據(jù)實際制造中的調(diào)研信息，本文將不同種工件在不同機器上的換模時間設(shè)定為STijk=0.5×(PTik+PTjk)。全部算法采用C++編寫算法程序，運行于3.10 GHz CPU和512 GB RAM的計算機，求解數(shù)學(xué)規(guī)劃采用CPLEX求解軟件。

3.1 隨機性與靈敏度分析

場景數(shù)S過小會影響隨機問題的求解精度，過大會使模型的求解難度大大提升，故在進行數(shù)值實驗前需對S的選取進行靈敏度分析，結(jié)果如表1所示。其中算例1～算例7規(guī)模相同，工件1與工件2的來料時間服從指數(shù)分布，工件3的來料時間服從三角分布，PTik來自于算例Tail，所有算例的輸入?yún)?shù)均相同，僅對場景數(shù)S做出區(qū)分。

表1 場景數(shù)的靈敏度分析

由表1可知，S的選取對于求解精度與計算時間均有較為明顯的影響。由于用Benders分解算法求解MIPSAA的子問題時需求解S組對偶問題，隨著場景數(shù)的增多，Benders分解的計算時間也隨之增加，當(dāng)S過大時，求解子問題會耗費大量時間，而調(diào)度方案在S大于75后已不再改變，且結(jié)果也基本穩(wěn)定，故S過大不僅對求解精度提升不高，還會影響問題的求解速度。當(dāng)場景數(shù)達(dá)到100時,需近4 h才能夠完成求解，但是S過小時,由于場景數(shù)過少,則對于隨機性的處理不夠，解的穩(wěn)定性不高，結(jié)果的波動較大，從靈敏度的分析結(jié)果來看，當(dāng)場景數(shù)達(dá)到75時，算例對最優(yōu)加工順序與分批情況的求解結(jié)果已基本趨向收斂，可以確定S取75以上可滿足實驗需求。

為了驗證隨機性對求解結(jié)果的影響，通過實驗對工件來料時間的隨機性進行分析，其中算例8～12規(guī)模相同，PTik來自于算例Tail，場景數(shù)為75，工件1的來料時間服從指數(shù)分布，期望為200，工件2與工件3均服從正態(tài)分布，期望分別為200與150，算例8～13中3種工件的來料時間期望、加工時間、換模時間均相同，工件2與工件3的來料時間方差不同，各算例詳細(xì)信息如表2所示。由表2可知，由于不同算例下工件來料時間期望相同，求解結(jié)果在一定范圍內(nèi)波動，但隨著方差的增大，結(jié)果的波動逐漸增大，且各算例的最優(yōu)加工順序與最優(yōu)分批情況有明顯不同，來料時間的隨機性對于工件的最優(yōu)調(diào)度方案的求解有明顯影響。

表2 來料時間隨機性影響分析

3.2 精確算法求解實驗

為驗證本文設(shè)計的SAA&BD算法的準(zhǔn)確性與求解效率，進一步采用benchmark數(shù)據(jù)進行數(shù)值實驗，并將結(jié)果與直接用CPLEX求解MIPSAA的解對比來進行SAA&BD算法效能的評估。以規(guī)模為3-3-5、3-3-8和3-3-10的20組數(shù)據(jù)作為算例13～32，其中算例13～22的加工時間數(shù)據(jù)來自于數(shù)據(jù)集Tail，算例23～32的加工時間數(shù)據(jù)來自于數(shù)據(jù)集Mari，SAA&BD算法預(yù)設(shè)cut數(shù)P=500，鄰域搜索觸發(fā)參數(shù)LS=20，場景數(shù)S=75,且不同算例下工件來料時間所服從的分布均不同，同一分布下不同算例的工件來料時間均值及方差不等，設(shè)置CPLEX計算時間為6 h，如表3所示為算例13～32的計算結(jié)果。其中GapC=(UBC-LBC)/UBC·100%,Gap=(UBC-UBB)/UBB·100%。

表3 小規(guī)模算例benchmark數(shù)值實驗結(jié)果

續(xù)表3

由表3可知，在場景數(shù)等于75時,對于不同benchmark均能在3.5 h內(nèi)完成求解，算法求解時間較為穩(wěn)定。同時算例13～32的工件來料時間均服從不同的隨機分布，Benders分解的求解時間并沒有出現(xiàn)明顯的差異性，穩(wěn)定性較高，可以認(rèn)為該算法足以應(yīng)對復(fù)雜的隨機環(huán)境。與CPLEX求解MIPSAA對比，即使是小規(guī)模算例，CPLEX也無法在6 h內(nèi)達(dá)到收斂，且上下界差距較大，而采用Benders分解，雖與CPLEX的結(jié)果相近，但Benders分解在求解時間上優(yōu)于CPLEX，且由于上下界相等，可以證明其得到的是最優(yōu)解。進一步地，在各隨機場景下Benders分解均能求得最優(yōu)解且求解時間穩(wěn)定，而CPLEX在6 h內(nèi)求得的解上下界Gapc在5%～60%之間波動，穩(wěn)定性較差，可見在隨機環(huán)境中，Benders分解有較高的穩(wěn)定性。

3.3 啟發(fā)式算法求解問題實驗

對于規(guī)模較大的問題，本文同樣采用兩組benchmark數(shù)據(jù)分別進行數(shù)值實驗，為了驗證算法的有效性,本文以收斂性為評價指標(biāo),將TS與遺傳算法、CPLEX直接求解方法進行比較和分析。采用算例規(guī)模為5-5-10的20組數(shù)據(jù)作為本文算例33～52，其中算例33～42的加工時間數(shù)據(jù)來自于數(shù)據(jù)集Tail，算例43～52的加工時間數(shù)據(jù)來自于數(shù)據(jù)集Mari。由于工件種類較多，工件的隨機來料情況更為復(fù)雜，取S=200。設(shè)置CPLEX計算時間為4 h，禁忌搜索與遺傳算法終止條件為迭代1 h，如表4所示為算例33～52的計算結(jié)果。RTS、RGA與RCPLEX分別為TS、遺傳算法(GA)以及CPLEX求得的最大完成時間，其中GAPGA與GAPCPLEX計算方式分別為(RGA-RTS)/RTS·100%、(RCPLEX-RTS)/RTS·100%。

表4 大規(guī)模算例benchmark數(shù)值實驗結(jié)果

續(xù)表4

遺傳算法的基本設(shè)計如下：首先隨機生成50個個體組成初始種群，個體編碼方式與TS一致，通過對編碼解碼求得每條染色體的適應(yīng)度值。在每一輪迭代中，基于適應(yīng)度采用輪盤賭方法對群體進行選擇，對被選擇的兩個個體進行交叉，交叉方式為隨機生成兩個不同的基因位進行單點交換，交叉完成后對個體進行變異操作，隨機選取5個基因位作變動。交叉概率為0.8，變異概率為0.1，一輪迭代后得到下一個種群，當(dāng)計算時間達(dá)到1 h，輸出進化過程中所得到的具有最優(yōu)適應(yīng)度個體作為最優(yōu)解輸出，終止計算。

由表4可知，與GA對比，同樣的計算時間內(nèi)，TS求得的最優(yōu)解均明顯優(yōu)于GA，且GAPGA在4%～11%。如圖6所示為算例33與53的收斂圖(橫坐標(biāo)為迭代次數(shù)，縱坐標(biāo)為當(dāng)前最優(yōu)解)，TS求解1 h迭代次數(shù)約為100次。由圖6可明顯看出，TS在100次迭代內(nèi)迅速收斂且已基本穩(wěn)定，而GA在500次迭代內(nèi)并不能快速找到一個優(yōu)秀的解并達(dá)到收斂。

與CPLEX求解混合整數(shù)規(guī)劃模型相比，TS在1 h內(nèi)得到的最優(yōu)解均優(yōu)于GA與CPLEX在4 h內(nèi)給出的最優(yōu)上界，結(jié)果的GAP值可達(dá)10%以上，且在10%～60%之間波動較大。對于部分算例，如算例27中CPLEX上界與TS最優(yōu)解的GAP大于57%，算例44與49中CPLEX在4 h內(nèi)均未找到可行解，可見禁忌搜索算法無論求解時間還是求解效果都優(yōu)于CPLEX。從TS自身的求解性能來看，算例33～42與算例43～52分別來自于兩組不同的數(shù)據(jù)集，但均能夠在100次迭代次數(shù)內(nèi)達(dá)到收斂，具有一定的穩(wěn)定性，且對于工件來料時間分布并不敏感，對于來料時間的隨機性不再限定分布必須服從特定的分布，增加了本文研究所針對問題的適應(yīng)性。

4 結(jié)束語

本文以工件來料時間隨機且考慮機器換模時間的柔性流水車間的批量流調(diào)度問題為研究對象，基于SAA方法建立了混合整數(shù)規(guī)劃模型，并設(shè)計了精確求解與啟發(fā)式求解兩套求解算法。針對小規(guī)模問題，采用Benders分解加速模型的求解，實現(xiàn)了小規(guī)模問題的精確求解。針對難以在有限時間內(nèi)實現(xiàn)精確求解的大規(guī)模問題，設(shè)計了禁忌搜索算法。最終，本文通過數(shù)值實驗對比所設(shè)計的算法與遺傳算法以及CPLEX的求解結(jié)果，驗證了模型的可行性，并證明了Benders分解在求解復(fù)雜混合整數(shù)規(guī)劃模型時的有效性與禁忌搜索算法在求解大規(guī)模問題上的優(yōu)勢。下一步，將針對更為復(fù)雜的混合車間批量流隨機調(diào)度問題進行研究。