孟磊磊 張超勇 詹欣隆 洪 輝 羅 敏
1.華中科技大學(xué)數(shù)字制造裝備與技術(shù)國家重點實驗室,武漢,4300742.湖北汽車工業(yè)學(xué)院電氣與信息工程學(xué)院,十堰,442002
據(jù)統(tǒng)計,我國擁有機床800多萬臺,每臺機床平均功率按10 kW算,總功率超過8 000萬kW,耗電量是總裝容量為2 250萬kW的三峽電站的3倍多[1]。然而,機床在實際生產(chǎn)中,大部分時間處于空載狀態(tài),據(jù)統(tǒng)計,機床大約80%的能量消耗在待機以及關(guān)機重啟過程中,真正用于加工的能量占比不足20%[2],GUTOWSKI等[3]給出的一條汽車自動生產(chǎn)線的能量效率只有14.8%。機床處于空閑等待狀態(tài)所消耗的能量為無效能耗,節(jié)能調(diào)度的主要目的就是減少這部分能耗。
并行機調(diào)度問題(parallel machine scheduling problem,PMSP)在實際生產(chǎn)制造過程中比較常見。不相關(guān)并行機調(diào)度問題(unrelated parallel machine scheduling problem, UPMSP) 是PMSP中最普遍的一類問題, 工件的加工時間取決于所分配的機床,UPMSP已經(jīng)被證明為NP-hard問題。目前,對UPMSP的研究主要集中在最小化最大完工時間、成本、交貨率等方面[4-6],而對車間能耗方面不太重視。
當機床連續(xù)待機時間較長時,對于允許中途關(guān)機的機床可以通過關(guān)機/重啟策略來減小機床待機能耗。關(guān)機/重啟策略最早由MOUZON等[2]提出,通過對非瓶頸機床實行關(guān)機/重啟策略,避免機床長時間處于空閑狀態(tài),可以節(jié)約80%的待機能耗。隨后關(guān)機/重啟策略逐步被應(yīng)用于單機[7]、并行機[8-11]調(diào)度研究。
針對UPMSP,文獻[12]以加權(quán)拖期成本和能耗成本之和為目標,提出了求解該問題的多種啟發(fā)式規(guī)則方法。文獻[13]研究了并行機批量調(diào)度問題,以最小化能耗和最大完工時間為目標,提出了一種基于Pareto的多目標蟻群優(yōu)化算法。針對加工時間可控的UPMSP,文獻[14]以最小化能耗和最大完工時間為目標,提出了求解該問題的多目標協(xié)同果蠅優(yōu)化算法;文獻[15]同時以最小化加工能耗和最大完工時間為目標,提出了一種高效的多目標memetic差分進化算法。
UPMSP可以通過精確方法,如混合整數(shù)規(guī)劃(mixed integer programming ,MIP)模型、分支定界等[16]方法求解,也可通過元啟發(fā)式方法(如遺傳算法[17]、蟻群算法[8]、分布估計算法[5]等)和啟發(fā)式規(guī)則[12]等方法求解。精確算法求解速度慢,對于中大規(guī)模問題甚至很難求出問題的可行解,但是針對小規(guī)模問題,精確算法可以求出調(diào)度問題的最優(yōu)方案,而啟發(fā)式方法和元啟發(fā)式方法為近似求解方法,雖然在較短時間內(nèi)可以得出問題的可行解,但解的質(zhì)量不能保證為最優(yōu)解。MIP模型是認識調(diào)度問題的基礎(chǔ),是探索與挖掘調(diào)度分配規(guī)則的關(guān)鍵所在[18]。高效的MIP模型可以在相同的時間內(nèi)得到更好的計算結(jié)果,因此構(gòu)建高效MIP模型具有重要意義。
本文結(jié)合UPMSP自身的特點,以節(jié)能為目標,考慮關(guān)機/重啟策略,提出5個解決該問題的MIP模型,并從建模過程、線性化過程、模型尺寸復(fù)雜度以及計算復(fù)雜度等方面對5個MIP模型進行對比分析。
參數(shù)定義見表1。其中,為保證每臺機床上都安排工件加工,每臺機床最多可安排n-m+1個工件,每臺機床最大位置數(shù)為n-m+1。
UPMSP可描述為:n個工件在m臺機床上加工,每個工件只有一道工序,且同一工件在不同機床上的加工時間是不同的。該問題滿足以下基本假設(shè): 不同工件的釋放時間、交貨期是不同的;所有工件必須在規(guī)定的交貨期內(nèi)完成加工;機床在0時刻處于可用狀態(tài);所有工件在所有機床上的加工時間、加工功率等是已知的;同一機床上不同工件間的轉(zhuǎn)換時間忽略不計;任一工件可在任一機床上加工;工件一旦開始加工則不可中斷;對于每個機床,在同一時刻最多只能加工一個工件;對于每個工件,在同一時刻最多只能在一臺機床上加工。
表1 參數(shù)定義
UPMSP的目的就是合理地將工件安排在機床上,且確定工件在機床上的加工順序以達到某些調(diào)度目標,如最大完工時間、成本、能耗等的最優(yōu)。本文研究UPMSP節(jié)能調(diào)度問題,目標為最小化車間總能耗。
車間能耗主要包括機床能耗、公共能耗兩部分,機床能耗進一步可分為加工能耗及待機能耗。
1.3.1加工能耗
(1)
總加工能耗
(2)
因為每個工件只能選擇在一臺機床上加工,故引入0-1變量Yi,k,如果工件i選擇在機床k上加工,則Yi,k=1,否則,Yi,k=0。
1.3.2待機能耗
待機能耗是指由于上一工件已經(jīng)加工完,下一個工件還沒有到達,從而使機床出現(xiàn)閑置狀態(tài)所消耗的能量,機床k的第t到t+1位置之間的待機能耗
(3)
式中,F(xiàn)k,t為機床k上第t個位置的結(jié)束時間;Sk,t為機床k上第t個位置的開始時間。
機床k的總待機能耗
(4)
所有機床總待機能耗
(5)
1.3.3公共能耗
公共能耗是指車間公共設(shè)施的能源消耗,是指為了維持車間正常運行所必須消耗的能源,主要包括照明、通風(fēng)、供暖、空調(diào)等耗能,其值為公共功率P0與最大完工時間Cmax的乘積,可表示為
EC=P0Cmax
(6)
1.3.4車間總能耗
車間總能耗ETC為總加工能耗、待機能耗以及公共能耗之和:
(7)
當機床的一個連續(xù)待機時間段Sk,t+1-Fk,t比較長時,可以實施關(guān)機/重啟策略,節(jié)約無用待機能耗,實現(xiàn)關(guān)機/重啟策略的最短待機時間,即空載平衡時間
(8)
空載平衡時間Tk,B表示機床在空載的過程中,當空載的時間不小于一次關(guān)機/重啟策略所需要的時間Tk,且機床空載所消耗的能耗不小于機床一次關(guān)機/重啟所需要的能耗Ek,min時才可實施關(guān)機/重啟策略。
當引入關(guān)機/重啟策略后,待機段能耗Ek,t可表示為
(9)
當引入關(guān)機/重啟策略后,車間總能耗
(10)
其中,Zk,t為0-1變量,用來控制是否實行關(guān)機/重啟策略,當Zk,t=1時,表示在機床k的第t到t+1位置之間實行關(guān)機/重啟策略,此時Ek,t為一次關(guān)機/重啟所需要的能耗Ek,min;否則,Zk,t=0,Ek,t根據(jù)具體待機時間而定。
目前關(guān)于UPMSP的MIP模型較多, 但是大部分是針對時間目標的[16-17,19]。本文針對UPMSP自身的特點,以節(jié)能為目標,同時考慮關(guān)機/重啟策略,基于機床位置的Wagner建模思想[20],建立5個考慮關(guān)機/重啟策略的MIP模型,其中2個為基于空閑時間的非線性模型,3個為線性模型(2個線性模型通過將2個非線性模型引入中間決策變量進行線性化處理得到,第3個線性模型基于空閑能耗的建模思想直接構(gòu)建)。
基于機床位置的Wagner建模方法建立在“機床位置”概念上,即將一個機床按照時間先后分成若干個段,每一段即成為一個位置,一個位置最多只能安排一個工件,因此,只要確定了工件與機床位置的對應(yīng)關(guān)系,工件在機床上的排序即可得到。
5個MIP模型按照建模思路可進一步細分為兩類:第一類是基于空閑時間的建模方法,包括非線性模型1、線性模型1與非線性模型2、線性模型2;第二類是基于空閑能耗的建模方法,包括線性模型3?;诳臻e時間的建模方法是指機床待機能耗通過待機段時間與待機功率來計算,而基于空閑能耗的建模方法直接定義空閑段能耗決策變量。
因為基于Wagner建模方法的MIP模型必須引入機床位置占用變量Xi,k,t,Xi,k,t與機床選擇變量Yi,k存在對應(yīng)關(guān)系,因此本文5個模型將不再使用Yi,k,以減少0-1決策變量數(shù),降低模型復(fù)雜度。對應(yīng)關(guān)系如下:
(11)
模型1包括非線性模型1以及線性模型1。其中線性模型1通過引入中間決策變量將非線性模型1的目標函數(shù)進行線性化處理之后得到。
非線性模型1有Xi,k,t、Zk,t、Sk,t及Cmax4個決策變量,其中,Xi,k,t用于確定工件的機床位置選擇,Zk,t用于確定在機床相應(yīng)位置是否實行關(guān)機/重啟操作,Sk,t和Cmax為連續(xù)決策變量。Xi,k,t和Zk,t的表達式分別如下:
為了將非線性目標函數(shù)進行線性化處理,引入中間決策變量Ai,k,t、Uk,t和Wk,t,其中Uk,t、Wk,t為連續(xù)決策變量,Ai,k,t為0-1中間決策變量,因此線性模型1共有Xi,k,t、Zk,t、Sk,t、Cmax、Ai,k,t、Uk,t和Wk,t7個決策變量。
非線性模型1的目標函數(shù)為
(12)
式(12)中,等號右邊第一項表示機床總空閑能耗(待機或關(guān)機/重啟),第二項為機床總加工能耗,第三項為公共能耗。由式(12)可以看出,目標函數(shù)是非線性的,存在決策變量相乘的情況,如(1-Zk,t)Sk,t+1、(1-Zk,t)Sk,t以及(1-Zk,t)Xi,k,t,非線性模型求解非常復(fù)雜。本文引入中間決策變量Uk,t、Wk,t、Ai,k,t,用Uk,t+1代替(1-Zk,t)Sk,t+1、Wk,t代替(1-Zk,t)Sk,t、Ai,k,t代替(1-Zk,t)Xi,k,t,通過添加相應(yīng)的約束,實現(xiàn)非線性模型到線性模型的轉(zhuǎn)換。
線性模型1的目標函數(shù)為
(13)
約束條件為
(14)
(15)
(16)
ri≤Sk,t+M(1-Xi,k,t) ?t∈T
(17)
(18)
Sk,t≥0 ?t∈T
(19)
(20)
(21)
(22)
(23)
Uk,t+1≥Sk,t+1-MZk,t?t∈T1
(24)
Uk,t+1≤Sk,t+1+MZk,t?t∈T1
(25)
Uk,t+1≤M(1-Zk,t) ?t∈T1
(26)
Uk,t≥0 ?t∈T
(27)
Wk,t≥Sk,t-MZk,t?t∈T1
(28)
Wk,t≤Sk,t+MZk,t?t∈T1
(29)
Wk,t≤M(1-Zk,t) ?t∈T1
(30)
Wk,t≥0 ?t∈T1
(31)
Ai,k,t≤Xi,k,t?t∈T1
(32)
Ai,k,t≤1-Zk,t?t∈T1
(33)
Ai,k,t≥Xi,k,t-Zk,t?t∈T1
(34)
式(14)~式(34)中,?i∈I,?k∈K。式(14)表示任一工件只能選擇在一個機床上加工,且只占用該機床的一個位置;式(15)表示任一機床的任一位置最多只能安排一個工件加工;式(16)表示任一機床的位置按照先后順序安排工件加工;式(17)表示工件只能在到達之后才可以開始加工;式(18)表示工件必須在交貨期內(nèi)完成加工;式(19)約束決策變量范圍;式(20)表示最大完工時間約束;式(21)、式(22)表示關(guān)機/重啟策略約束,約束當機床k的t到t+1位置間存在關(guān)機/重啟策略時,即Zk,t=1, 第t+1位置的開始時間與第t位置的開始時間加上該位置的加工時間之差必不小于機床k的空載平衡時間Tk,B,否則,不存在關(guān)機/重啟策略;由于平時加工過程中,機床是不允許頻繁關(guān)機重啟的(頻繁關(guān)機重啟對機床電器元器件壽命影響很大),因此引入式(23)來限制關(guān)機重啟次數(shù); 式(24)~式(27)約束Uk,t+1=(1-Zk,t)Sk,t+1恒成立,其中,成對約束(式(24)與式(25))保證了當機床k的t到t+1位置間不存在關(guān)機/重啟策略時,即Zk,t=0時,保證Uk,t+1=(1-Zk,t)Sk,t+1成立;式(26)、式(27)保證了Zk,t=1時,Uk,t+1=(1-Zk,t)Sk,t+1成立,同理,式(28)~式(31)保證了Wk,t=(1-Zk,t)Sk,t恒成立,式(32)~式(34)保證了Ai,k,t=(1-Zk,t)Xi,k,t恒成立。
模型2包括非線性模型2和線性模型2。與模型1相比,模型2引入連續(xù)決策變量Fk,t,從而減少目標函數(shù)中的非線性項、中間決策變量以及相應(yīng)約束的數(shù)量。
非線性模型2的目標函數(shù)中含有非線性項(1-Zk,t)(Sk,t+1-Fk,t),通過引入中間決策變量線性化后即為線性化模型2的目標函數(shù)。線性化過程與模型1相似,引入中間決策變量Uk,t、Vk,t,分別代替非線性項(1-Zk,t)Sk,t+1、(1-Zk,t)Fk,t,實現(xiàn)非線性目標函數(shù)到線性目標函數(shù)的轉(zhuǎn)換。
模型2中的決策變量為:Xi,k,t、Zk,t、Sk,t、Cmax、Uk,t(含義與模型1中的相同);機床k上第t個位置的結(jié)束時間Fk,t;用來代替非線性項(1-Zk,t)Fk,t的中間決策變量Vk,t。
非線性模型2的目標函數(shù)為
(35)
線性模型2的目標函數(shù)為
(36)
約束條件為
(37)
Fk,t-M(1-Xi,k,t)≤di?t∈T
(38)
Fk,n-m+1≤Cmax
(39)
Sk,t+1-Fk,t≤Tk,B+MZk,t?t∈T1
(40)
Sk,t+1-Fk,t≥Tk,BZk,t?t∈T1
(41)
Vk,t≥Sk,t-MZk,t?t∈T1
(42)
Vk,t≤Sk,t+MZk,t?t∈T1
(43)
Vk,t≤M(1-Zk,t) ?t∈T1
(44)
Vk,t≥0 ?t∈T1
(45)
其中,?i∈I,?k∈K。式(37)表示機床位置結(jié)束時間等于其開始時間加上安排在該位置工件的加工時間;式(38)~式(41)的含義分別與式(18)、式(20)~式(22)約束含義相同;式(42)、式(45)用來約束Vk,t=(1-Zk,t)Fk,t恒成立,其中式(42)~式(43)用來約束當Zk,t=0時,Vk,t=Fk,t=(1-Zk,t)Fk,t成立,式(44)、式(45)用來保證當Zk,t=1時,Vk,t=0=(1-Zk,t)Fk,t成立。
模型3的目標函數(shù)為
(46)
約束條件為
(47)
(48)
?k∈K,t∈T1
模型對比從尺寸復(fù)雜度與計算復(fù)雜度進行,尺寸復(fù)雜度主要包括0-1決策變量數(shù)、約束數(shù)以及連續(xù)決策變量數(shù)3個方面。計算復(fù)雜度通過規(guī)定時間內(nèi)可求最優(yōu)解總數(shù)Q來判定,包括目標函數(shù)值的容差Gap=0最優(yōu)解個數(shù)Q0與Gap≠0最優(yōu)解個數(shù)Q1。當Q相同時,對比Q0;當Q與Q0都相同時,對比Q1;Q、Q0與Q1越大,模型越好。當Q、Q0與Q1都相同時,對比求解時間ttotal,ttotal也是一個重要評價指標,越小越好。Gap可定義為|SC-SB|/|SC|,其中SC表示目前為止可以找到的最優(yōu)解,SB表示可能的最優(yōu)解,是當前所有解的下限??梢?,Gap值越小越好,當Gap=0時,則得到問題的最優(yōu)解,程序會自動停止。由此,Gap值也常作為評價MIP模型求解效果的一個指標。
本文所有混合整數(shù)線性模型都由商業(yè)軟件IBM ILOG CPLEX12.7.1求解,編程語言采用CPLEX Studio自帶OPL語言編寫。CPLEX可以求解線性模型以及非線性模型。模型運行最長時間設(shè)置為3 600 s,所有實例獨立運行3次,最終結(jié)果取平均值。所有實例在Dell Vostro 3900臺式機上運行,該電腦配置為:Win7系統(tǒng)、i5-4460 3.20GHz四核處理器、8G內(nèi)存。
如果模型在3 600 s之內(nèi)自行停止,則可得到最優(yōu)解且可證明所得到的解為最優(yōu)解,即Gap=0,如果到截止時間3 600 s,程序強制停止,此時有可能得到Gap≠0最優(yōu)解,但是Gap≠0,是因為雖然求出了最優(yōu)解,但是在規(guī)定時間內(nèi)不能證明該解為最優(yōu)的。本文總共測試10組規(guī)模不同的實例。因為目前沒有關(guān)于UPMSP能耗相關(guān)的對比實例,因此本文隨機生成10組規(guī)模不同的實例。
表2 每個約束方程的約束數(shù)目
表2給出了上文所列每個約束方程的約束數(shù)量,結(jié)合每個模型所含約束方程,進而可以得每個模型的約束數(shù)(表3)。由各個MIP模型的建模過程可得出每個模型的尺寸復(fù)雜度,見表3。模型針對具體實例的尺寸復(fù)雜度見表4。其中10*2表示工件數(shù)*機床數(shù),其他類推。
表3 所有模型尺寸復(fù)雜度
表4 所有模型針對10組實例的尺寸復(fù)雜度
通過表3和表4可知,在0-1決策變量數(shù)方面,非線性模型1最多,其他模型相同。這是因為其他模型只含有Xi,k,t、Zk,t兩個0-1決策變量,而模型1除含有Xi,k,t、Zk,t兩個0-1決策變量外,還含有Ai,k,t。
在約束數(shù)方面,由表3和表4 可以看出,從多到少依次為線性模型1、線性模型2、非線性模型2、非線性模型1、模型3。這是由于線性模型3基于待機能耗的建模思想,不需要線性化處理目標函數(shù)以及相關(guān)的約束方程,約束最少。線性模型1需要引入3個中間決策變量(Ai,k,t、Uk,t以及Wk,t)以及相關(guān)的約束,而線性模型2只需要引入2個中間決策變量(Uk,t與Vk,t)及相關(guān)約束,因此,線性模型1約束數(shù)多于線性模型2。非線性模型1與非線性模型2因為沒有引入中間決策變量及相關(guān)約束,故約束數(shù)分別少于線性模型1與線性模型2。
在連續(xù)決策變量數(shù)方面,由表3和表4可以看出,連續(xù)決策變量數(shù)從多到少依次為線性模型2、線性模型1、非線性模型2、模型3、非線性模型1。
決策變量數(shù)以及約束數(shù)嚴重影響MIP模型的求解效率。對MIP模型影響程度從大到小依次為0-1決策變量數(shù)、約束數(shù)以及連續(xù)決策變量數(shù)[21]。除此之外,非線性特征對MIP模型的求解效率影響也很大。
表5 所有模型針對10組實例的求解結(jié)果
5種模型針對10組實例的求解結(jié)果見表5。其中,SO表示最優(yōu)解(可通過模型3長時間求解得到),SC表示當前解(如果在3 600 s內(nèi)程序自動停止,則為最優(yōu)解,否則對應(yīng)3 600 s所能得到的最好解),帶“*”的解表示此解為可行解但非最優(yōu)解。
由表5可知,在10個實例中,在規(guī)定時間3 600 s內(nèi),非線性模型1及線性模型1求解效果最差,只可求出10個實例中的5個,包括3個Gap=0最優(yōu)解以及2個Gap≠0最優(yōu)解;線性模型1由于引入3個中間決策變量(Ai,k,t、Uk,t、Wk,t)及相關(guān)的約束,決策變量數(shù)以及約束數(shù)遠多于非線性模型1,這嚴重影響了線性模型1的求解效率,因此求解時間長于非線性模型1。
非線性模型2可以求得10個實例中的7個最優(yōu)解,其中6個Gap=0 最優(yōu)解以及1個Gap≠0最優(yōu)解,效果差于線性模型2(其可求得8個Gap=0 最優(yōu)解),且求解時間長于線性模型2,這是由非線性模型2的非線性特征影響的。
非線性模型2和線性模型2由于添加了機床位置結(jié)束時間決策變量Fk,t,非線性模型2的非線性項少于非線性模型1,線性模型2的0-1決策變量以及約束數(shù)遠少于線性模型1,因此求解效果好于非線性模型1和線性模型1。
模型3基于空閑能耗的建模思想,目標函數(shù)是線性的,不需要對其進行線性化處理以及相應(yīng)的中間決策變量、約束等,從而使模型3具有最少的0-1決策變量、連續(xù)決策變量以及較少的約束數(shù),從而模型3求解效果最好。
除此之外,模型的求解效果還依賴具體實例的加工數(shù)據(jù),如針對實例20*5,非線性模型2求解時間為224.25 s,線性模型2求解時間為320.00 s,線性模型3求解時間為451.23 s,非線性模型2的求解效果好于線性模型2以及線性模型3,線性模型2效果好于線性模型3。
實例25*3的具體數(shù)據(jù)如表6和表7所示。圖1為實例25*3的最優(yōu)甘特圖。對圖1進行分析,由于考慮了不同工件釋放時間以及交貨期,因此每個工件必須在自身釋放時間-交貨期時間窗內(nèi)進行加工。工件1~8被分配到機床2上進行加工,機床2在0時刻開機,在212時刻關(guān)機。工件10、12、17、18、21、22、25被分配到機床1上進行加工,機床1在230時刻才開機而不是在0時刻開機,從而可以減少待機能耗。工件12在291時刻加工完后,此時下一工件17還沒有到達(其釋放時間為320),此時機床將處于長時間待機狀態(tài)。工件17在363時刻開始加工,在388時刻完成加工,完工時間小于交貨期400。時刻291~363長達72的待機時間將導(dǎo)致機床浪費大量待機能耗,在該待機段內(nèi)對機床1實行關(guān)機/重啟節(jié)能策略,可以節(jié)約能耗。機床3同樣不是在0時刻開機,而是在199時刻開機以節(jié)約能耗,所加工工件依次為9、11、13、14、15、16、19、20、23以及24。由于本文將公共功率設(shè)置為20,公共能耗占總能耗比重較大,從而所求能耗最優(yōu)解中的Cmax也達到了最優(yōu)(最小)。
表6 實例25*3加工時間數(shù)據(jù)
表7 實例25*3加工功率、能耗數(shù)據(jù)
圖1 實例25*3甘特圖Fig.1 Gantt chart of instance 25*3
可見,在保證最大完工時間不變以及滿足交貨期的情況下,通過延遲開機以及關(guān)機/重啟策略可以減少機床待機能耗,達到節(jié)能的目的。對節(jié)能措施的挖掘有利于以后指導(dǎo)元啟發(fā)式方法的設(shè)計,從而可以有效求解大規(guī)模問題。
本文以能耗最小為目標,提出了5個考慮關(guān)機/重啟策略的不相關(guān)并行機MIP模型。試驗結(jié)果表明:基于不同建模思路的MIP模型尺寸復(fù)雜度、計算復(fù)雜度差別很大,基于空閑能耗的線性化MIP模型求解效果最好,在今后的應(yīng)用過程中應(yīng)優(yōu)先選用。
雖然元啟發(fā)式方法、啟發(fā)式規(guī)則等在短時間內(nèi)可以得到問題較好的解,但是即使是對于小規(guī)模問題,也不一定能夠找到最優(yōu)解,也無法證明所得到的解是最優(yōu)的,而通過MIP模型求解則可得到并證明最優(yōu)解。MIP的求解結(jié)果可以作為啟發(fā)式以及元啟發(fā)式等近似方法的參照標準。由于本文直接使用CPLEX對問題進行求解,沒有加入針對問題特性的求解算法等,因此求解效率較低,下一步研究工作將針對本文問題特性,設(shè)計相應(yīng)的求解規(guī)則以嵌入到CPLEX,提高求解效率。