應(yīng)用MS Excel求解三參數(shù)威布爾分布函數(shù)的參數(shù)估計(jì)

劉子娟， 鄭學(xué)斌， 郭小軍

（1.中國(guó)湖南航空動(dòng)力機(jī)械研究所，湖南 株洲412002；2.長(zhǎng)沙升華微電子材料有限公司，長(zhǎng)沙410000）

0 引 言

在機(jī)械零件的疲勞試驗(yàn)數(shù)據(jù)統(tǒng)計(jì)分布中，廣泛應(yīng)用正態(tài)分布或?qū)?shù)正態(tài)分布[1]。但是，當(dāng)機(jī)械零件失效概率很小時(shí)，用正態(tài)分布評(píng)估的疲勞壽命趨近于0[2]。威布爾分布是瑞典物理學(xué)家Waloddi Weibull在分析材料強(qiáng)度時(shí)，根據(jù)實(shí)際經(jīng)驗(yàn)推導(dǎo)出來(lái)的分布形式[3]，與常用的正態(tài)分布或是對(duì)數(shù)正態(tài)分布相比，三參數(shù)威布爾分布往往能更準(zhǔn)確地描述零件疲勞強(qiáng)度或疲勞壽命的概率分布，這是因?yàn)槿齾?shù)威布爾分布中有個(gè)位置參數(shù)，其表征機(jī)械零件的最小壽命或最低疲勞強(qiáng)度，這更能反映機(jī)械零件疲勞特性的實(shí)際。因此，三參數(shù)威布爾分布擬合精度更高，在機(jī)械零件強(qiáng)度及疲勞壽命評(píng)價(jià)中得到了越來(lái)越廣泛的應(yīng)用[4-5]。

能用于三參數(shù)威布爾分布參數(shù)估計(jì)的方法大多計(jì)算復(fù)雜[6]，需要通過(guò)Matlab或其他計(jì)算機(jī)語(yǔ)言編程來(lái)計(jì)算，對(duì)于不掌握這些編程知識(shí)的人員來(lái)說(shuō)應(yīng)用難度很大。而MS Excel是大多數(shù)普通人員常用到的軟件，本文運(yùn)用Excel自帶的強(qiáng)大數(shù)學(xué)運(yùn)算、統(tǒng)計(jì)分析程序，對(duì)三參數(shù)威布爾分布進(jìn)行參數(shù)估計(jì)，以一組數(shù)據(jù)為例來(lái)說(shuō)明如何使用Excel中的數(shù)學(xué)運(yùn)算函數(shù)、規(guī)劃求解功能，以實(shí)現(xiàn)快速、高效地求解三參數(shù)威布爾分布的參數(shù)估計(jì)問(wèn)題。

1 三參數(shù)威布爾分布的壽命分布函數(shù)

分布函數(shù)為

式中：N為零件壽命，N≥γ；m為形狀參數(shù)，m＞0；η為尺度參數(shù)，η＞0；γ為位置參數(shù)，在機(jī)械零件壽命分析中，γ作為零件的最小壽命值，表示產(chǎn)品的使用壽命在達(dá)到γ值以前不會(huì)失效。因此對(duì)于實(shí)際零件產(chǎn)品來(lái)說(shuō)，γ≥0（特別地當(dāng)γ=0時(shí)，三參數(shù)威布爾分布函數(shù)退化成二參數(shù)威布爾分布函數(shù)）。

對(duì)式（1）兩邊取2次自然對(duì)數(shù)處理后得到

設(shè)零件壽命估計(jì)的可靠度為R，根據(jù)失效率F與可靠度R（即存活率）的關(guān)系

對(duì)式（5）作適當(dāng)變換得

然后對(duì)式（6）兩邊取自然對(duì)數(shù)后整理得

式中，NR為三參數(shù)威布爾分布在可靠度下R的壽命估計(jì)。

2 位置參數(shù)γ的求解方法

先給一個(gè)初始值γ0（γ0≤Nmin，Nmin為試驗(yàn)得到的零件壽命數(shù)據(jù)組中最小壽命），利用最小二乘法，可求的m值、B值。

式（8）、式（9）中，n為樣本容量。

式中：

當(dāng)線性相關(guān)符合度r值的絕對(duì)值趨近1時(shí)，說(shuō)明兩個(gè)變量之間線性相關(guān)良好。為了使r值的絕對(duì)值越接近于1，需要優(yōu)化γ0值。

3 經(jīng)驗(yàn)分布的累計(jì)失效概率選取

