基于混合威布爾分布的水稻插秧機的可靠性分析及剩余壽命預測

文昌俊, 陳洋洋*, 何永豪, 陳凡

(1.湖北工業(yè)大學機械工程學院, 武漢 430068; 2.湖北省現(xiàn)代制造質量工程重點實驗室, 武漢 430068)

水稻插秧機作為水稻插秧移栽的主要機械設備,在水稻秧苗移栽期間,由于工作環(huán)境差、沖擊、疲勞磨損等原因,勢必會引起機械性能退化,并隨著運行時間的累積最終導致機器失效。因此,正確計算出水稻插秧機的可靠性指標是現(xiàn)代化農業(yè)耕作的必然要求,也是水稻插秧機預測性維護所需的重要參考依據(jù)。此外,剩余壽命預測對水稻插秧機的維護而言,也具有重要的意義,通過對水稻插秧機進行剩余壽命預測,根據(jù)預測結果對水稻插秧機進行預測性維護。

在可靠性工程領域,威布爾分布是最常用的概率分布模型[1],通過形狀參數(shù)的改變,威布爾模型可以描述失效率上升和失效率下降的產品壽命[2]。蔡文斌等[3]利用三參數(shù)威布爾分布,建立抽油桿疲勞壽命預測模型,基于三參數(shù)威布爾分布模型和正態(tài)分布Basquin模型對HL和HY兩種抽油桿進行疲勞壽命預測,結果表明三參數(shù)威布爾分布模型比正態(tài)分布模型有更高的擬合精度。王誼等[4]針對當前國網用電信息采集系統(tǒng)中電能表可靠性評價不夠準確的問題,利用威布爾分布模型,計算出了電能表的不可靠度,實現(xiàn)了威布爾分布理論在電能表評價中的實際應用。喬宏霞等[5]為改善西部鹽湖地區(qū)混凝土耐久性能,采用威布爾函數(shù)建立CaCO3改性混凝土的耐久性退化模型,通過威布爾分布模型計算得出概率密度函數(shù)和可靠度函數(shù),確定出納米CaCO3改性混凝土在硫酸鹽境下的最長使用壽命。繆吉昌等[6]基于嬰兒培養(yǎng)箱風機的失效數(shù)據(jù)建立威布爾分布模型,隨后利用極大似然法估計出模型參數(shù),最后在其模型基礎上對風機進行剩余壽命預測。然而上述文獻都是基于單一威布爾分布模型進行計算的,但是水稻插秧機通常包含多個零部件,結構復雜,在水稻插秧期間,受多種應力作用,多種故障機制并存,這種情況下,用單一的威布爾分布模型對其壽命進行擬合會出現(xiàn)較大誤差,此外文獻也沒提及在混合威布爾的基礎上對設備進行剩余壽命的預測。

為了解決上述問題,在PZ型水稻插秧機故障數(shù)據(jù)的基礎上,現(xiàn)建立兩參數(shù)混合威布爾分布模型,其次,采用改進粒子群算法與非線性最小二乘法相結合對其模型的參數(shù)進行求解,利用Kolmogorov-Smirnov test以及對比單一威布爾模型和混合威布爾模型與失效數(shù)據(jù)的擬合效果證明該模型的合理性,然后通過該模型得到水稻插秧機可靠性指標和預防性維修周期,最后在混合威布爾模型的基礎上對水稻插秧機進行剩余壽命預測。

1 混合威布爾分布模型

由于混合威布爾模型包含更多參數(shù),所以能夠針對較為復雜情況進行數(shù)據(jù)擬合[7]。假設一組數(shù)據(jù)樣本n服從混合威布爾分布,對應各子分布均服從相同分布,但參數(shù)不同,各子分布的權重為w1,w2,…,wm,權重之和為1,各子分布的概率密度函數(shù)分別為f1(t),f2(t),…,fm(t),混合威布爾分布的故障概率密度函數(shù)f(t)、可靠度函數(shù)R(t)和故障概率分布函數(shù)F(t)分別為

