基于排隊網絡模型的混流制造系統服務率優化*

張巖巖，白金花，武 福，李忠學

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

1 引 言

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

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

2 排隊理論及排隊網絡

排隊理論又稱為隨機服務系統理論，它通過對服務對象的到達模式及服務時間的統計分析，得到等待時間、排隊長度、忙期長短等相關參數的統計規律，根據這些參數的統計規律來改進服務系統的結構或者重組服務系統，使系統的特定指標達到最優。

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

3 系統建模

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

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

圖1 不同類型產品的聚合

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

(1)

(2)

(3)

其中:式(1)可由r類產品在排隊網絡中的轉移概率矩陣確定。

產品在系統中逗留會產生一定的損失費用，提高工作單元的服務率雖然可以減少產品的逗留時間，但也意味著運行成本的增加，不利于企業實現高效益的目標，如圖3所示。

圖3 排隊系統的期望費用圖

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

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

(4)

約束條件：

μir>λir

(5)

其中:Fsi為工作單元i單位時間內提高一個單位的服務率時所增加的服務費用，Fwi為混合產品在工作單元i內逗留單位時間的損失費用，μi為工作單元i對混合產品的平均服務率，E(Li)為工作單元i單位時間內的平均產品數。

4 模型分析

4.1 E(μi)的計算

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

(6)

因此，得到工作單元i對混合產品的平均服務率為：

(7)

4.2 E(Li)的計算

Little定理 設產品按照泊松流到達系統，參數為λ，服務機構的服務時間服從參數為μ的負指數分布。排隊系統中單位時間內的平均產品數為E(L)，平均逗留時間為E(T)，那么兩者之間的關系為：

E(L)=λE(T)

(8)

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

(9)

(10)

由于工作單元對產品的服務規則是先到先服務，任何產品都不具有優先權，所有的產品在同一個排隊隊列中等待接受服務，因此混合產品在工作單元i內的平均逗留時間為：

(11)

根據Little定理，得到工作單元i單位時間內的平均產品數為：

(12)

1-qki}]

(13)

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

(14)

(15)

(16)

根據變異系數的定義得到工作單元i對r類產品加工時間的變異系數平方為：

(17)

(18)

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

(19)

則工作單元i對混合產品加工時間的變異系數平方為：

(20)

因此，最終的目標函數可以記為：

(21)

約束條件：

μir>λir

(22)

4.5 E(Li)的簡化處理

(23)

因此，目標函數式(21)可以修正記為：

(24)

約束條件：

μir>λir

(25)

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

5 算 例

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

圖4 兩類產品的加工路徑

現假設兩類產品在各個工作單元接受的服務時間服從負指數分布且所有的到達過程和服務過程相互獨立。為了在減小產品逗留可能性的同時又不至于因其服務率過高而增加運行成本，對此混流制造系統各個工作單元的服務率進行控制優化。表1所示為每一個工作單元單位時間內每一類產品的逗留費用與服務費用。

表1 每個工作單元單位時間內每類產品 的逗留費用與服務費用

5.1 確定混合產品訪問每個工作單元平均到達率

(1) 兩類產品在排隊網絡中的轉移概率矩陣P1和P2

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

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

根據式(1)，結合轉移概率矩陣Pr可以得到兩類產品訪問每一個工作單元的平均到達率為：

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

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

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

(0.0033，0.0031，0.0031，0.0035)

5.2 確定混合產品在制造系統中的轉移概率矩陣P

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

P=

5.3 確定每一個工作單元對混合產品加工時間的變異系數平方

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

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

在本算例中，每一類產品在系統中各個工作單元接受的服務時間服從負指數分布，因此：

=(1.0000，1.0000，1.0000，1.0000)

5.4 確定混合產品到達每一個工作單元的平均間隔時間的變異系數平方

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

(0.0033，0.0031，0.0031，0.0035)

由于每一類產品都是按照泊松流到達制造系統的，而相互獨立的泊松流的合成流仍然為泊松流，因此：

(1.0000，1.0000，1.0000，1.0000)

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

(1.0056，0.9946，0.9999，0.9996)

5.5 建立目標函數

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

minF(μ11，μ12)

(26)

約束條件：

μ11>0.0016

(27)

μ12>0.0017

(28)

5.6 目標函數求解算法設計

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

(1) 內點懲罰函數法的實現步驟

(29)

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

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

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

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

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

(2) 坐標輪換法的實現步驟

(3) 黃金分割法的實現步驟

① 在進退法確定的初始搜索區間[a，b]內取兩個試算點α1和α2，計算函數值φ(α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)

在新區間內取一個新點α1，計算函數值φ(α1)：

α1=a+0.382(b-a)

(32)

以縮小后的新區間作為原區間，重復步驟②。

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

置換符號：

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

在新區間內取一個新點α2，計算函數值φ(α2)：

α2=a+0.618(b-a)

(33)

以縮小后的新區間作為原區間，重復步驟②。

(4) 進退法的實現步驟

① 給定初始點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))，轉向步驟④，否則轉向步驟⑤；

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

⑤ 若k2=1，轉向步驟⑥，否則轉向步驟⑦；

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

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

5.7 求解目標函數

根據內點懲罰函數法，引入懲罰因子，將目標函數式(26)轉化為無約束優化問題。構造的懲罰函數為：

(34)

算法程序由MATLAB語言編寫，求解過程中的所有程序在MATLAB7.5.0環境中運行。對工作單元1服務率優化的最終結果為：

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

同樣地，可以求出其他工作單元加工每一類產品的最優服務率分別為：

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

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

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

即完成對算例中混流制造系統的服務率優化工作。

6 結 語

應用排隊理論對混流制造系統的服務率優化問題建立了物理模型和數學模型，并利用最優化方法和MATLAB程序對所建數學模型進行了優化求解。在此過程中，將混流制造系統等效為了開放式排隊網絡模型，每一類產品到達網絡的時間間隔以及在各個工作單元接受的服務時間都是服從負指數分布的。事實上，在每一類產品的到達時間間隔以及服務時間都服從一般分布的情形下，文中所建模型和相關計算過程同樣適用。算例計算結果表明了應用排隊理論對混流制造系統進行服務率優化的方法是可行有效的，從而為制造系統服務率優化問題的解決提供了一種新的思路。

參考文獻：

[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].控制理論與應用，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.