基于排隊網(wǎng)絡(luò)模型的混流制造系統(tǒng)服務(wù)率優(yōu)化*

張巖巖，白金花，武 福，李忠學(xué)

(蘭州交通大學(xué) 機(jī)電工程學(xué)院，甘肅 蘭州 730070)

1 引 言

隨著市場競爭激烈化以及顧客需求個性化的增加，混流生產(chǎn)方式成為制造業(yè)普遍采用的一種生產(chǎn)組織方式，它可在基本不改變設(shè)施布局的前提下，在同一制造系統(tǒng)內(nèi)加工出多種不同型號和數(shù)量的產(chǎn)品，具有很高的靈活性，在汽車、家電、玩具等產(chǎn)品制造企業(yè)中有著廣闊的應(yīng)用前景。擁有較好的服務(wù)率水平，是制造企業(yè)鞏固和發(fā)展自身優(yōu)勢的一種持久能力，也是其最為關(guān)心和追求的效用指標(biāo)，如何對制造系統(tǒng)的服務(wù)率水平進(jìn)行優(yōu)化也就成為人們普遍關(guān)心的問題。但在實(shí)際生產(chǎn)過程中，制造系統(tǒng)的服務(wù)率水平受到許多外部及內(nèi)部復(fù)雜因素的制約，影響著生產(chǎn)計劃、成本、供求關(guān)系等。在混流制造系統(tǒng)的實(shí)際運(yùn)行中，不同類型的產(chǎn)品同時投入生產(chǎn)：對于同一類型的產(chǎn)品來說，不同工作單元的平均加工時間不同；而對于不同類型的產(chǎn)品來說，同一工作單元的平均加工時間又存在差異。因此對混流制造系統(tǒng)，服務(wù)率不合理造成資源浪費(fèi)以及生產(chǎn)效率低下的問題較突出。

排隊是制造系統(tǒng)中的一個常見現(xiàn)象。排隊理論應(yīng)用廣泛，是研究一切服務(wù)系統(tǒng)的有效工具，這使得利用排隊理論對混流制造系統(tǒng)優(yōu)化問題進(jìn)行研究成為可能。在已有的基于排隊理論優(yōu)化服務(wù)率的研究中， Grassmann，Chen和Kashyap[1]對具有狀態(tài)無關(guān)到達(dá)率的M/G/1排隊模型進(jìn)行了研究，提出了確定其在穩(wěn)定狀態(tài)下最優(yōu)服務(wù)率的方法；Jo[2]利用拉格朗日方法得到了開放式Jackson排隊網(wǎng)絡(luò)中的每一個節(jié)點(diǎn)在相鄰兩個穩(wěn)定狀態(tài)下最優(yōu)服務(wù)率的相互關(guān)系，之后給出了計算每個節(jié)點(diǎn)最優(yōu)服務(wù)率的有效算法；Song，Xing和Sun[3]研究了隨機(jī)需求不可靠制造系統(tǒng)的最優(yōu)服務(wù)率控制問題等。而基于排隊理論對更為復(fù)雜的混流制造系統(tǒng)進(jìn)行服務(wù)率優(yōu)化的相關(guān)研究則顯得不足。鑒于混流制造系統(tǒng)的效益性、現(xiàn)實(shí)性以及排隊理論的實(shí)用性、有效性，本文應(yīng)用排隊理論對混流制造系統(tǒng)的服務(wù)率優(yōu)化問題進(jìn)行研究。

2 排隊理論及排隊網(wǎng)絡(luò)

排隊理論又稱為隨機(jī)服務(wù)系統(tǒng)理論，它通過對服務(wù)對象的到達(dá)模式及服務(wù)時間的統(tǒng)計分析，得到等待時間、排隊長度、忙期長短等相關(guān)參數(shù)的統(tǒng)計規(guī)律，根據(jù)這些參數(shù)的統(tǒng)計規(guī)律來改進(jìn)服務(wù)系統(tǒng)的結(jié)構(gòu)或者重組服務(wù)系統(tǒng)，使系統(tǒng)的特定指標(biāo)達(dá)到最優(yōu)。

