考慮復(fù)雜隨機(jī)來(lái)料的批量流作業(yè)調(diào)度

章 旸，王鋮愷，劉 冉

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

0 引言

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

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

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

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

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

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

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

1.1 問(wèn)題描述

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

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

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

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

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

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

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

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

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

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

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

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

yij為0-1變量，nij>0時(shí)，yij=1;nij=0時(shí)，yij=0。

(2)隨機(jī)變量

ATis為場(chǎng)景s中工件i的來(lái)料時(shí)間。

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

2 求解算法

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

dTy+(b-By)Tup≤w。

(24)

式(24)稱為optimalitycut。

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

(b-By)Tup≤0。

(25)

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

(26)

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

allthefeasibilitycuts}。

(27)

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

2.1 基于Benders分解的精確算法

2.1.1 主問(wèn)題與子問(wèn)題構(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的對(duì)偶問(wèn)題DSPSAA-s如下，其中：aijghks、bijghks、cijs和dijs為SPSAA-s的對(duì)偶變量;DSPSAA由S個(gè)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)

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

minw。

(41)

s.t.

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

(42)

其中

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

(43)

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

2.2.2 算法加速策略

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

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

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

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

1)交換：從搜索起點(diǎn)的第一位開(kāi)始交換Sc0相鄰數(shù)對(duì)。

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

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

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

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

(44)

(45)

(46)

(47)

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

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

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

(48)

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

為了表達(dá)出每種工件在整個(gè)加工順序中的排序，引入ind變量作為輔助變量。ind無(wú)范圍約束，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中。改進(jìn)后的Benders分解流程圖如圖3所示。

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

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

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

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

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

2.3.1 禁忌搜索算法

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

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

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

禁忌搜索算法的基本過(guò)程如下：

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

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

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

步驟4更新禁忌列表TL。

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

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

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

(1)輸入變量

L為編碼長(zhǎng)度，子批編號(hào)l∈L;

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

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

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

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

tlks為連續(xù)變量，第s個(gè)場(chǎng)景下第l個(gè)子批在第k臺(tái)機(jī)器上的加工完成時(shí)間。

(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個(gè)場(chǎng)景下，最后一個(gè)子批在最后一臺(tái)機(jī)器上的完成時(shí)間均值。約束(52)表示同一臺(tái)機(jī)器上相鄰子批在該機(jī)器上的加工完成時(shí)間之間的關(guān)系;約束(53)表示同一子批在相鄰兩臺(tái)機(jī)器上的加工完成時(shí)間之間的關(guān)系;約束(54)表示所有工件均在來(lái)料后才可開(kāi)始加工;約束(55)是對(duì)各種工件各子批總量的約束;約束(56)是對(duì)所有變量取值范圍的約束。

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的對(duì)偶問(wèn)題DSPTS-s如下，其中λlks，γlks和δls分別為MIPTS-s的對(duì)偶變量。

(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)為變量的取值約束。

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

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

3 數(shù)值實(shí)驗(yàn)

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

3.1 隨機(jī)性與靈敏度分析

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

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

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

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

表2 來(lái)料時(shí)間隨機(jī)性影響分析

3.2 精確算法求解實(shí)驗(yàn)

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

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

續(xù)表3

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

3.3 啟發(fā)式算法求解問(wèn)題實(shí)驗(yàn)

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

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

續(xù)表4

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

由表4可知，與GA對(duì)比，同樣的計(jì)算時(shí)間內(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)并不能快速找到一個(gè)優(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%之間波動(dòng)較大。對(duì)于部分算例，如算例27中CPLEX上界與TS最優(yōu)解的GAP大于57%，算例44與49中CPLEX在4 h內(nèi)均未找到可行解，可見(jiàn)禁忌搜索算法無(wú)論求解時(shí)間還是求解效果都優(yōu)于CPLEX。從TS自身的求解性能來(lái)看，算例33～42與算例43～52分別來(lái)自于兩組不同的數(shù)據(jù)集，但均能夠在100次迭代次數(shù)內(nèi)達(dá)到收斂，具有一定的穩(wěn)定性，且對(duì)于工件來(lái)料時(shí)間分布并不敏感，對(duì)于來(lái)料時(shí)間的隨機(jī)性不再限定分布必須服從特定的分布，增加了本文研究所針對(duì)問(wèn)題的適應(yīng)性。

4 結(jié)束語(yǔ)

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