馮海林, 郭 甜

(西安電子科技大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 西安 710071)

在可靠性壽命測試中觀測得到的數(shù)據(jù)形式有兩種, 分別是組成系統(tǒng)各元件的壽命和系統(tǒng)的壽命[1]. 對結(jié)構(gòu)確定的系統(tǒng), 實(shí)驗(yàn)中只能觀測到系統(tǒng)的壽命, 即當(dāng)系統(tǒng)工作至失效時(shí), 無法確定哪一個(gè)元件導(dǎo)致系統(tǒng)失效. 另一方面, 對單獨(dú)元件的壽命實(shí)驗(yàn), 將各元件配置在系統(tǒng)中并對該系統(tǒng)進(jìn)行壽命測試, 不僅可節(jié)約實(shí)驗(yàn)時(shí)間, 同時(shí)所測試的壽命數(shù)據(jù)也包含了組成系統(tǒng)各元件之間的相關(guān)信息. 因此, 基于系統(tǒng)壽命數(shù)據(jù)對元件壽命分布的統(tǒng)計(jì)推斷有一定的意義[2].

在關(guān)聯(lián)系統(tǒng)的背景下, Samaniego[3]提出了Signature的概念, 它不僅能給出系統(tǒng)的結(jié)構(gòu)信息, 并且可以將系統(tǒng)壽命的可靠性函數(shù)表示為元件次序壽命的可靠性函數(shù)與Signature向量的混合. 在元件壽命的參數(shù)估計(jì)方面, Bhattacharya等[1]基于系統(tǒng)失效壽命數(shù)據(jù), 討論了經(jīng)驗(yàn)元件分布的參數(shù)估計(jì), 并推斷了元件的重要度； Ng等[4]基于Signature已知的系統(tǒng)壽命數(shù)據(jù), 針對元件壽命為比例失效率(PHR)模型的情況, 給出了元件壽命分布的參數(shù)估計(jì). 在元件壽命分布的非參數(shù)推斷方面, Balakrishnan等[5]使用系統(tǒng)壽命的次序統(tǒng)計(jì)量, 研究了元件壽命分布分位點(diǎn)的精確非參數(shù)置信區(qū)間, 并討論了元件壽命分布分位點(diǎn)的置信下限和雙側(cè)置信限； Volterman等[6]考慮了系統(tǒng)失效數(shù)據(jù)的融合, 將兩個(gè)或兩個(gè)以上具有不同Signature的系統(tǒng)聯(lián)立, 基于大量系統(tǒng)的壽命數(shù)據(jù), 給出了元件壽命的準(zhǔn)確聯(lián)合非參數(shù)推斷. 以上統(tǒng)計(jì)推斷均基于完全壽命實(shí)驗(yàn), 即進(jìn)行到所有的系統(tǒng)完全失效, 為節(jié)省時(shí)間和成本, 可采用刪失壽命數(shù)據(jù)進(jìn)行統(tǒng)計(jì)推斷. Balakrishnan等[7]通過考慮Signature已知的可靠性系統(tǒng)的Ⅱ型右刪失壽命數(shù)據(jù), 研究了當(dāng)元件壽命屬于一般尺度參數(shù)族和位置-尺度參數(shù)族時(shí), 元件壽命分布參數(shù)的最優(yōu)線性無偏估計(jì), 并進(jìn)一步討論了系統(tǒng)失效時(shí)間的最優(yōu)線性無偏預(yù)測； 文獻(xiàn)[8]將其拓展到線性統(tǒng)計(jì)推斷中, 對大量Monte-Carlo(MC)仿真的Ⅱ型右刪失系統(tǒng)失效數(shù)據(jù), 給出了Signature已知時(shí)元件壽命分布的位置-尺度參數(shù)的極大似然估計(jì)和基于回歸的估計(jì), 結(jié)果表明, 這些估計(jì)優(yōu)于前面的最優(yōu)線性無偏估計(jì)； 文獻(xiàn)[9]結(jié)合結(jié)構(gòu)已知的完全系統(tǒng)壽命和刪失系統(tǒng)壽命數(shù)據(jù), 使用隨機(jī)最大期望(EM)算法獲得了相關(guān)參數(shù)的極大似然估計(jì)(MLE), 并證明了該算法對多種元件壽命分布族都適用.