(3)

式中:m為混合威布爾分布模型的重數(shù);wi為第i個子分布權重;βi為第i個威布爾分布的形狀參數(shù);ηi為第i個威布爾分布的尺度參數(shù);γi為第i個子分布的位置參數(shù)。由于水稻插秧機從零時刻開始使用便有可能出現(xiàn)故障,所以位置參數(shù)γ=0,此時為兩參數(shù)混合威布爾分布。

以PZ型水稻插秧機歷史維修數(shù)據(jù)為例進行試驗,共收集到282組故障數(shù)據(jù),將282組故障間隔時間按照從小到大順序進行排序,如表1所示。

表1 水稻插秧機故障數(shù)據(jù)

將故障時間t與函數(shù)值F(t)之間的關系通過式(4)轉化為x與y之間的關系,將各點繪制在x-y坐標系中,得到數(shù)據(jù)點分布圖,如圖1所示。

圖1 失效概率分布變換后的曲線Fig.1 The transformed curve of failure probability distribution

(4)

由圖1可知,該曲線存在明顯轉折,轉折前后有明顯不同的斜率,如果用一條直線去擬合圖中的數(shù)據(jù)點,誤差將會很大,因此不能用單一威布爾分布去擬合這些故障數(shù)據(jù),假設數(shù)據(jù)樣本服從兩重兩參數(shù)混合威布爾分布。

2 參數(shù)估計及擬合優(yōu)度檢驗

2.1 非線性最小二乘法

混合威布爾分布模型較復雜,無法將t和F(t)變換成線性關系,所以不能用線性回歸的方法去估計未知參數(shù)[8]。非線性最小二乘法通過最小化殘差平方和的方法估計參數(shù),在數(shù)據(jù)的擬合分析方面有著廣泛應用。利用非線性最小二乘法,建立兩重兩參數(shù)混合威布爾分布的非線性最小二乘優(yōu)化模型。將容量為282的水稻插秧機故障數(shù)據(jù)樣本按從小到大的順序進行排列(t1,t2,…,t282),采用中位秩公式[式(5)]計算可靠度的觀測值,即

(5)

(6)

以殘差平方總和最小為目標,構建混合威布爾分布參數(shù)估計數(shù)學模型為

(7)

式(7)中:g(θ)為混合威布爾分布參數(shù)估計模型的目標函數(shù);βupi、βdowni分別為第i個子分布形狀參數(shù)的最大值和最小值;ηupi、ηdowni分別為第i個子分布尺度參數(shù)的最大值和最小值。

式(7)為超越方程組,采用常規(guī)方法求解存在諸多困難,最佳途徑是使用優(yōu)化算法,找到一組最優(yōu)的θ使對應目標函數(shù)g(θ)的值最小,此時的θ作為模型中未知參數(shù)的估計值。利用改進粒子群算法對式(7)的模型進行求解,將圖解法的估計值作為粒子群算法的搜索空間,圖解法粗略估計子分布權重w1在[0.3,0.4],尺度參數(shù)η1和η2的范圍分別在[50,100]、[200,300],子分布形狀參數(shù)的初始范圍并不影響計算過程,先假定子分布形狀參數(shù)的范圍在[0,10],若計算出的子分布形狀參的值在取值區(qū)間邊界,則放寬該邊界范圍并重新計算。

2.2 改進粒子群算法

粒子群算法是一種基于鳥類群體覓食而提出的一種啟發(fā)式全局優(yōu)化算法[9]。粒子通過個體最優(yōu)解pi和全局最優(yōu)解pg不斷迭代調整自身的位置和速度,從而不斷更新候選解,在一個K維的目標搜索空間中,有N個粒子組成一個群體,粒子通過式(8)和式(9)進行更新迭代。

(9)