根據(jù)經(jīng)驗(yàn)，設(shè)有n個(gè)產(chǎn)品進(jìn)行壽命試驗(yàn)數(shù)據(jù)，按其失效壽命值Ni（i=1、2、3……n）由小到大排序?yàn)镹1＜N2＜N3＜……＜Nn，對(duì)應(yīng)的累積失效概率（經(jīng)驗(yàn)分布函數(shù)）為F（N1）＜F（N2）＜……＜F（Nn）。其中第i個(gè)失效產(chǎn)品的累積失效概率F（Ni）可用中位秩算法求得

或用平均秩算法：

4 運(yùn)用MS Excel進(jìn)行三參數(shù)威布爾分布的參數(shù)估計(jì)

MS Excel具有強(qiáng)大的統(tǒng)計(jì)和計(jì)算功能，其內(nèi)含的“規(guī)劃求解”功能更是求解最優(yōu)化問(wèn)題的便捷計(jì)算工具。因此，在MS Excel中編制好相應(yīng)的計(jì)算公式后，可以很容易求解出三參數(shù)威布爾分布函數(shù)中位置參數(shù)γ、形狀參數(shù)m和尺度參數(shù)η的最優(yōu)值，并且可依據(jù)得到的3個(gè)參數(shù)值，求出指定某個(gè)可靠度R值下對(duì)應(yīng)的零件壽命值。

現(xiàn)就以下面一個(gè)例子來(lái)說(shuō)明，如何使用MS Excel表格來(lái)對(duì)機(jī)械零件進(jìn)行三參數(shù)威布爾分布的參數(shù)求解和壽命估計(jì)。

設(shè)某零件經(jīng)疲勞試驗(yàn)，測(cè)得其疲勞壽命循環(huán)次數(shù)分別為：29 613、35 862、42 267、46 500、51 233、54 600、60 970、66 820、71 608、75 088、79 604、82 200、87 368、92 570、100 790，共15個(gè)值（單位：次）。然后按圖1的樣式在MS Excel中設(shè)置表格。

1）在單元格A2～A16中由小到大升序方式輸入上述試驗(yàn)測(cè)得的疲勞壽命循環(huán)次數(shù)值。

2）在單元格G18中指定一個(gè)初始值γ0，此處取29 000（注意，γ0≤Nmin）。

3）在單元格B2～B16中依次輸入序號(hào)1～15。

5）在單元格D2～D16中，計(jì)算Ni-γ0值。在D2單元格中輸入公式“=A2-$G$18”，同樣用下拉填充方式填充D3～D16單元格，得到D3～D16單元格中的計(jì)算值。

6）在單元格E2～E16中，計(jì)算1-F（Ni）值。在E2單元格中輸入公式“=1-C2”，用下拉填充方式填充E3～E16單元格，在對(duì)應(yīng)單元格中得到計(jì)算值。

7）在單元格F2～F16中，計(jì)算Xi=ln（Ni-γ0）值。在F2單元格中輸入公式“=ln（D2）”，用下拉填充方式填充F3～F16單元格，得到對(duì)應(yīng)單元格的計(jì)算值。

這樣，我們就完成了計(jì)算表格制作的第1步。

接著，根據(jù)所選取的初始值γ0值，對(duì)圖1表格中的第F列（X值）和第G列（Y）進(jìn)行線性規(guī)劃擬合計(jì)算,然后求出線性擬合方程的斜率m值、截距B值和線性相關(guān)系數(shù)r值。

在此調(diào)用MS Excel自帶的函數(shù)命令進(jìn)行計(jì)算以上所述的m值、B值、r值。

在單元格G19中輸入“=SLOPE（G2:G16,F2:F16）”，調(diào)用MS Excel中的求斜率函數(shù)命令，求式（3）所指的線性擬合方程中的斜率m值；在單元格G20中輸入“=INTERCEPT（G2:G16,F2:F16）”，調(diào)用MS Excel中的求截距函數(shù)命令，求式（3）所指的線性擬合方程中的-B值（即得截距B值）；在單元格G21中輸入“=CORREL（G2:G16,F2:F16）”，調(diào)用MS Excel中的求相關(guān)系數(shù)函數(shù)命令，用以反映［Xi,Yi］線性擬合的相關(guān)程度，即線性相關(guān)性符合度r值；在單元格G22中輸入“=G19”，即為形狀參數(shù)m值；在單元格G23中輸入“=EXP（-G20/G22）”，即為尺度參數(shù)η值；在單元格G24中輸入“=INT（G18）”，即取位置參數(shù)γ值為整數(shù)。設(shè)置好上述初始條件（已經(jīng)對(duì)位置參數(shù)賦予初始值），并輸入相關(guān)計(jì)算公式后，得到了一組初始數(shù)據(jù)，如圖1所示。