對于一個(gè)實(shí)際運(yùn)行的系統(tǒng), 由于條件所限或記錄延遲等原因, 可觀測到的壽命數(shù)據(jù)更多是Ⅱ型雙刪失樣本數(shù)據(jù), 即前(l-1)個(gè)和后(m-r)個(gè)數(shù)據(jù)刪失, 此時(shí)觀測到的壽命數(shù)據(jù)次序統(tǒng)計(jì)量為

Tl ∶m

這類數(shù)據(jù)在實(shí)際生活中廣泛存在[10]. 因此, 針對Ⅱ型雙刪失壽命數(shù)據(jù)元件壽命分布的統(tǒng)計(jì)推斷更有實(shí)際意義. 本文針對Signature已知的系統(tǒng), 基于Ⅱ型雙刪失系統(tǒng)壽命數(shù)據(jù), 在元件壽命分布為位置-尺度分布族情形下, 給出元件壽命分布的統(tǒng)計(jì)推斷結(jié)果, 并以Weibull壽命分布做驗(yàn)證, 即通過MC仿真實(shí)驗(yàn)?zāi)M產(chǎn)生系統(tǒng)壽命樣本, 分別設(shè)計(jì)N-M算法、 SQP算法和免疫克隆選擇算法實(shí)現(xiàn)點(diǎn)估計(jì)的結(jié)算, 最后給出參數(shù)的漸近置信區(qū)間估計(jì).

1 基于Signature的系統(tǒng)可靠性函數(shù)

(1)

引理1表明, 當(dāng)系統(tǒng)的Signature已知時(shí), 可用系統(tǒng)壽命數(shù)據(jù)實(shí)現(xiàn)對系統(tǒng)元件壽命分布參數(shù)的推斷. 由引理1, 系統(tǒng)壽命T的概率密度可表示為

(2)

特別地, 當(dāng)元件壽命X的概率密度函數(shù)為

(3)

時(shí), 可得系統(tǒng)壽命的概率密度函數(shù)為

(4)

2 基于Signature的系統(tǒng)元件壽命參數(shù)推斷

2.1 元件壽命參數(shù)的似然估計(jì)

首先給出元件壽命分布參數(shù)的點(diǎn)估計(jì). 假設(shè)對m個(gè)系統(tǒng)進(jìn)行壽命測試實(shí)驗(yàn), 每個(gè)系統(tǒng)由獨(dú)立同分布的n元件構(gòu)成, 前(l-1)個(gè)系統(tǒng)壽命數(shù)據(jù)丟失, 觀測得到第l～r個(gè)系統(tǒng)的失效時(shí)間, 即得到一個(gè)Ⅱ型雙刪失樣本

tl ∶m

其中l(wèi),r分別為左、 右刪失數(shù)量,l/n和r/n為樣本雙側(cè)刪失比. 基于上述系統(tǒng)壽命數(shù)據(jù), 令θ=(μ,σ)′為參數(shù)向量, 則θ的似然函數(shù)為

相應(yīng)的對數(shù)似然函數(shù)為

其中c是一個(gè)不依賴于參數(shù)的常數(shù). 將式(6)中系統(tǒng)壽命的PDF,CDF,SF分別由等式(1),(2)替換可得

為了給出參數(shù)θ=(μ,σ)′的極大似然估計(jì), 先對似然函數(shù)式(5)或?qū)?shù)似然函數(shù)式(6)取最大值, 再對其求解. 式(6)的一階、 二階偏導(dǎo)分別為