式中:vit和xit分別為第i個粒子第t次迭代時的速度矢量和位置矢量;w為慣性權重;c1和c2為學習因子;r1和r2是介于[0,1]的隨機數(shù)。

在標準的粒子群優(yōu)化算法中,慣性權重w、自身認知c1以及社會認知c2這三者在收斂效果中起著關鍵作用[10]。慣性權重w代表上一代粒子的速度對當前粒子速度的影響程度,在搜索初期階段,可以將慣性權重w設置地較大些,保證粒子能夠以較快的速度搜索到較好的區(qū)域,搜索后期,較小的慣性權重可以增加算法的局部搜索能力,保證粒子能在最優(yōu)解區(qū)域做更精確的搜索。采用線性遞減的慣性權重,其表達式為

(10)

式(10)中:wmax、wmin分別為w的最大值、最小值;M為最大迭代次數(shù)。

在算法初期階段,自身認知部分應該占優(yōu)勢,采取較大的自身認知c1,較小的社會認知c2,使得算法初期在整個空間進行搜索,保證算法的全局收斂。隨著迭代次數(shù)的不斷增加,自身認知部分c1應逐漸減小,同時,社會認知部分c2應該逐漸擴大,使粒子更關注于社會信息,保持快速收斂。設計的自適應學習因子表達式為

(11)

式(11)中:cmax、cmin分別為c的最大值、最小值。

將遺傳算法的選擇、交叉和變異操作引入粒子群的優(yōu)化算法中,得到改進粒子群算法,該算法具有遺傳算法全局收斂和粒子群算法收斂速度快的特點,改進粒子群算法的主要步驟如下。

步驟1初始化參數(shù)。包括種群個數(shù)N,最大迭代次數(shù)M,粒子速度和位置的界限,遺傳算法和粒子群算法運行時所需參數(shù)等。

步驟2計算每個粒子的適應度值。將每個粒子的當前位置作為個體最優(yōu)位置pi,當前粒子中適應度最小的粒子位置作為全局最優(yōu)位置pg。

步驟3通過式(8)和式(9)更新粒子的速度和位置,然后檢驗其速度和位置是否超過邊界范圍,并計算更新之后每個粒子的適應度值。

步驟4計算平均適應度值,將平均適應度值的一半作為界限,小于此界限的粒子劃分到子種群P1里面,其余粒子劃分到P2里面。

步驟5子種群P1的粒子適應度值較小,直接復制進入下一代,子種群P2的粒子總維數(shù)的3/5分別與個體極值pi和全局極值pg進行交叉、變異操作,子種群P2通過選擇交叉、變異之后得到的新子種群P3,計算P3中各個粒子的適應度值。將P2和P3粒子的適應度值進行比較,取適應度值小的一半進入到下一代。

步驟6更新個體極值和全局極值。計算新子代的適應度值,與個體極值和全局極值進行比較,如果新粒子的適應度值比個體極值pi的適應度值小,則更新pi,將當前新粒子的位置賦值給pi。如果新粒子的適應度值比群體全局極值pg的適應度值小。則更新pg,將當前新粒子的位置賦值給pg。

步驟7不斷重復步驟3～步驟6,直到目標函數(shù)達到收斂精度或者達到最大迭代次數(shù)。

2.3 參數(shù)估計

在表1中給出的故障數(shù)據(jù)的基礎上,利用改進粒子群算法對式(7)進行求解,改進粒子群算法基本參數(shù):粒子群數(shù)目N=30,迭代次數(shù)M=50,cmax和cmin分別為2.0與0.5,wmax和wmin分別為0.9與0.4。粒子群算法的基本參數(shù):種群和迭代次數(shù)和改進粒子群算法的保持一致,c1和c2為1.494 45,w為0.6,粒子群算法和改進粒子群算法的迭代過程如圖2所示,參數(shù)估計結果如表2所示。