一個排隊系統(tǒng)中存在多項服務(wù)，每項服務(wù)都有自己的服務(wù)機(jī)構(gòu)。在這種情況下，系統(tǒng)中也就存在多個排隊隊列，整個排隊系統(tǒng)便形成排隊網(wǎng)絡(luò)。排隊網(wǎng)絡(luò)是一個比較復(fù)雜的服務(wù)系統(tǒng)，但其應(yīng)用廣泛。排隊網(wǎng)絡(luò)研究的首個突破是20世紀(jì)五六十年代的Jackson網(wǎng)絡(luò)模型[4]，使得計算機(jī)系統(tǒng)建模工作得到了很大的簡化。隨后，針對更加復(fù)雜的系統(tǒng)建模，如引入不同類型的顧客和新的服務(wù)規(guī)則，Baskett，Chandy，Muntz和Palacios[5]等學(xué)者提出了BCMP定理。而對于更加一般的排隊網(wǎng)絡(luò)，由于對其性能參數(shù)的評價難以獲得準(zhǔn)確解，許多學(xué)者研究并提出了一些如均值分析、分解等近似方法。本文則基于開放式排隊網(wǎng)絡(luò)理論對混流制造系統(tǒng)建立物理模型和數(shù)學(xué)模型。

3 系統(tǒng)建模

在混流制造系統(tǒng)中，各種類型的產(chǎn)品混合連續(xù)地到達(dá)各個工作單元接受服務(wù)，不同類型的產(chǎn)品在各個工作單元接受的服務(wù)時間不同。依據(jù)開放式排隊網(wǎng)絡(luò)模型的特征(顧客由外部進(jìn)入網(wǎng)絡(luò)，接受完所有的預(yù)定服務(wù)后離開網(wǎng)絡(luò))，一個混流制造系統(tǒng)可以等效為一個開放式排隊網(wǎng)絡(luò)模型。

對于混流排隊網(wǎng)絡(luò)模型，在計算過程中不同類型的產(chǎn)品之間相互干擾，使得單獨(dú)針對每一類產(chǎn)品求解排隊系統(tǒng)的相關(guān)指標(biāo)之后再線性相加的策略成為不可能，因此在建立數(shù)學(xué)模型之前，聚合到達(dá)每一個工作單元的所有各類產(chǎn)品為一類混合產(chǎn)品，并將其看作單一類型的產(chǎn)品，如圖1所示，則整個混流排隊網(wǎng)絡(luò)就簡化成為了單一產(chǎn)品類型的排隊網(wǎng)絡(luò)，其中的每一個工作單元都等效為G/G/1排隊模型?？紤]工作單元i，進(jìn)而可以根據(jù)其典型隊列模型圖和式(1)～(3)求得每一類產(chǎn)品和混合產(chǎn)品到達(dá)每一個工作單元的平均到達(dá)率λir和λi以及混合產(chǎn)品在制造系統(tǒng)中的轉(zhuǎn)移概率pij。圖2所示為i工作單元的典型隊列模型圖。

圖1 不同類型產(chǎn)品的聚合

圖2 i工作單元的典型隊列模型圖

(1)

(2)

(3)

其中:式(1)可由r類產(chǎn)品在排隊網(wǎng)絡(luò)中的轉(zhuǎn)移概率矩陣確定。

產(chǎn)品在系統(tǒng)中逗留會產(chǎn)生一定的損失費(fèi)用，提高工作單元的服務(wù)率雖然可以減少產(chǎn)品的逗留時間，但也意味著運(yùn)行成本的增加，不利于企業(yè)實(shí)現(xiàn)高效益的目標(biāo)，如圖3所示。

圖3 排隊系統(tǒng)的期望費(fèi)用圖

綜合考慮上述因素，如果取目標(biāo)函數(shù)F(μi1，μi2，…，μiR)為工作單元i單位時間內(nèi)的服務(wù)費(fèi)用與產(chǎn)品逗留損失費(fèi)用之和的期望值，則基于排隊理論建立混流制造系統(tǒng)的期望費(fèi)用模型如下：