為計(jì)算系統(tǒng)壽命分布的偏導(dǎo)函數(shù)式(8),(9), 需先計(jì)算下列系統(tǒng)壽命分布的PDF和SF對相應(yīng)位置-尺度參數(shù)的局部偏導(dǎo):

2.2 元件壽命參數(shù)的漸近置信區(qū)間估計(jì)

在點(diǎn)估計(jì)結(jié)果的基礎(chǔ)上, 可給出參數(shù)的漸近置信區(qū)間估計(jì). 由式(9)計(jì)算得到局部Fisher信息矩陣為

(14)

參數(shù)極大似然估計(jì)的方差-協(xié)方差矩陣為

(15)

(16)

(17)

其中,zq為標(biāo)準(zhǔn)正態(tài)分布的第q個(gè)上分位點(diǎn). 由于參數(shù)σ≥0, 故可修正置信區(qū)間為

(18)

2.3 參數(shù)估計(jì)的優(yōu)化算法設(shè)計(jì)

理論上, 參數(shù)的極大似然估計(jì)需要求解下列似然方程組：

(19)

顯然式(19)是超越方程構(gòu)成的似然方程組, 對復(fù)雜的元件壽命分布如Weibull分布, 其直接求解相當(dāng)復(fù)雜甚至不可能. 按照極大似然原理, 似然函數(shù)極大化問題即為下列函數(shù)的極小化問題:

min{-l(θ;tl ∶m,t2∶m,…,tr ∶m)},

(20)

式(20)是一個(gè)復(fù)雜的多維非線性無約束問題. 本文設(shè)計(jì)3種優(yōu)化算法求解該問題.

1) Nelder-Mead算法(簡稱N-M算法), 即沿著使目標(biāo)函數(shù)值下降的方向逐步逼近最優(yōu)解；

2) SQP算法(或稱序列二次規(guī)劃法), 即為了減少計(jì)算量, 將對數(shù)似然函數(shù)的梯度信息加入到優(yōu)化求解中, 用一系列的線性規(guī)劃或二次規(guī)劃逐次逼近原非線性規(guī)劃問題的最優(yōu)解；

3) 免疫克隆選擇算法, 即為避免初始值對迭代結(jié)果的影響, 將抗體親和力對應(yīng)于似然函數(shù)模型(7), 抗體即為參數(shù)估計(jì)值. 設(shè)定一定的抗體群規(guī)模, 抗體初始化時(shí), 隨機(jī)產(chǎn)生一代抗體, 計(jì)算親和力大小. 在一代進(jìn)化中, 根據(jù)親合力大小, 對抗體進(jìn)行選擇、 克隆和變異操作, 得到一組子抗體群, 并對其進(jìn)行再選擇, 篩選出新一代抗體, 設(shè)定一定迭代次數(shù)后, 將最后一代抗體群中的最優(yōu)抗體作為最優(yōu)解輸出.

3 MC仿真實(shí)驗(yàn)

仿真的軟件環(huán)境為Windows7, 64位操作系統(tǒng)和MATLAB R2014a； 硬件環(huán)境為Intel Core i5-3230處理器, 4 GB內(nèi)存. 本文針對Navarro等[12]提出的3種不同Signature系統(tǒng), 進(jìn)行MC(Monte Carlo)仿真研究, 仿真實(shí)驗(yàn)設(shè)置中, 根據(jù)m的大小, 將仿真分為小樣本和大樣本兩種情況. 其中： 系統(tǒng)1是一個(gè)四元件串并聯(lián)Ⅲ系統(tǒng), 系統(tǒng)Signature為s1=(1/4,1/4,1/2,0)； 系統(tǒng)2是一個(gè)四元件混合并聯(lián)Ⅰ系統(tǒng), 系統(tǒng)Signature為s2=(0,1/2,1/4,1/4)； 系統(tǒng)3是一個(gè)三元件并串聯(lián)系統(tǒng), 系統(tǒng)Signature為s3=(0,2/3,1/3).