圖1 Excel計(jì)算表設(shè)計(jì)

最后，使用MS Excel 2003中的“規(guī)劃求解”功能估計(jì)位置參數(shù)的最優(yōu)值。先在MS Excel中加載“規(guī)劃求解”宏。加載完成后，選擇“工具→規(guī)劃求解”功能，打開(kāi)規(guī)劃求解參數(shù)對(duì)話框，如圖2所示。

先將圖2中的“設(shè)置目標(biāo)單元格”設(shè)為單元格“$G$21”，且目標(biāo)值“等于最大值”，即要求該單元格所代表的相關(guān)系數(shù)r值為最大值；再將圖2中的“可變單元格”選為單元格“$G$18”，即表示位置參數(shù)γ0可發(fā)生變化。接著設(shè)置約束條件，由于零件的疲勞壽命有0≤γ≤N1（N1即為收集的零件壽命值中的最小值），點(diǎn)擊圖2中“約束”條件偏右的“添加”按鈕，添加約束條件：$G$18＜=$A$2 和$G$18＞=0。

為了保證精度，點(diǎn)擊規(guī)劃求解參數(shù)對(duì)話框右側(cè)的“選項(xiàng)”按鈕，彈出對(duì)話框，如圖3所示，并按照?qǐng)D3所示進(jìn)行參數(shù)調(diào)整，然后點(diǎn)擊“確定”，返回到圖2所示的對(duì)話框。點(diǎn)擊圖2中的“求解”按鈕，MS Excel將自動(dòng)規(guī)劃求解計(jì)算出最終結(jié)果，求解結(jié)束后彈出的對(duì)話框如圖4所示，點(diǎn)擊圖4中的“保存規(guī)劃求解結(jié)果”后“確定”。

最終得到規(guī)劃求解后的新結(jié)果見(jiàn)圖5所示，此時(shí)顯示相關(guān)參數(shù)r值的最大值是0.996 240 275；位置參數(shù)值γ =10 292、形狀參數(shù)值m=2.535 947 749、尺度參 數(shù) 值 η =62 220.315 37。其他參數(shù)的新結(jié)果可參看圖5所示。

圖2 規(guī)劃求解參數(shù)對(duì)話框

圖3 規(guī)劃求解“選項(xiàng)”對(duì)話框

另外在單元格G25中輸入指定可靠度R值，如輸入0.9。在單元格G26中輸入“=G24+G23*LN（1/G25）^（1/G22）”，即為式（7）所指的指定可靠度下的壽命值NR，此處NR=35 909，即經(jīng)求出三參數(shù)威布爾分布函數(shù)的3個(gè)最優(yōu)參數(shù)值后，零件的壽命服從該三參數(shù)值的威布爾分布，且當(dāng)可靠度R值等于0.9的情況下，零件的疲勞壽命循環(huán)次數(shù)為35 909次。如圖5所示。

圖4 規(guī)劃求解結(jié)果

圖5 規(guī)劃求解后的新數(shù)據(jù)結(jié)果

從上述運(yùn)用MS Excel進(jìn)行三參數(shù)威布爾分布的參數(shù)估計(jì)過(guò)程及形成的計(jì)算結(jié)果，可以看出，通過(guò)Excel制作相應(yīng)的表格后，能快速、高效地求出威布爾分布的3個(gè)參數(shù)，并進(jìn)一步評(píng)估出指定可靠度條件下零件的疲勞壽命值，而無(wú)需通過(guò)復(fù)雜的編程來(lái)求解。

5 結(jié) 語(yǔ)

1）本文根據(jù)線性相關(guān)系數(shù)最大化法則，利用MS Excel進(jìn)行數(shù)據(jù)處理并求解三參數(shù)威布爾分布的3個(gè)參數(shù)，MS Excel表格制作簡(jiǎn)便，計(jì)算操作方便快捷，容易上手操作，可以在工程中得到推廣和應(yīng)用。

2）在使用MS Excel的“規(guī)劃求解”進(jìn)行參數(shù)估計(jì)時(shí)，要注意初始值的選擇和求解精度的設(shè)置。初始值應(yīng)盡量選擇與最小失效壽命（N1值）接近，否則在進(jìn)行規(guī)劃求解后很可能會(huì)得不到滿意的最優(yōu)值。

3）由于累計(jì)失效概率有中位秩、平均秩等計(jì)算方式，采用不同的累積失效概率算法，三參數(shù)威布爾分布的參數(shù)估計(jì)和可靠性壽命估計(jì)的計(jì)算結(jié)果將稍有不同。