minF(μi1，μi2，…，μiR)=Fsiμi+FwiE(Li)

(4)

約束條件：

μir>λir

(5)

其中:Fsi為工作單元i單位時間內(nèi)提高一個單位的服務(wù)率時所增加的服務(wù)費(fèi)用，F(xiàn)wi為混合產(chǎn)品在工作單元i內(nèi)逗留單位時間的損失費(fèi)用，μi為工作單元i對混合產(chǎn)品的平均服務(wù)率，E(Li)為工作單元i單位時間內(nèi)的平均產(chǎn)品數(shù)。

4 模型分析

4.1 E(μi)的計算

當(dāng)混流制造系統(tǒng)處于穩(wěn)態(tài)時，工作單元i加工r類產(chǎn)品的概率和平均加工時間分別為λir/λi和E(tir)，則工作單元i對混合產(chǎn)品的平均加工時間E(ti)為：

(6)

因此，得到工作單元i對混合產(chǎn)品的平均服務(wù)率為：

(7)

4.2 E(Li)的計算

Little定理 設(shè)產(chǎn)品按照泊松流到達(dá)系統(tǒng)，參數(shù)為λ，服務(wù)機(jī)構(gòu)的服務(wù)時間服從參數(shù)為μ的負(fù)指數(shù)分布。排隊系統(tǒng)中單位時間內(nèi)的平均產(chǎn)品數(shù)為E(L)，平均逗留時間為E(T)，那么兩者之間的關(guān)系為：

E(L)=λE(T)

(8)

對于一般服務(wù)時間分布或一般到達(dá)的排隊系統(tǒng)，并非在任何時刻都有馬爾可夫性，對其隊列長度的平穩(wěn)分布進(jìn)行求解相當(dāng)困難，也沒有明顯可利用的結(jié)果。由于G/G/1排隊模型中僅知道均值和變異系數(shù)的平方兩個參數(shù)，本文利用Kingman[6-7]近似式(9)計算混合產(chǎn)品在接受服務(wù)前的平均等待時間E(Tqi)：

(9)

(10)

由于工作單元對產(chǎn)品的服務(wù)規(guī)則是先到先服務(wù)，任何產(chǎn)品都不具有優(yōu)先權(quán)，所有的產(chǎn)品在同一個排隊隊列中等待接受服務(wù)，因此混合產(chǎn)品在工作單元i內(nèi)的平均逗留時間為：

(11)

根據(jù)Little定理，得到工作單元i單位時間內(nèi)的平均產(chǎn)品數(shù)為：

(12)

1-qki}]

(13)

ωi=[1+4(1-ρi)2(υi-1)]-1

(14)

(15)

(16)

根據(jù)變異系數(shù)的定義得到工作單元i對r類產(chǎn)品加工時間的變異系數(shù)平方為：

(17)

(18)

因此，工作單元i對混合產(chǎn)品加工時間的二次矩為：

(19)

則工作單元i對混合產(chǎn)品加工時間的變異系數(shù)平方為：

(20)

因此，最終的目標(biāo)函數(shù)可以記為：

(21)

約束條件：

μir>λir

(22)

4.5 E(Li)的簡化處理

(23)

因此，目標(biāo)函數(shù)式(21)可以修正記為：

(24)

約束條件：

μir>λir

(25)

選擇合適的方法求解目標(biāo)函數(shù)式(24)，就使得計算過程在很大程度上得到簡化。

5 算 例

一混流制造系統(tǒng)同時生產(chǎn)2種類型的產(chǎn)品，r(r=1，2)類產(chǎn)品按照參數(shù)為λar=0.0014(件/s)的泊松過程從外部到達(dá)制造系統(tǒng)，加工路徑如圖4所示[11]，其中圓形節(jié)點(diǎn)表示制造系統(tǒng)中的工作單元，節(jié)點(diǎn)中的數(shù)字標(biāo)注表示工作單元編號，工作單元之間轉(zhuǎn)移路徑上的數(shù)字標(biāo)注表示產(chǎn)品在一個工作單元服務(wù)完成后進(jìn)入下一個工作單元接受服務(wù)的概率。