仿真實(shí)驗(yàn)實(shí)例選取系統(tǒng)的元件壽命W為Weibull分布

(21)

其中:β>0是尺度參數(shù);α>0是形狀參數(shù). 則隨機(jī)變量X=lnW服從最小極值(SEV)分布, CDF和PDF分別為

(22)

(23)

(24)

(25)

仿真過程如下：

1) 產(chǎn)生系統(tǒng)Signature為s=(s1,s2,…,sn)、 元件壽命為Weibull分布的m組系統(tǒng)壽命數(shù)據(jù)T1∶m

① 由一個(gè)離散n點(diǎn)分布、 以概率P(V=v)=sv(v=1,2,…,n)產(chǎn)生一個(gè)隨機(jī)數(shù)v；

② 由次序統(tǒng)計(jì)量的性質(zhì), 從正態(tài)分布產(chǎn)生第v個(gè)次序統(tǒng)計(jì)量U(v), 其中U(v)服從參數(shù)為v和m+1-v的Beta分布；

2) 設(shè)定左右刪失量l和r(1

3) 分別用N-M算法、SQP算法、 免疫克隆選擇算法求解式(20), 得到參數(shù)的點(diǎn)估計(jì)與區(qū)間估計(jì).

4) 設(shè)定不同大小的樣本容量和左右刪失數(shù), 重復(fù)過程1)～3).

5) 重復(fù)步驟1)～4)共10 000次, 并計(jì)算點(diǎn)估計(jì)值的偏差(Bias)和均方誤差(MSE), 以及各區(qū)間估計(jì)的95%置信區(qū)間的覆蓋率(CR)和平均覆蓋寬度(ACI).

3.1 小樣本仿真

設(shè)壽命實(shí)驗(yàn)中有m=10個(gè)樣本系統(tǒng), 元件壽命的Weibull分布中參數(shù), β=3, α=2, 則對數(shù)轉(zhuǎn)換后對應(yīng)的位置、 尺度參數(shù)分別為μ=ln3, σ=1/2. 設(shè)壽命樣本的左刪失率p=10%, 即l=2, 令右刪失率q分別為10%,20%,…,70%, 即r=m(1-q)=9,8,…,2, 則可得3種不同Signature下參數(shù)μ和σ的點(diǎn)估計(jì)與區(qū)間估計(jì)的仿真結(jié)果分別如圖1和圖2所示.

圖1 l=2時(shí)元件壽命位置參數(shù)μ的點(diǎn)估計(jì)與區(qū)間估計(jì)結(jié)果 Fig.1 Point estimation and interval estimation results of component lifetime location parameter μ when l=2

由圖1和圖2可見:N-M算法區(qū)間估計(jì)的平均寬度可穩(wěn)定在較小的值, 但對區(qū)間估計(jì)的覆蓋率較低, 在對σ的估計(jì)中, 區(qū)間估計(jì)的覆蓋率明顯低于其他兩種算法；SQP方法點(diǎn)估計(jì)的均方誤差最小, 區(qū)間估計(jì)覆蓋率大, 但隨著刪失率的增大, 還能保證較高的覆蓋率, 同時(shí)區(qū)間的平均寬度也較大； 免疫克隆選擇算法在對μ的點(diǎn)估計(jì)中, 偏差均為負(fù)且可以收斂到0, 區(qū)間刪失率較低時(shí)能保證較高的覆蓋率, 但隨著刪失率增大, 覆蓋率下降很快.

在算法迭代過程中,N-M算法對初值的選擇較敏感, 數(shù)值迭代過程中需要選擇適合的初始值； 免疫選擇算法隨機(jī)生成初始參數(shù)值, 通過數(shù)值迭代選擇, 逐步尋找到最優(yōu)的估計(jì)值, 對初始值不敏感, 魯棒性強(qiáng)；SQP算法給出的位置參數(shù)和尺度參數(shù)的點(diǎn)估計(jì)和區(qū)間估計(jì)效果最好, 精度更高.