圖2 粒子群算法和改進粒子群算法的優(yōu)化過程Fig.2 Optimization process of particle swarm optimization and improved particle swarm optimization

表2 粒子群算法和混合粒子群算法的優(yōu)化結果

在粒子群迭代過程中加入遺傳算法的交叉、變異,從而增強了種群的多樣性,避免算法陷入局部最優(yōu)。在收斂速度上,從圖2中可發(fā)現(xiàn)粒子群算法在迭代至45代左右達到最優(yōu)值,而改進粒子群算法在迭代至20代左右達到最優(yōu)值,表明改進粒子群算法比粒子群算法收斂速度更快,在收斂精度上,由表2可知,粒子群算法的收斂精度為0.015 2,而改進粒子群算法的收斂精度為0.014 5,表明改進粒子群算法具有更高的收斂精度,將改進PSO算法求出的參數(shù)值:w1=0.326 2,η1=72.821 6,m1=0.889 9,w2=0.673 8,η2=228.839 6,m2=2.050 0作為兩參數(shù)混合威布爾模型的參數(shù)值。

2.4 擬合優(yōu)度檢驗

為了驗證所選分布模型的合理性,需要進行擬合優(yōu)度檢驗。常用的模型檢驗方法[11]有χ2檢驗和K-S檢驗。相較于χ2檢驗,K-S檢驗更加精細[12],應用更為廣泛。因此采用K-S檢驗法來檢驗模型的擬合優(yōu)度。設總體分布為F(t),F0(t)為某個已知的連續(xù)型分布函數(shù),假設檢驗問題H0:F(t)=F0(t)。

現(xiàn)假設:水稻插秧機的故障數(shù)據(jù)服從兩重兩參數(shù)混合威布爾分布。

(12)

代入水稻插秧機的282組故障時間數(shù)據(jù),詳細數(shù)據(jù)如表3所示。

表3 插秧機壽命分布檢驗

由表3結果可知,柯爾莫哥洛夫的檢驗統(tǒng)計量Dn=max{|F0(t)-i/n|,|F0(t)-(i-1)/n|}=0.022 8,當故障數(shù)據(jù)樣本量n=282時,取顯著水平a=0.05,查K-S臨界值表可得d(285,0.05)=0.080 4,d(280,0.05)=0.081 2,通關比較可知Dn=0.022 8

3 可靠性分析及剩余壽命預測

3.1 可靠性分析

可靠性分析是進行剩余壽命預測的基礎[13]。上述結果表明,水稻插秧機故障數(shù)據(jù)服從兩參數(shù)混合威布爾分布,將改進PSO算法求出的模型參數(shù)代入式(1)、式(2)得到水稻插秧機的故障概率密度函數(shù)f(t)、可靠度函數(shù)R(t)。繪制故障概率密度函數(shù)f(t)與t的關系圖,如圖3所示。

圖3 水稻插秧機失效概率密度曲線Fig.3 Failure probability density curve of rice transplanter

(13)

(14)

在表1數(shù)據(jù)的基礎上利用最小二乘法計算出單一威布爾模型的可靠度函數(shù)為

(15)

單一威布爾模型、混合威布爾模型的可靠度曲線與失效數(shù)據(jù)擬合效果如圖4所示。

圖4 水稻插秧機可靠度曲線Fig.4 Reliability curve of rice transplanter

從圖4可以看出,混合威布爾分布與水稻插秧機失效數(shù)據(jù)樣本的擬合效果要好于單一威布爾分布,能更準確描述水稻插秧機的失效規(guī)律,證明了所選兩參數(shù)混合威布爾模型的合理性。水稻插秧機的可靠性評價指標包括平均壽命(MTBF)、中位壽命t(0.5)、特征壽命t(e-1)。

(16)

式(16)中:Γ為Gamma函數(shù)。

