有限小批量定制生產伯努利三機器裝配系統實時性能分析

劉暢 賈之陽

裝配系統是生產實踐中最基本的生產系統結構之一.系統中最終的產品通常由兩個或兩個以上組件裝配構成(例如:汽車、家電、消費電子).相比在穩態分析研究方面取得的大量成果,裝配系統的暫態過程仍然未被深入研究.這主要是由于不同零件生產線的相互作用導致了對裝配系統的分析要比傳統的串行線復雜得多.裝配系統的早期研究通常只考慮多隊列單服務器的情況,即幾種類型的零件到達一個裝配機器,從而被執行組裝操作[1].后來,Lipper等[2]和Kuo等[3]研究了有限隊列容量的三機器系統.在這些研究中,兩個服務器代表組件的生產,而另一個服務器代表裝配操作.此外,Manitz[4]對基于排隊模型的裝配系統進行了研究.具有不可靠機器和有限緩沖區容量的裝配系統的穩態性能評估在一些論文中進行了討論[5?10].

需要注意的是,雖然在過去幾十年里存在大量的串行線和裝配生產線的研究工作,但是其中的絕大部分都是假設系統運行在穩定狀態之下的.與此同時,具有有限緩沖區容量和不可靠機器的生產系統的暫態性能只在最近的一些論文中得到初步的研究,其中大部分集中在基于伯努利機器的串行線系統研究[11?14].此外,Meerkov等[15]研究了使用幾何可靠性機器模型的串行線的暫態性能.當一個生產系統根據客戶訂單或需求預測,每次只生產一個批次(或批量)的同類產品時,我們稱之為有限小批量生產運行,基于此,Jia等[16?18]近兩年研究了串行生產線和閉環生產系統的實時系統性能.

基于暫態的裝配系統的分析研究僅在以下幾篇論文中出現,Alexander等[19]研究了一類具有無限隊列容量的單機器馬爾科夫類裝配排隊系統的暫態吞吐量.Jia等[20?21]研究了復雜裝配系統的暫態性能,但是假設了系統具有無限的原材料供應量.與此同時,值得關注的是,近年來智能制造技術的發展對生產系統的暫態和動態特性研究提出了更高的要求,這對于研究相應的實時生產控制算法也至關重要.因此,本文的目的是研究有限小批量定制化生產運行下,具有有限緩沖區容量的三機裝配系統的基于暫態的性能評價.

本文余下的部分組織如下:第1節介紹了本文所研究的系統的前提假設并定義了所關注的系統性能指標.第2節詳細討論了系統數學模型的建立以及相關性能指標的推導過程.然后,提出的一種基于分解的用于近似原始系統性能的方法在第3節中被詳細介紹.所提出方法的準確性通過數值實驗進行了驗證并分析展示.最后,第4節給出了總結和未來工作展望.

1 數學模型和性能指標

1.1 數學模型

考慮如圖1所示的一個三機裝配系統,其中圓形表示機器,矩形表示緩沖區.

圖1 三機裝配系統Fig.1 Assembly production system with three Bernoulli machines

系統根據以下假設來進行定義:

1)系統的最終產品(F0)需要兩個組件.一個組件(R1)由機器m1處理.我們稱系統的這一部分(從機器m1到b1)為零件生產線1.類似地,另一個組件(R2)由機器m2處理.系統的這一部分(從機器m2到b2)稱為零件生產線2.

2)機器m0從零件生產線1和零件生產線2各取一個完成的零件裝配組成一個成品.

3)機器mi,i=0,1,2,擁有恒定且相同的周期時間τ.以一個加工周期τ為一段,將時間軸分段.所有機器在一個新的生產批次開始時運行.小批量定制生產下的每個批次具有有限的產量,每個生產批次的規模為B.每臺機器在加工完規定數量的工件后立即停止工作.

4)機器遵循伯努利可靠性模型,即,機器mi,i=0,1,2,如果既沒有被阻塞也沒有饑餓,在一個時間間隙(即加工周期)里加工處理一個工件的概率是pi,未能加工處理一個工件的概率是1?pi.參數pi∈(0,1)稱為機器mi的效率.

5)每一個在制品緩沖區,bi,i=1,2,可以用其容量Ni來表征,0

6)如果機器m0在時間間隙n內處于工作狀態,緩沖區b1或者b2在時間間隙開始時為空,則機器m0在時隙n內會饑餓.機器m1和m2在一個批次生產結束前不會出現饑餓的情況.