圖4 兩類產(chǎn)品的加工路徑

現(xiàn)假設(shè)兩類產(chǎn)品在各個工作單元接受的服務(wù)時間服從負(fù)指數(shù)分布且所有的到達(dá)過程和服務(wù)過程相互獨(dú)立。為了在減小產(chǎn)品逗留可能性的同時又不至于因其服務(wù)率過高而增加運(yùn)行成本，對此混流制造系統(tǒng)各個工作單元的服務(wù)率進(jìn)行控制優(yōu)化。表1所示為每一個工作單元單位時間內(nèi)每一類產(chǎn)品的逗留費(fèi)用與服務(wù)費(fèi)用。

表1 每個工作單元單位時間內(nèi)每類產(chǎn)品 的逗留費(fèi)用與服務(wù)費(fèi)用

5.1 確定混合產(chǎn)品訪問每個工作單元平均到達(dá)率

(1) 兩類產(chǎn)品在排隊網(wǎng)絡(luò)中的轉(zhuǎn)移概率矩陣P1和P2

本算例中，轉(zhuǎn)移概率矩陣Pr(r=1，2)中的元素pijr(i=1，2，…，4；j=1，2，…，4)表示r類產(chǎn)品在工作單元i服務(wù)完成后進(jìn)入工作單元j接受服務(wù)的概率，p0ir表示r類產(chǎn)品從網(wǎng)絡(luò)外部到達(dá)工作單元i接受服務(wù)的概率，pi5r表示r類產(chǎn)品在工作單元i服務(wù)完成后離開網(wǎng)絡(luò)的概率。

(2) 計算兩類產(chǎn)品訪問每一個工作單元的平均到達(dá)率

根據(jù)式(1)，結(jié)合轉(zhuǎn)移概率矩陣Pr可以得到兩類產(chǎn)品訪問每一個工作單元的平均到達(dá)率為：

(λ11，λ21，λ31，λ41)=(0.0016，0.0012，0.0015，0.0000)

(λ12，λ22，λ32，λ42)=(0.0017，0.0019，0.0016，0.0035)

因此混合產(chǎn)品訪問每一個工作單元的平均到達(dá)率為：(λ1，λ2，λ3，λ4)=

(0.0033，0.0031，0.0031，0.0035)

5.2 確定混合產(chǎn)品在制造系統(tǒng)中的轉(zhuǎn)移概率矩陣P

由式(3)計算可得混合產(chǎn)品在制造系統(tǒng)中的轉(zhuǎn)移概率矩陣P為：

P=

5.3 確定每一個工作單元對混合產(chǎn)品加工時間的變異系數(shù)平方

在滿足約束μir>λir的前提下，為方便計算，選擇μir的初始值為：

(2) 計算每一個工作單元對混合產(chǎn)品加工時間的變異系數(shù)平方

在本算例中，每一類產(chǎn)品在系統(tǒng)中各個工作單元接受的服務(wù)時間服從負(fù)指數(shù)分布，因此：

=(1.0000，1.0000，1.0000，1.0000)

5.4 確定混合產(chǎn)品到達(dá)每一個工作單元的平均間隔時間的變異系數(shù)平方

(ρ1，ρ2，ρ3，ρ4)=

(0.0033，0.0031，0.0031，0.0035)

由于每一類產(chǎn)品都是按照泊松流到達(dá)制造系統(tǒng)的，而相互獨(dú)立的泊松流的合成流仍然為泊松流，因此：

(1.0000，1.0000，1.0000，1.0000)

由式(13)～(16)可計算出混合產(chǎn)品到達(dá)每一個工作單元的平均間隔時間的變異系數(shù)平方為：

(1.0056，0.9946，0.9999，0.9996)

5.5 建立目標(biāo)函數(shù)

將上面所得參數(shù)代入式(24)就可得到各個工作單元最小化期望費(fèi)用的目標(biāo)函數(shù)。以工作單元1為例，其目標(biāo)函數(shù)為：