將表2中計算出的模型參數(shù)代入式(16)得到該水稻插秧機的平均壽命MTBF=161.75 h,通過式(14)可計算出:t(0.5)=147.14 h,t(e-1)=191.31 h。為了保證插秧機的可靠性,取可靠度R(t)為0.6來計算預防性維修周期,取不同的可靠度代入式(14)進行計算,得到不同可靠度下的水稻插秧機預防性維修周期如表4所示。

表4 不同可靠度下的插秧機預防性維修周期

根據(jù)表4可知,水稻插秧機的預防性維修周期隨著R(t)的減小而增大,在R(t)為0.6時,預防性維修周期為115.19 h,水稻插秧機的平均壽命MTBF=161.75 h, 115.19 h小于插秧機的平均壽命,表明規(guī)定的預防性維修周期較為合理。

3.2 剩余壽命預測

在可靠性理論中,具有年齡t的設備從t時刻開始繼續(xù)使用下去直到失效為止所經歷的時間稱為具有年齡t的設備的剩余壽命[14],壽命模型能否準確描述產品的壽命分布,剩余壽命是關鍵,通過對剩余壽命進行預測,可以得知運行時長與剩余壽命之間的具體關系。

假設T是水稻插秧機的失效時間,水稻插秧機已經運行了時間t,那么剩余使用壽命時間為條件隨機變量x,其中t

(17)

式(17)中:x為剩余壽命時間。

將式(12)代入式(17),計算出水稻插秧機的剩余使用壽命x的可靠度函數(shù)Rt(x)為

(18)

現(xiàn)以A、B、C、D四臺正在運行的水稻插秧機為例進行剩余壽命預測,A、B、C、D四臺水稻插秧機的工作時長t分別為50、100、200、300 h,將t代入式(18)中,可得水稻插秧機剩余壽命可靠度曲線如圖5所示。

從圖5中可以看出,隨著水稻插秧機運行時間的增加,相同可靠度下的剩余壽命相應得減少;運行時間越長,其可靠度曲線越陡,與實際情況相符合,水稻插秧機運行時間越長,受到的磨損、沖擊越嚴重,那么插秧機剩余使用壽命就越低,表明水稻插秧機剩余壽命預測模型的有效性。根據(jù)式(18)可計算得到A、B、C、D水稻插秧機的中位壽命分別為129.47、107.21、73.72、54.11 h,即進行插秧機進行預測性維護時,優(yōu)先對工作年齡達到300 h的D水稻插秧機進行維護。

4 結論

針對單一威布爾模型不能準確評估含多種失效模式的水稻插秧機,利用兩重兩參數(shù)混合威布爾分布對失效數(shù)據(jù)進行擬合,通過改進粒子群算法和最小二乘法相結合求解出模型參數(shù)。然后通過K-S檢驗以及比較單一威布爾和混合威布爾的可靠度曲線與失效數(shù)據(jù)的擬合效果證明了混合威布爾模型的合理性,并利用該模型計算得到相應的可靠性指標,最后在混合威布爾模型基礎上對水稻插秧機進行剩余壽命預測,得到水稻插秧機運行時間和剩余壽命之間的對應關系。

研究過程及結果如下。

(1)針對單一威布爾模型不能準確描述多種故障機制并存的水稻插秧機的失效規(guī)律,利用混合威布爾模型進行評估,提高評估的準確性。

(2)與傳統(tǒng)的粒子群算法相比,利用改進粒子群算法計算得到的參數(shù)值更加準確,得到了更準確分布模型。在此模型基礎上計算出了水稻插秧機的可靠性指標以及預防性維修周期。

(3)在更準確的混合威布爾模型的基礎上計算得到插秧機的剩余使用壽命,便于使用和維修人員安排合理的修理計劃。

后續(xù)研究過程中,將水稻插秧機分為多個子系統(tǒng),針對故障多發(fā)的子系統(tǒng)分別進行可靠性分析,使得整機的可靠性分析更加準確。