7)如果機器mi,i=1,2,在時間間隙n內處于工作狀態,緩沖區bi在時間間隙開始時有Ni個在制品工件,并且裝配機器m0沒能從其中取走一個工件進行處理(由于故障或源自另一條零件生產線的饑餓情況),則機器mi,i=1,2,在時隙n內被阻塞.即,加工前阻塞機制.同時假設機器m0任何時候都不會被阻塞.

注1.值得注意的是,在許多生產系統中,機器周期時間幾乎是恒定或接近恒定的.這樣的情況大多見于汽車、電子、電器等行業的生產系統.還需注意到,伯努利可靠性機器模型是適用于平均故障時間接近機器的加工周期的情況(參見使用伯努利模型為實例的文獻[22?24]).具有其他可靠性機器模型(例如:幾何型、指數型、威布爾型、對數正態型等)的生產系統將在今后的工作中進一步研究.

注2.基于批次的生產廣泛用于各種制造系統(小規模,中等規模,甚至大規模生產,單型或多類型產品生產等).一個批次有時被稱為一個分組或一個訂單.

注3.由于通常定制化生產下每個批次生產數量是有限的,整個生產過程部分或完全是在暫態下進行的.因此,嚴格來說,穩定狀態分析不再適用,而基于暫態的系統分析是必要的.

注4.上述的模型僅僅包括兩條零件生產線和一個裝配操作機器.每條零件生產線僅包含了一臺機器和一個緩沖區.每條零件生產線擁有多臺機器和緩沖區,以及擁有多條零件生產線的復雜裝配系統具有類似的假設,并且這樣的裝配系統會在未來工作中進一步研究.

1.2 性能指標

在上述定義的模型框架下,我們感興趣的性能指標包括:

1)生產率PR(n):在時間間隙n+1里,機器m0生產工件個數的期望;

2)消耗率CRi(n):在時間間隙n+1里,機器mi,i=1,2,消耗原材料個數的期望;

3)在制品庫存水平,WIPi(n):在時間間隙n里緩沖區bi,i=1,2,中的在制品個數的期望;

4)阻塞率BLi(n):機器mi,i=1,2,在時間間隙n+1里被阻塞的概率.

由于機器m0可能由于任意一條零件生產線而饑餓,我們定義機器饑餓率為:

一種通過遞歸聚合來估計這些穩態性能值的方法在文獻[22]中被提出.在本文中,我們提出了在有限量定制生產運行下評估這些暫態性能指標的方法.

此外,使ct表示機器m0完成生產B個產品的時間.將其均值表示為:

2 系統性能精確分析

2.1 性能分析

用fi(n)表示機器mi在時間間隙n結束時已經生產的工件總數量,用hi(n)表示在時間間隙n結束時緩沖區內的在制品工件數量.顯而易見,

那么,不失一般性,系統可以用一個狀態為(h1(n),h2(n),f0(n))的馬爾科夫鏈來表征,其中,

顯然,此馬爾科夫鏈的最大系統狀態數為

需要注意,有一些系統狀態是不可達到的,比如,(1,1,B),因為機器m1和m2在加工好B個工件后立刻停止了運作.換句話說,在任意一個時間間隙里,h1+f0≤B,并且h2+f0≤B.

為了計算這一馬爾科夫鏈中的狀態間轉移概率,我們首先如表1排列系統的狀態.

表1 系統狀態排序Table1 Arrangement of the system states

其中,

同時也需要注意,在每個時間間隙中,系統狀態的樣本空間是由機器23種的工作狀態所組成的.那么,

因此,在每一個時間間隔開始時,對系統的每一個可達狀態i,i∈{1,···,Q},如果,并且,可以枚舉所有的23種機器狀態的組合,根據系統動態性式(4)來確定相應的在這一時間間隔結束時的結果狀態j,j∈{1,···,Q}.然后,對于得到相同結果狀態的機器狀態組合情況,使用式(5)來計算相應的轉移概率,并將這些概率相加,最終得到一個時間間隔里,從起始的系統狀態i到結果狀態j的轉移概率.對于所有符合條件的系統狀態重復這一步驟.

系統的實時性能可以通過下式計算:

其中

其中,01,k和J1,k分別代表1×k的零矩陣和元素全為1的矩陣.與此同時,i×j維矩陣Ci×j=[Ii···Ii]表示由j/i個單位矩陣Ii組成的矩陣.

3 基于分解的性能評估

