蔣仁言, 張碧雯

(長沙理工大學 汽車與機械工程學院,湖南 長沙 410114)

0 引言

威布爾分布是建模非負、連續(xù)隨機變量使用最廣泛的分布模型[1]；威布爾分布的更新函數(shù)(簡稱威布爾更新函數(shù))在可靠性和可用度分析、維修決策優(yōu)化、質(zhì)保費用分析、備件庫存管理、排隊系統(tǒng)、交通流分析等領(lǐng)域有廣泛的應(yīng)用[2～8]。威布爾更新函數(shù)沒有解析表達式，這給解各種涉及更新函數(shù)的優(yōu)化問題帶來極大不便。因此，開發(fā)威布爾更新函數(shù)的近似式已經(jīng)吸引了持續(xù)的注意。

根據(jù)近似式在時間軸t的適用范圍，可將各種威布爾更新函數(shù)的近似式分為以下兩大類。第一類近似式適合于中偏小的t值[9]，另一類近似式適合于整個t值范圍(即0到∞)。本文著重討論第二類近似式。

根據(jù)模型的構(gòu)建方式，第二類近似式可進一步分為三個子類：分段函數(shù)近似[10～13]，級數(shù)近似[14～16]和混合模型形式的近似[9,17]。一個好的近似式應(yīng)同時滿足精確性和簡單性兩個要求。分段模型在結(jié)構(gòu)上簡單，但其精度通常偏低。對于某些分布函數(shù)(如正態(tài)分布和伽瑪分布)，級數(shù)模型的精度隨項數(shù)的增加而增大，但過多的項數(shù)導(dǎo)致近似更復(fù)雜[18]?；旌闲问降慕撇捎靡粋€權(quán)函數(shù)將兩個極限關(guān)系平滑地連接成一個函數(shù)；在結(jié)構(gòu)上相對地簡單，且其精度比分段函數(shù)近似要高，即此類近似能較好地滿足精確性和簡單性要求。但是，現(xiàn)有混合形式的近似當威布爾形狀參數(shù)較大時，其精度仍有一定的改進空間[17]。

為了改進現(xiàn)有近似式在形狀參數(shù)較大時精度欠高的問題，本文提出一個新的混合形式的威布爾更新函數(shù)近似式，可供威布爾形狀參數(shù)是大的(>3.65)情況時使用。此外，通過對比現(xiàn)有各近似式在不同形狀參數(shù)取值范圍內(nèi)的精度，識別好的近似式，確定其適用范圍，從而確立一個基于形狀參數(shù)的近似式選擇方法。

本文結(jié)構(gòu)如下。第1節(jié)介紹現(xiàn)有威布爾更新函數(shù)的近似式，分析其精度和適用范圍。第2節(jié)提出新的近似式，分析其精度和適用范圍。第3節(jié)通過一個維修政策優(yōu)化的數(shù)例例證提出的近似式的精確性和有用性。第4節(jié)總結(jié)全文。

1 現(xiàn)有近似式及其精度

1.1 更新函數(shù)精確值的計算

計算更新函數(shù)精確值是評價一個近似式精確性的先決條件。令F(t)記非負隨機變量T的分布函數(shù)。本文著重考慮F(t)是威布爾分布的情況，其表達式為:

F(t)=1-exp[-(t/η)β]

(1)

這里，β是形狀參數(shù)，η是尺度參數(shù)。在可靠性和維修應(yīng)用領(lǐng)域，β一般大于1，且很少超過4[12]。平均(μ)、標準偏差(σ)和變異系數(shù)(cv)分別為

μ=ηΓ(1+1/β),

σ=η[Γ(1+2/β)-Γ2(1+1/β)]0.5,cv=σ/μ

(2)

這里Γ(.)是伽瑪函數(shù)。

令F(k)(t)記F(t)的k重卷積(F(1)(t)=F(t))，N(t)記時間區(qū)間(0,t]內(nèi)的更新事件數(shù)。更新函數(shù)(記為M(t))是離散隨機變量N(t)的數(shù)學期望，可表達為下面的級數(shù)形式或積分形式[15,19]:

(3)

對于大多數(shù)分布函數(shù)，包括威布爾分布，上述的卷積和積分沒有解析表達式。因此，開發(fā)更新函數(shù)的近似式成為一個重要的研究方向。最簡單、最著名的兩個近似關(guān)系為:

(4)

前者適合小t的情況，后者適合于大t的情況。

對于一個給定的近似式(記為M∞(t)),一個重要的問題是確定它在時間軸t上的適用范圍。這個問題等價于評價近似式的精度。為此，需要知道一些β和t的組合下的更新函數(shù)的精確值。對于一個給定的β和t的組合，其精確值可通過數(shù)值積分法解式(3)中的積分方程獲得[20,21]。文獻[17]的補充材料(Supplementary Material)顯示了η=1，β=1.0(0.5)4.5 (即從1.0到4.5，步長為0.5)和t=0.05(0.05)3.00(即從0.05到3.00，步長為0.05)條件下的威布爾更新函數(shù)精確值。這樣，近似式Ma(t)的精確性可由下式評價:

ε(t;β)=|1-Ma(t;β)/M(t;β)|

(5)

對于一個給定的β，令εM記上述時間范圍內(nèi)ε(t;β)值中的最大值，用它作為評價近似式精確性的測度。

方程(4)的第一個關(guān)系的適用范圍可定義為(0,τ1)[17]，其中

τ1=sup{t;ε(t;β)<0.01}

(6)

方程(4)的第二個關(guān)系的適用范圍可定義為(τ2,∞ )[17]，其中

τ2=inf{t:ε(t;β)≤0.01}

(7)

這樣，“小的t”意味著t<τ1；“中等大小的t”意味著τ1τ2。

表1 近似式的εM值

1.2 現(xiàn)有近似式及其精度

以下主要介紹分段形式和混合形式的近似式；級數(shù)形式的近似式因缺乏簡單性不作介紹。

最早的分段形式的近似式可追溯到Spearman[10]，其表達式為

Ma(t)=max[F(t),M2(t)]

(8)

(9)

對于威布爾分布，F(xiàn)(t)和M2(t)有一個唯一的交點[12]。因此，式(8)實際上是一個兩重分段模型。表1第2列顯示了式(8)的εM值。從中可以看出，該近似的精度偏低。

Parsa和Jin[13]提出一個可積的三重分段線性近似:

(10)

最近，Jiang[11]提出下面的兩重分段近似式

(11)

其中，ts是分界點，根據(jù)該點處的連續(xù)和平滑條件予以確定

M1(t)=H(t)/[1+αH(t)],α≥0

(12)

H(t)是累積風險函數(shù)，α是ts的一個已知函數(shù)。表1第4列給出了這個近似式的εM值。從中可以看出，這個近似式比前述兩個近似式更簡單、更精確。但是，隨著β的增大，εM值快速增大。

Jiang[12]修改式(12)為

M1(t)-F(t)/[1-αF(t)],α≥0

(13)

分界點ts以類似的方式確定，α是ts的已知函數(shù)。表1第5列給出了該近似式的εM值。從中可以看出，當β較大時，該近似式比由式(11)和(12)所構(gòu)成的近似式更精確。然而，當β較大時，εM值仍然偏大。

Jiang和Chen[17]提出下面的混合形式的近似式:

Ma(t)=w(t)M1(t)+[1-w(t)]M∞(t)

(14)