minF(μ11，μ12)

(26)

約束條件：

μ11>0.0016

(27)

μ12>0.0017

(28)

5.6 目標(biāo)函數(shù)求解算法設(shè)計

目標(biāo)函數(shù)式(26)是一個不等式約束非線性規(guī)劃問題，可利用內(nèi)點(diǎn)懲罰函數(shù)法將原約束優(yōu)化問題轉(zhuǎn)化為無約束優(yōu)化問題，然后選取坐標(biāo)輪換法對目標(biāo)函數(shù)進(jìn)行求解，在坐標(biāo)輪換法的每一輪迭代過程中利用黃金分割法確定各個步長，而黃金分割法實(shí)施的初始搜索區(qū)間[a，b]則通過進(jìn)退法確定。

(1) 內(nèi)點(diǎn)懲罰函數(shù)法的實(shí)現(xiàn)步驟

(29)

② 在可行域之內(nèi)選取初始點(diǎn)x(0)，令k=1；

③ 選取適當(dāng)大小的初始懲罰因子r(0)；

④ 從x(k-1)點(diǎn)出發(fā)，利用坐標(biāo)輪換法求解minφ(x，r)的極小點(diǎn)xk1(r(k))；

⑤ 判斷精度|xk1(r(k))-xk1(r(k-1)) |≤ε1，若滿足，停止迭代計算，并以xk1(r(k))作為原問題的近似極小值，否則轉(zhuǎn)向步驟⑥；

⑥ 取r(k+1)=cr(k)(c≈0.1～0.7)，x(k-1)=xk1(r(k))，k=k+1，轉(zhuǎn)向步驟④。

(2) 坐標(biāo)輪換法的實(shí)現(xiàn)步驟

(3) 黃金分割法的實(shí)現(xiàn)步驟

① 在進(jìn)退法確定的初始搜索區(qū)間[a，b]內(nèi)取兩個試算點(diǎn)α1和α2，計算函數(shù)值φ(α1)和φ(α2)：

α1=a+0.382(b-a)

(30)

α2=a+0.618(b-a)

(31)

② 比較φ(α1)和φ(α2)

第一種情況，φ(α1)<φ(α2)：

置換符號：

α2→b，α1→α2，φ(α1)→φ(α2)

在新區(qū)間內(nèi)取一個新點(diǎn)α1，計算函數(shù)值φ(α1)：

α1=a+0.382(b-a)

(32)

以縮小后的新區(qū)間作為原區(qū)間，重復(fù)步驟②。

第二種情況，φ(α1)>φ(α2)：

置換符號：

α1→a，α2→α1，φ(α2)→φ(α1)

在新區(qū)間內(nèi)取一個新點(diǎn)α2，計算函數(shù)值φ(α2)：

α2=a+0.618(b-a)

(33)

以縮小后的新區(qū)間作為原區(qū)間，重復(fù)步驟②。

(4) 進(jìn)退法的實(shí)現(xiàn)步驟

① 給定初始點(diǎn)a(0)和初始步長h0，令h0→h，a(0)→a(1)，0→k2；

② 令a(4)=a(1)+h，k2+1→k2；

③ 比較φ(a(4))和φ(a(1))，若φ(a(4))<φ(a(1))，轉(zhuǎn)向步驟④，否則轉(zhuǎn)向步驟⑤；

④ 令a(1)→a(2)，a(4)→a(1)，φ(a(1))→φ(a(2))，φ(a(4))→φ(a(1))，h=2h，轉(zhuǎn)向步驟②；

⑤ 若k2=1，轉(zhuǎn)向步驟⑥，否則轉(zhuǎn)向步驟⑦；

⑥ 令h=-h，a(4)→a(2)，φ(a(4))→φ(a(2))，轉(zhuǎn)向步驟②；

⑦ 令a(2)→a(3)，a(1)→a(2)，a(4)→a(1)，終止計算，極小點(diǎn)所在的兩個可能初始搜索區(qū)間為[a，b]=[a(1)，a(3)]或[a，b]=[a(3)，a(1)]。

5.7 求解目標(biāo)函數(shù)