上面描述的精確分析可以擴展到更大的系統,即每個零件生產線中有多臺機器的系統.然而,隨著機器數量M,緩沖區容量,和生產規模B的增長,馬爾科夫鏈狀態的數量呈指數型增長,這將導致對大型的復雜裝配系統的分析變得不可能.因此,本節提出了一種基于分解的算法,并將其應用于三臺伯努利機器的小型裝配系統.相應的研究結果將在未來的工作中擴展到更通用的大型系統中.

3.1 基于分解的概念

文獻[8]提出一種分解方法,將原系統分解為一對串行線:上線和下線,研究了基于無限原材料供應的裝配系統的穩態性能.此外,我們以前的工作[20?21]解決了這類系統的暫態性能研究的問題.與此同時,當考慮到小批量有限量生產運行下的串行線,基于暫態的系統性能近似評估也在我們以前的工作[16?17]中進行了討論.在這一節中,我們將基于有限量生產運行下系統的暫態性能分析擴展到三臺機器的裝配系統性能分析研究中.對由多臺機器組成的零件生產線或多條零件生產線以及多個裝配操作的復雜裝配系統,將在今后的研究中進行分析.

圖2 輔助裝配系統Fig.2 Auxiliary assembly system

具體而言,引入三種輔助系統/生產線來分析此類系統.輔助裝配系統(圖2所示)首先被引入,這一輔助裝配系統具有所有原始的機器和緩沖區,但假設具有無限的原材料供應.只有原來的裝配操作中裝配機器m0處于工作狀態并且緩沖區b1非空的情況下,才可能處于工作狀態.因此,讓hi(n)表示在時間間隙n結束時緩沖區bi中的在制品零件數,可以通過式(9)估算:

圖3 輔助雙機串行線Fig.3 Auxiliary two-machine lines

用于分析這種輔助兩機串行線的方法在文獻 [11]中被提出.具體來說,用和,分別來表示輔助兩機生產線中b1和b2在時隙n結束時,有i個工件的概率.令. 根據文獻 [11],的演化可表示為:

為了研究這一輔助裝配系統的暫態性能,分析方法可以參考我們以前的工作(參閱文獻[20]).具體而言,使用輔助兩機線(圖3)來近似分析.

上生產線通過移除輔助裝配系統中的機器m2和緩沖區b2來構造.考慮到這種修改,組裝機器m0由效率隨時間變化的虛擬機器(圖3(a))來代替.同樣,下生產線可以通過移除機器m1和緩沖區b1,同時使用效率隨時間變化的虛擬機器來構造.

最后,引入有限量生產運行下的輔助單機生產線(見圖4).

圖4 輔助單機生產線Fig.4 Auxiliary one-machine lines

由于無論是m1還是m2,能夠生產一個工件的前提條件都是當且僅當它處于工作狀態且不被阻塞,同時無論或者是,能夠生產的條件是當且僅當其處于工作狀態且不會饑餓,我們定義隨時間變化的輔助單機機器效率如下:

為了分析輔助單機生產線,注意,它們每個都可由一個馬爾科夫鏈來表征,其中,系統狀態為已被這臺機器加工過的工件數量(參閱文獻[17]).讓,其中表示在時隙n結束時已經加工了j個工件的概率.的演化可以通過以下線性時變方程給出:

其中初始狀態是

其中初始狀態是

綜上,為了分析圖1中的有限量運行下的三機裝配線的暫態性能,我們將原始系統的動態特性進行分解和簡化,通過分析一系列分解后相互影響的動態特性更加簡單的系統,來近似評估原始系統的實時性能.具體來說,對于原始系統(圖1),其動態特性包括兩方面:緩沖區中在制品數量的演化和在每臺機器上已完成加工處理的工件數量.首先引入使用原始系統機器和緩沖區參數的輔助裝配系統,同時假設無限原材料(圖2).在這個系統中,我們只關注系統中的緩沖區在制品數量的演化.為了分析圖2所示系統,進一步引入輔助雙機串行線(圖3),其中,為了考慮移除相應機器和緩沖區所帶來的影響,上生產線和下生產線中裝配機器所在位置分別使用相應的參數時變的虛擬機器來替代.因此,通過分析輔助雙機串行線(圖3),事實上可以得到輔助裝配系統(圖2)中系統狀態(緩沖區在制品數量)的實時分布情況.最后,引入輔助單機生產線(圖4)來分析在相應機器上完成加工處理工件數量的動態特性.而每一臺單機生產線的時變參數都是在考慮了輔助雙機串行線中的系統狀態的影響下,近似推導得出的.

3.2 性能評估近似公式