圖2 l=2時(shí)元件壽命位置參數(shù)σ的點(diǎn)估計(jì)與區(qū)間估計(jì)結(jié)果Fig.2 Point estimation and interval estimation results of component lifetime location parameter σ when l=2

3.2 大樣本仿真

設(shè)壽命實(shí)驗(yàn)中有m=50和m=100個(gè)系統(tǒng), 本文列舉了左右刪失數(shù)不同7種情況下3種不同算法對不同Signature的仿真結(jié)果, 點(diǎn)估計(jì)與區(qū)間估計(jì)結(jié)果分別列于表1～表4.

表1 大樣本下N-M算法、 SQP算法及免疫克隆選擇算法對參數(shù)μ的點(diǎn)估計(jì)結(jié)果

表2 大樣本下N-M算法、 SQP算法及免疫克隆選擇算法對參數(shù)σ的點(diǎn)估計(jì)結(jié)果

表3 大樣本下N-M算法、 SQP算法及免疫克隆選擇算法對參數(shù)μ的區(qū)間估計(jì)結(jié)果

表4 大樣本下N-M算法、 SQP算法及免疫克隆選擇算法對參數(shù)σ的區(qū)間估計(jì)結(jié)果

表5 10 000次仿真下3種算法的平均運(yùn)行時(shí)間(s)

4 實(shí)例分析

表6 Signature分別為s1,s2,s3, 元件壽命服從Weibull (3,2)分布的系統(tǒng)壽命結(jié)果

設(shè)左、 右刪失數(shù)分別為l=2和r=9, 通過3種不同優(yōu)化算法得到元件壽命分布參數(shù)的點(diǎn)估計(jì)與區(qū)間估計(jì)結(jié)果列于表7. 由表7可見： 由SQP算法得到的點(diǎn)估計(jì)方差最小, 進(jìn)一步得到的置信區(qū)間估計(jì)最精確; 免疫優(yōu)化算法不用設(shè)定初始值, 對不同壽命數(shù)據(jù)下參數(shù)的優(yōu)化結(jié)果值差異不大. 實(shí)驗(yàn)結(jié)果驗(yàn)證了大樣本和小樣本情形下的參數(shù)估計(jì)結(jié)果.

表7 對表6中系統(tǒng)壽命數(shù)據(jù)的估計(jì)結(jié)果

綜上, 本文在Ⅱ型雙刪失系統(tǒng)壽命數(shù)據(jù)下, 研究了基于系統(tǒng)Signature的元件壽命分布參數(shù)的統(tǒng)計(jì)推斷問題, 使用極大似然估計(jì)及其漸近理論給出了參數(shù)的點(diǎn)估計(jì)和漸近置信區(qū)間估計(jì), 并設(shè)計(jì)了優(yōu)化計(jì)算似然函數(shù)的N-M算法、 SQP算法和免疫克隆選擇算法. 以元件壽命為Weibull分布為例, 分別給出了小樣本和大樣本實(shí)驗(yàn)數(shù)據(jù)下的估計(jì)結(jié)果, 表明SQP算法給出的參數(shù)估計(jì)效果最好. 隨著樣本量的增大, 基于Ⅱ型雙刪失樣本的似然函數(shù)偏導(dǎo)表達(dá)式更復(fù)雜, 雖然估計(jì)精度和效果會隨著樣本數(shù)據(jù)的增大而提高, 但也會使點(diǎn)估計(jì)的方差和協(xié)方差計(jì)算量更大, 進(jìn)而使區(qū)間估計(jì)的計(jì)算變得繁瑣, 針對此問題, 本文設(shè)計(jì)的3種優(yōu)化算法均能有效解決, 并得到了較準(zhǔn)確的點(diǎn)估計(jì)和區(qū)間估計(jì)結(jié)果.