`其中，由下式確定[9]：

M1(t)=pF(t)+(1-p)H(t)
p=1-exp{-{(β-1)/0.873}0.9269}

(15)

權(quán)函數(shù)是β的函數(shù)：

w(t)=1-Φ(t;μw,σw)

(16)

這里，Φ(t;μw,σw)是參數(shù)為μw和σw的正態(tài)分布函數(shù)，μw=(0.9139+0.2020β)η,σw=(|0.6302β-2.0001|+0.1226)η/6。

2 提出的近似式及其精度

2.1 提出的近似式

提出的近似關(guān)系取式(14)的形式，具有不同的M1(t)；權(quán)函數(shù)w(t)有不同的參數(shù)。

Baker[22]提到，當t>2μ時，式(4)中的M∞(t)一般地是精確的。這意味，M1(t)應(yīng)該恰當?shù)卮_定使其在時間區(qū)間(0,2μ)內(nèi)是精確的。下式可能滿足這個要求：

(17)

為簡單起見，使用一個分布函數(shù)G2(t)近似F(2)(t)；使用另一個分布函數(shù)G3(t)近似F(3)(t)。這樣，M1(t)被定義為

M1(t)=F(t)+G2(t)+G3(t)

(18)

考慮威布爾、伽瑪、正態(tài)和對數(shù)正態(tài)分布作為G2(t)和G3(t)的備選模型。對于一個給定的β和η的組合，使用式(2)計算F(t)的均值μ和方差σ2。則F(k)(t)(k=2,3)的平均為kμ，方差為kσ；Gk(t)的參數(shù)用矩法確定。例如，如果Gk(t)是形狀參數(shù)為uk、尺度參數(shù)為v的伽瑪分布函數(shù)，則

uk=kμ2/σ2,v=σ2/μ

(19)

表2 備選模型的εA值

為了對比備選模型的性能，類似于計算威布爾更新函數(shù)的精確值，使用數(shù)值積分方法計算F(k)(t)的精確值。對于一個給定的備選模型和β值，計算在t/η=0.05(0.05)2.00處的相對誤差εg(t)，εA令記這些相對誤差的平均，用它作為評價備選模型性能的測度。采用相對誤差(而非絕對誤差)的目的是為了強調(diào)備選模型的左尾特征，因為M1(t)適用于中偏小的t值。

表2顯示了備選模型的值εA。它清楚地表明，伽瑪分布是最好的。因此，G2(t)和G3(t)是伽瑪分布，其參數(shù)由式(19)計算?，F(xiàn)在確定權(quán)函數(shù)的參數(shù)：μw和σw。仔細分析M1(t)和M∞(t)之間的交點發(fā)現(xiàn)有以下三種情況：

·對于β=1, 它們相交在原點。

·當β大于且接近于1時，有一個交點(參見圖1，ΔM(t)=M1(t)-M∞(t))。

·當β較大時，M1(t)和M∞(t)之間有多個交點(參見圖2)。用t1記最小的交點，用t2記最大的交點。在這兩點之間，隨著t的增大，和越來越接近。

圖1 ΔM(t)=M1(t)-M∞(t)的圖形(β=1.5)

對于β=1.0(0.5)4.5，表3第2和第3行給出了對應(yīng)的t1/η和t2/η的值。平均參數(shù)μw取t1和t2的加權(quán)和，即

μw=w1t1+w2t2

(20)

為使μw更接近于t2，假定wk(k=1,2)正比于tk。這得wk=tk/(t1+t2)。式(20)成為

(21)

表3 t1/η、t2/η、μw/η和σw/η的值

標準偏差參數(shù)σw取為

σw=(t2-μw)/Φ-1(0.999;0,1)=(t2-μw)/3.0902

(22)

這里Φ-1(.)是正態(tài)分布的逆函數(shù)。當t1=0時，μw=t2，取σw=0.01η。表3最后兩行給出了μw/η和σw/η的值。

2.2 本文近似的精度

表1倒數(shù)第2列給出了本文提出的近似式的εA值；最后一列給出最好的近似式。從中可以看出，在β=1.5和2.0之間，最好的近似式由文獻[11]的近似變?yōu)槲墨I[17]的近似。使用二次多項式插值法推斷兩者在β=1.69處有相同的εA值。類似地，在β=3.5和4.0之間，最好的近似式由文獻[17]的近似變?yōu)楸疚奶岢龅慕?，并且推斷兩者在?3.65處有相同的εA值。因此，我們可以得出以下結(jié)論：

·對于β<1.69，由Jiang[11]所提出的近似式提供最好的精度；

·對于β=1.69-3.65，由Jiang和Chen[17]所提出的近似式有最高的精度；

·對于β>3.65，本文所提出的近似式是最精確的。

·按照以上β值范圍選擇模型，最大相對誤差將小于2%。這個精度滿足一般應(yīng)用的要求。

3 示例

當一個產(chǎn)品含有幾個相同的元件(如軸承)，作為一個基于時間的預(yù)防維修政策，批更換政策(block replacement policy)可能適用于這些元件的預(yù)防更換[23]。在這個政策下，每隔時間周期Tp，元件被預(yù)防地更換；如果元件在此之前發(fā)生失效，則只更換失效的元件。

在時間區(qū)間((k-1)Tp,kTp)(k=1,2,…) 內(nèi)，每個元件的失效事件形成一個更新過程，平均失效數(shù)為更新函數(shù)M(t)。令cp記每個元件的平均預(yù)防更換費用，cf記失效更換費用。通常，cf比cp大得多。在上述時間區(qū)間內(nèi)的期望費用率為[23]

J(Tp)=[cp+cfM(Tp)]/Tp

(23)

如果采用失效更換政策，費用率將是J(∞)=cf/μ。批更換政策的有效性可用維修費用節(jié)省率描述：

(24)

此值越大越好。

設(shè)想元件壽命服從β=4.0和平均壽命為700小時的威布爾分布；cf和cp的值分別為100和250費用單位。使用精確的更新函數(shù)求得的最優(yōu)解顯示在表4第一列。如果采用失效更換政策，費用率將是0.3571。這表明，批更換政策能帶來ηe=21.42%的維修費用節(jié)省。因此，對于本文數(shù)例，批更換政策的有效性十分明顯。

表4 基于不同近似式的最優(yōu)解和相對誤差值

為進一步例證本文所提出的近似式的精確性，表4最后五列顯示了從其它五個近似式獲得的最優(yōu)解和相對誤差值。從中可以看出，只有從文獻[17]的近似式得到的最優(yōu)解的精度是可以接受的。這例證了基于β值選擇近似式的必要性。

表5 本文近似式的參數(shù)

3 總結(jié)和結(jié)論

本文系統(tǒng)地總結(jié)了威布爾更新函數(shù)的近似計算式，針對現(xiàn)有近似式在β較大時精度不高的問題，提出了一個新的近似式。精度分析發(fā)現(xiàn)，當β<1.69時，由Jiang[11]所提出的近似式是最精確的；當β∈[1.69,3.65]時，由Jiang和Chen[17]所提出的近似式是最精確的；當β>3.65時，本文所提出的近似式是最精確的。若按此選擇近似式，最大相對誤差將小于2%。

本文數(shù)例例證了所提出的近似式的精確性，基于β值選擇近似式的必要性，以及維修決策優(yōu)化的有效性。

注意到用伽瑪分布G2(t)和G3(t)近似F(2)(t)和F(3)(t)的精度仍有改進的空間(見表2)，故尋找更好的分布模型供近似F(2)(t)和F(3)(t)是一個值得進一步研究的課題。