根據(jù)內(nèi)點(diǎn)懲罰函數(shù)法，引入懲罰因子，將目標(biāo)函數(shù)式(26)轉(zhuǎn)化為無約束優(yōu)化問題。構(gòu)造的懲罰函數(shù)為：

(34)

算法程序由MATLAB語言編寫，求解過程中的所有程序在MATLAB7.5.0環(huán)境中運(yùn)行。對工作單元1服務(wù)率優(yōu)化的最終結(jié)果為：

(μ11，μ12)=(0.0 288，0.0 358)

同樣地，可以求出其他工作單元加工每一類產(chǎn)品的最優(yōu)服務(wù)率分別為：

(μ21，μ22)=(0.0 303，0.0 314)

(μ31，μ32)=(0.0 340，0.0 364)

(μ41,μ42)=(0，0.0 314)

即完成對算例中混流制造系統(tǒng)的服務(wù)率優(yōu)化工作。

6 結(jié) 語

應(yīng)用排隊理論對混流制造系統(tǒng)的服務(wù)率優(yōu)化問題建立了物理模型和數(shù)學(xué)模型，并利用最優(yōu)化方法和MATLAB程序?qū)λ〝?shù)學(xué)模型進(jìn)行了優(yōu)化求解。在此過程中，將混流制造系統(tǒng)等效為了開放式排隊網(wǎng)絡(luò)模型，每一類產(chǎn)品到達(dá)網(wǎng)絡(luò)的時間間隔以及在各個工作單元接受的服務(wù)時間都是服從負(fù)指數(shù)分布的。事實(shí)上，在每一類產(chǎn)品的到達(dá)時間間隔以及服務(wù)時間都服從一般分布的情形下，文中所建模型和相關(guān)計算過程同樣適用。算例計算結(jié)果表明了應(yīng)用排隊理論對混流制造系統(tǒng)進(jìn)行服務(wù)率優(yōu)化的方法是可行有效的，從而為制造系統(tǒng)服務(wù)率優(yōu)化問題的解決提供了一種新的思路。

參考文獻(xiàn)：

[1] Grassmann W K，Chen X, Kashyap B R K. Optimal Service Rates for the State-Dependent M/G/1 Queues in Steady State[J]. Operations Research Letters，2001，29(2)：57-63.

[2] Jo K Y. A Lagrangian Algorithm for Computing the Optimal Service Rates in Jackson Queuing networks[J]. Computers and Operations Research，1989，16(5)：431-440.

[3] Song D P，Xing W, Sun Y X. Optimal Service Rate Allocation Policy of an Unreliable Manufacturing System with Random Demands[J].控制理論與應(yīng)用，1998，15(4)：621-626.

[4] Jackson J R. Networks of Waiting Lines[J]. Operations Research，1957，5(4)：518-521.

[5] Baskett F，Chandy K M，Muntz R R, et al. Open，Closed，and Mixed Networks of Queues with Different Classes of Customers[J].Journal of the Association for Computing Machinery，1975，22(2)：248-260.

[6] Kingman J F C. On Queues in Heavy Traffic[J]. Journal of the Royal Statistical Society. Series B：Methodological，1962，24(2)：383-392.

[7] Kingman J F C. The Heavy Traffic Approximation in the Theory of Queues[C]. Proceedings of the Symposium on Congestion Theory. Chapel Hill：The University of North Carolina Press，1964：137-169.

[8] Whitt W. The Queueing Network Analyzer[J]. The Bell System Technical Journal，1983，62(9)：2779-2815.

[9] Bitran G R, Tirupati D. Tradeoff curves，Targeting and Balancing in Manufacturing Queueing Networks[J].Operations Research，1989，37(4)：547-564.

[10] Bitran G R, Morabito R. State-of-the-art survey：Open Queueing Networks：Optimization and Performance Evaluation Models for Discrete Manufacturing Systems[J]. Production and Operations Management，1996，5(2)：163-193.

[11] Curry G L, Feldman R M. Manufacturing Systems Modeling and Analysis[M]. USA: Springer，2010.