基于上述構造的輔助生產線或生產系統,我們提出了近似原系統性能指標的計算公式.首先,有限量生產運行下一個批次的生產完成時間通過使用輔助虛擬單機線中的任意一個來近似估算.不失一般性,使用,同時令表示原始系統中機器m0在時隙n結束時處理加工完整個批次所有工件的近似概率.那么,

其次,原系統的生產率和各個零件生產線的消耗率可由輔助單機生產線的生產率來近似:

為了估算WIPi(n),BLi(n)和ST0,i(n),兩種輔助生產線需要結合起來.具體來說,這些性能評估使用相應的輔助兩機生產線來近似估算,同時考慮相對應的機器在輔助單機生產線上還沒有完成加工整個批次所有產品的概率:

最后,一個批次的完成時間期望可以被近似為:

其中,T滿足以下條件:

綜上,基于分解的計算方法流程圖如圖5所示.

圖5 分解算法流程圖Fig.5 Flow chart of the calculation procedure

對于所提出的性能近似的方法的精確程度,我們通過對10000條參數隨機而均勻地從式(21)所示的集合或者區間中選取的三機伯努利裝配系統,進行基于精確解析和基于分解的近似性能評估分析,來驗證所提近似方法的精確性.

對于每一條參數隨機產生的裝配系統,我們分別通過精確分析式(8)和基于分解的近似分析式(17)～(20)來計算其各項性能指標.結果顯示,對于這10000條裝配系統的各項性能指標的平均相對誤差,它們的中值都在1%以下.

作為一個例子,考慮圖6所顯示的裝配系統.每個機器(圓形表示)上的數字表示其效率,而每個緩沖區(矩形)中的數字表示其容量.這些參數是隨機生成的.在本例中,所有緩沖區都被假設在起始狀態時是空的.首先需要注意的是,使用精確分析方法,根據式(2),系統的狀態數量為1620;而經過分解后,我們只需要分析六個相對較小但相互影響的系統:一條雙機上生產線,一條雙機下生產線,兩條上單機生產線,兩條下單機生產線.六個較小的馬爾科夫鏈的總狀態數為333.在保證精確度的基礎上,相較精確分析,基于分解的近似分析使系統狀態數量有了極大的降低.與此同時,從計算時間的角度來看,使用MATLAB軟件在同一臺電腦配置為因特爾酷睿i7-6700的CPU和16GB的RAM上,基于精確分析和基于分解的近似分析,所需要的運算時間分別為13.35秒和0.11秒,近似算法在計算高效性上也顯示出了極大的優勢.系統的暫態性能如圖7所示,從圖中可以看出,整個生產運行過程分為三個階段.在第一階段,產品開始進入空系統.在此期間,生產率和在制品數量都從0上升到穩態值.同時,由于更多的工件進入系統,零件生產線1(或者零件生產線2)的消耗率從p1(或者p2)開始逐漸減小.在第二階段中,系統運行接近穩定狀態,所有暫態性能指標都或多或少地處于平穩狀態.最后,當生產運行接近完成時,所有性能指標開始下降,最終達到0.基于該分解算法的高精度也可以從圖中清晰地看到.需要注意的是,雖然精確的分析在這種小型裝配系統中仍然可以被推導出來,然而隨著系統參數的增長,精確分析也變得越來越不可能實現.基于分解思想的性能近似評估方法的計算高效性將在這樣的大型裝配系統中體現出來.深入的相關研究將在未來的工作中被進一步討論.

圖6 三機伯努利裝配系統的數值實例Fig.6 Example of an assembly system with three Bernoulli machines

4 結束語

本文研究了具有三臺伯努利機器,有限緩沖區容量和有限量生產運行下的裝配系統的暫態性能評估問題.具體地,首先推導了系統性能評價的精確數學模型和解析公式.然后,提出了一種基于分解的性能評估算法,通過將系統轉換成一系列相互作用的輔助串行線來近似評估原始系統的暫態性能.論文推導了基于分解的三機裝配系統實時性能估計公式,并通過數值實驗驗證了算法的準確性和計算高效性.

圖7 分解近似與精確分析的三機伯努利裝配系統暫態實時性能評估對比Fig.7 Comparison of decomposition-based approxiamtion and exact analysis for transientperformance evaluation in assembly system with three Bernoulli machines

今后在這方面的工作包括將算法擴展到每個零件生產線具有多臺機器和多個緩沖區的系統,或多條零件生產線和多裝配操作的復雜裝配系統.此外,還會將研究結果推廣到具有其他機器可靠性模型(幾何型、指數型、威布爾型等)的裝配系統中.