楊皓潔, 呂建偉, 徐一帆

(海軍工程大學(xué)管理工程與裝備經(jīng)濟(jì)系, 湖北 武漢 430033)

0 引 言

具有復(fù)雜冗余結(jié)構(gòu)、時(shí)間相關(guān)故障的多狀態(tài)系統(tǒng)是當(dāng)前可靠性領(lǐng)域的熱點(diǎn)研究對象之一。復(fù)雜冗余結(jié)構(gòu)中的元件常表現(xiàn)出動態(tài)故障行為,相應(yīng)元件在啟用前后具有不同的失效率;而時(shí)間相關(guān)故障意味著故障率具有時(shí)變的特點(diǎn),狀態(tài)轉(zhuǎn)移時(shí)間服從非指數(shù)分布。在該條件下故障率是時(shí)間的函數(shù),因此故障與時(shí)間相關(guān)。

考慮這些因素會使系統(tǒng)故障演化規(guī)律和任務(wù)過程更為復(fù)雜,增大了系統(tǒng)可靠性及任務(wù)成功性評估的難度。系統(tǒng)可靠性評估的研究方法大致可分為解析法和仿真法。其中,解析法有組合模型法[1-2]、狀態(tài)空間法[3-4]、概率近似法[5]等,具有計(jì)算精度高、計(jì)算結(jié)果理想可靠等特點(diǎn)。組合模型法中,以故障樹分析為基礎(chǔ)的多值決策圖[6-8](布爾擴(kuò)展模型法)、序列二元決策圖等[9-10]是解決多狀態(tài)系統(tǒng)和動態(tài)系統(tǒng)可靠性的較好方法,但缺點(diǎn)是單獨(dú)使用此類方法時(shí),不能處理可修系統(tǒng),對于系統(tǒng)中的可修部分必須結(jié)合狀態(tài)空間法。另一類主流方法是通用生成函數(shù)[11-12],該方法通過引入z變換函數(shù)、通用生成算子等工具,再結(jié)合系統(tǒng)元件的狀態(tài)概率分布,計(jì)算得到系統(tǒng)的可靠性指標(biāo)。通用生成函數(shù)是解決多狀態(tài)系統(tǒng)可靠性問題的高效方法,但當(dāng)系統(tǒng)動態(tài)規(guī)律較為復(fù)雜時(shí),仍存在建模和求解困難的問題。當(dāng)面對時(shí)間相關(guān)故障問題時(shí),需要利用半馬爾可夫模型進(jìn)行建模,但該模型在建模和計(jì)算求解上都較為復(fù)雜。此外,解析法在面對時(shí)間相關(guān)故障問題時(shí)同樣存在建模和求解復(fù)雜的問題。仿真法作為一類評估復(fù)雜系統(tǒng)可靠性及任務(wù)成功性的統(tǒng)計(jì)類方法[13-15],主要包括基于離散事件仿真[16]、貝葉斯理論[17-18]、petri網(wǎng)[19]等建模方法。其優(yōu)點(diǎn)是易于建模且模型通用性較強(qiáng)。但由于仿真抽樣效率在針對具體問題時(shí)各不相同,計(jì)算精度和仿真次數(shù)之間的矛盾難以兼顧。目前主要通過特定的方差減小方法[20-23]來提高計(jì)算結(jié)果的精度,但相關(guān)研究中均以元件狀態(tài)轉(zhuǎn)移時(shí)間服從指數(shù)分布為前提，該條件下不需考慮狀態(tài)轉(zhuǎn)移時(shí)間的剩余分布問題,并且系統(tǒng)狀態(tài)轉(zhuǎn)移時(shí)間也服從指數(shù)分布,容易實(shí)現(xiàn)系統(tǒng)狀態(tài)轉(zhuǎn)移時(shí)間的直接抽樣,在簡化問題的同時(shí)限制了方法的使用范圍。一旦考慮元件狀態(tài)轉(zhuǎn)移時(shí)間不服從指數(shù)分布,系統(tǒng)狀態(tài)轉(zhuǎn)移時(shí)間的概率分布將變得相當(dāng)復(fù)雜,導(dǎo)致在對其進(jìn)行抽樣時(shí),相關(guān)計(jì)算較為困難。

為了在考慮時(shí)間相關(guān)故障、多狀態(tài)性能輸出、系統(tǒng)冗余儲備等因素的同時(shí),改進(jìn)抽樣效率和估值精度,本文提出了一種基于蒙特卡羅離散事件仿真,冗余多狀態(tài)系統(tǒng)可靠性及任務(wù)成功性的評估方法,用以針對故障率具有時(shí)變特點(diǎn)的時(shí)間相關(guān)故障。通過設(shè)計(jì)不服從指數(shù)分布狀態(tài)轉(zhuǎn)移隨機(jī)抽樣的數(shù)值算法,實(shí)現(xiàn)故障事件的隨機(jī)抽樣。在此基礎(chǔ)上,結(jié)合受迫轉(zhuǎn)移、失效偏倚等方差減小方法[24-27],改善小概率事件抽樣效率和估值精度,為該問題提出一種較為完善的解決方案。

1 問題描述及分析

某多狀態(tài)系統(tǒng)S由n臺不同功能的設(shè)備組成,其中設(shè)備j有mj種性能狀態(tài)。任意設(shè)備的狀態(tài)變化都會引起系統(tǒng)狀態(tài)轉(zhuǎn)移,因此在任務(wù)過程中的某一時(shí)刻t,系統(tǒng)狀態(tài)S(t)及設(shè)備性能狀態(tài)Gj(t)的關(guān)系如下:

S(t)={G1(t),G2(t),…,Gn(t)}

(1)

(2)

(3)

式中:s為系統(tǒng)狀態(tài)空間;gj為設(shè)備的性能狀態(tài)空間。實(shí)際上,系統(tǒng)狀態(tài)就是由各設(shè)備狀態(tài)組成的向量,但不同的系統(tǒng)狀態(tài)可能有相同的性能輸出表現(xiàn),因此定義系統(tǒng)性能狀態(tài)G(t)為

G(t)=Φ(S(t))∈g={g1,g2,…,gm}

(4)

式中:Φ(·)為某具體系統(tǒng)的結(jié)構(gòu)和運(yùn)行方式確定的構(gòu)型函數(shù);g為系統(tǒng)的性能狀態(tài)空間。設(shè)備的性能狀態(tài)轉(zhuǎn)移關(guān)系如圖1所示。

圖1 設(shè)備性能狀態(tài)轉(zhuǎn)移圖

圖1中各類狀態(tài)及對應(yīng)編號如下:完好為n、完好到故障間的中間狀態(tài)為(n-1)到2、一般故障為1、致命故障為0。其中完好對應(yīng)性能輸出為100%。中間狀態(tài)對應(yīng)性能輸出為0%～100%。一般故障、致命故障對應(yīng)性能輸出均為0%,其區(qū)別在于一般故障為可修復(fù)故障,致命故障為無法修復(fù)或任務(wù)時(shí)間內(nèi)不具有維修條件的故障,也標(biāo)為吸收態(tài)。λi,i′和μi,i′分別為該設(shè)備從狀態(tài)i′→i的故障率和修復(fù)率。

為達(dá)到提高系統(tǒng)可靠性的目的,系統(tǒng)中存在一些冗余結(jié)構(gòu),包括冷、溫、熱儲備。由冗余結(jié)構(gòu)和系統(tǒng)運(yùn)行產(chǎn)生設(shè)備的工作狀態(tài)大致可分為3類:工作、儲備待機(jī)和維修。設(shè)備工作狀態(tài)根據(jù)系統(tǒng)結(jié)構(gòu)、備用設(shè)備啟動規(guī)則和設(shè)備當(dāng)前狀態(tài)來確定。假設(shè)x、y分別為一個(gè)冗余結(jié)構(gòu)中的主設(shè)備和備用設(shè)備,圖2展示了冷、溫儲備結(jié)構(gòu)下的主設(shè)備與備用設(shè)備工作狀態(tài)的轉(zhuǎn)移關(guān)系。熱儲備結(jié)構(gòu)較為簡單(即均視為工作),在此不再贅述。

圖2 冷、溫儲備設(shè)備工作狀態(tài)轉(zhuǎn)移規(guī)則圖

圖2中Gx(t)和Gy(t)表示設(shè)備x和y在t時(shí)刻的性能狀態(tài),數(shù)值為圖1中性能狀態(tài)的編號;Sysstate表示系統(tǒng)狀態(tài);down表示系統(tǒng)宕機(jī),available表示系統(tǒng)可用。當(dāng)主設(shè)備性能狀態(tài)低于備用設(shè)備或發(fā)生一般故障時(shí)應(yīng)當(dāng)轉(zhuǎn)入維修,同時(shí)備用設(shè)備轉(zhuǎn)入工作狀態(tài);主設(shè)備維修完成后立即轉(zhuǎn)入工作狀態(tài),同時(shí)備用設(shè)備轉(zhuǎn)入待機(jī)狀態(tài);當(dāng)系統(tǒng)中其他部分性能退化導(dǎo)致系統(tǒng)宕機(jī)時(shí),主、備單元均轉(zhuǎn)入待機(jī)狀態(tài)。根據(jù)此規(guī)則并結(jié)合可靠性框圖即可確定所有設(shè)備的工作狀態(tài)。在動態(tài)系統(tǒng)中，某設(shè)備的狀態(tài)由工作狀態(tài)和性能狀態(tài)共同定義。其意義在于確定系統(tǒng)中所有的狀態(tài)轉(zhuǎn)移速率,也是蒙特卡羅隨機(jī)事件抽樣的前提。

2 基于蒙特卡羅仿真的模型描述

2.1 仿真流程設(shè)計(jì)

本文所采用的蒙特卡羅離散事件仿真方法流程如圖3所示。在單次任務(wù)仿真過程中,需進(jìn)行多次離散事件抽樣循環(huán)。一次完整的離散事件抽樣步驟如下:

步驟 1在系統(tǒng)最底層(設(shè)備級)進(jìn)行系統(tǒng)狀態(tài)轉(zhuǎn)移隨機(jī)抽樣,得到故障或維修事件發(fā)生的時(shí)間和轉(zhuǎn)入狀態(tài),這是仿真的核心步驟,如圖3中灰色部分所示;

圖3 總體仿真流程圖

步驟 2以底層設(shè)備狀態(tài)變化驅(qū)動上層系統(tǒng)狀態(tài)變化;

步驟 3根據(jù)系統(tǒng)狀態(tài)變化更新任務(wù)相關(guān)變量;

步驟 4根據(jù)任務(wù)相關(guān)變量和系統(tǒng)狀態(tài)判斷任務(wù)是否失敗。

在完成仿真后整理數(shù)據(jù)得到統(tǒng)計(jì)結(jié)果。設(shè)任務(wù)成功概率為D,系統(tǒng)瞬時(shí)可用度為A(t),則有

D=(N-Nd)/N

(5)

A(t)=[N-Nf(t)]/N

(6)

式中:Nd為任務(wù)失敗頻數(shù);Nf(t)為t時(shí)刻系統(tǒng)處于不可用狀態(tài)頻數(shù)。Nd和Nf(t)的計(jì)數(shù)與系統(tǒng)狀態(tài)轉(zhuǎn)移的權(quán)重wt有關(guān)。wt是一個(gè)與時(shí)間相關(guān)的變量,其含義為t時(shí)刻前所有隨機(jī)事件對之后事件發(fā)生的貢獻(xiàn)程度。因此,wt代表t時(shí)刻前每一次狀態(tài)轉(zhuǎn)移累積的權(quán)重,同時(shí)代表t時(shí)刻某隨機(jī)事件的仿真計(jì)數(shù)值。任務(wù)開始時(shí)令w0=1。具體wt的計(jì)算與系統(tǒng)狀態(tài)轉(zhuǎn)移隨機(jī)抽樣在第2.2節(jié)、第2.3節(jié)一并展開。當(dāng)任務(wù)失敗時(shí)更新任務(wù)失敗頻數(shù):

Nd←Nd+wT0

(7)

當(dāng)t時(shí)刻系統(tǒng)處于不可用狀態(tài)時(shí)更新系統(tǒng)不可用狀態(tài)頻數(shù):

Nf(t)←Nf(t)+wt

(8)

不同的抽樣方法所產(chǎn)生的隨機(jī)事件的權(quán)重可能是不同的,wt可能不為1,因此Nd和Nf(t)也可能不是整數(shù)。

2.2 系統(tǒng)狀態(tài)轉(zhuǎn)移隨機(jī)抽樣

蒙特卡羅方法通過生成隨機(jī)事件,模擬系統(tǒng)的動態(tài)行為。本研究使用的基本抽樣方法是蒙特卡羅間接仿真法,在此基礎(chǔ)上才能夠使用方差減小方法來提高仿真效率及減小統(tǒng)計(jì)誤差。蒙特卡羅間接仿真法的基本步驟是:先隨機(jī)生成系統(tǒng)當(dāng)前狀態(tài)的駐留時(shí)間(或狀態(tài)轉(zhuǎn)移的發(fā)生時(shí)間),再通過輪盤賭等方式隨機(jī)確定觸發(fā)事件對應(yīng)的設(shè)備和該設(shè)備轉(zhuǎn)移時(shí)進(jìn)入的狀態(tài),如圖4所示。

圖4 多狀態(tài)系統(tǒng)狀態(tài)轉(zhuǎn)移的蒙特卡羅間接仿真法抽樣過程

對于多狀態(tài)系統(tǒng)而言,分系統(tǒng)、設(shè)備存在多個(gè)性能狀態(tài),各狀態(tài)之間的轉(zhuǎn)移概率服從特定概率分布。系統(tǒng)在某一給定時(shí)刻發(fā)生狀態(tài)轉(zhuǎn)移并進(jìn)入新的狀態(tài)是由一個(gè)狀態(tài)概率轉(zhuǎn)移核(下文簡稱為“轉(zhuǎn)移核”)決定的,轉(zhuǎn)移核可以全面地描述系統(tǒng)在時(shí)域內(nèi)的隨機(jī)行為。轉(zhuǎn)移核本質(zhì)上是一個(gè)概率密度函數(shù)。多狀態(tài)系統(tǒng)的轉(zhuǎn)移核可以表示為

K((k,t)|(k′,t′))=f(t|k′,t′)q(k|k′,t)

(9)

這里用f(t|k′,t′)表示在時(shí)刻t′發(fā)生下一次狀態(tài)轉(zhuǎn)移時(shí)刻為t的概率密度函數(shù),用q(k|k′,t)表示在時(shí)刻t時(shí)狀態(tài)k′→k的轉(zhuǎn)移概率,且有

(10)

(11)

(12)

式中：γk′(t)表示當(dāng)前處于時(shí)刻t時(shí)系統(tǒng)轉(zhuǎn)移出當(dāng)前狀態(tài)k′的轉(zhuǎn)移速率;γkk′(t)表示當(dāng)前處于時(shí)刻t時(shí)狀態(tài)k′→k的轉(zhuǎn)移速率;k′和k均屬于系統(tǒng)狀態(tài)空間s。系統(tǒng)經(jīng)過多次狀態(tài)轉(zhuǎn)移形成了一條狀態(tài)的隨機(jī)游動序列ΓM:(k0,t0),(k1,t1),…,(kM,tM),以實(shí)現(xiàn)離散事件抽樣的循環(huán)。

在轉(zhuǎn)移核的基礎(chǔ)上就可以對轉(zhuǎn)移發(fā)生時(shí)刻t和轉(zhuǎn)移進(jìn)入狀態(tài)k進(jìn)行抽樣,蒙特卡羅間接仿真在正常狀態(tài)轉(zhuǎn)移空間的概率分布下進(jìn)行抽樣,因此系統(tǒng)狀態(tài)轉(zhuǎn)移權(quán)重為

wt=wt′

(13)

(1)轉(zhuǎn)移發(fā)生時(shí)刻t的抽樣

發(fā)生狀態(tài)轉(zhuǎn)移時(shí)刻t計(jì)算公式為

(14)

式中：ξ為0～1之間均勻分布的隨機(jī)數(shù)。為了描述的簡潔,這里將故障、降級、退化等系統(tǒng)狀態(tài)變化統(tǒng)稱為故障類隨機(jī)事件,將維護(hù)、修理等系統(tǒng)狀態(tài)變化統(tǒng)稱為修復(fù)性隨機(jī)事件。根據(jù)常規(guī)標(biāo)識習(xí)慣,當(dāng)轉(zhuǎn)移后系統(tǒng)發(fā)生故障或退化時(shí),轉(zhuǎn)移速率用λ表示;當(dāng)轉(zhuǎn)移后系統(tǒng)修復(fù)時(shí),轉(zhuǎn)移速率用μ表示;Ok′表示狀態(tài)k′的可達(dá)狀態(tài)集;用Λk′表示狀態(tài)k′可達(dá)的故障或降級狀態(tài)集,Ψk′表示可達(dá)的修復(fù)狀態(tài)集:

(15)

Ok′=Λk′∪Ψk′

(16)

由于本文假設(shè)設(shè)備當(dāng)前狀態(tài)的駐留時(shí)間均服從非指數(shù)分布,因此系統(tǒng)狀態(tài)轉(zhuǎn)移的抽樣帶有一定復(fù)雜性,具體內(nèi)容在下文中展開。

(2)轉(zhuǎn)移進(jìn)入狀態(tài)k的抽樣

用輪盤賭的方式可以確定轉(zhuǎn)移到達(dá)的具體狀態(tài):

(17)

式中：ζ為0～1之間均勻分布的隨機(jī)數(shù)。當(dāng)ζ落在式(17)區(qū)間內(nèi),則系統(tǒng)從狀態(tài)k′轉(zhuǎn)移至狀態(tài)k,也就是確定了發(fā)生狀態(tài)轉(zhuǎn)移的設(shè)備和該設(shè)備轉(zhuǎn)入的狀態(tài)。蒙特卡羅間接仿真法利用兩個(gè)隨機(jī)數(shù)就可以產(chǎn)生一次狀態(tài)轉(zhuǎn)移,而且不需要比較等步驟,一定程度上提高了計(jì)算效率,在此基礎(chǔ)上可使用方差減小方法進(jìn)行加速抽樣。

2.3 方差減小方法

(18)

(19)

(1)受迫轉(zhuǎn)移

(20)

(21)

可見受迫轉(zhuǎn)移法只改變轉(zhuǎn)移時(shí)間的抽樣空間,不改變轉(zhuǎn)移狀態(tài)的抽樣空間。然后按照蒙特卡羅間接仿真法進(jìn)行轉(zhuǎn)移發(fā)生時(shí)刻t的抽樣,此時(shí)式(14)變?yōu)?/p>

(22)

受迫轉(zhuǎn)移法使?fàn)顟B(tài)轉(zhuǎn)移時(shí)刻發(fā)生在任務(wù)結(jié)束時(shí)間內(nèi),客觀上起到了強(qiáng)制提升隨機(jī)事件發(fā)生概率的作用,從而增加了隨機(jī)事件的抽樣實(shí)現(xiàn)次數(shù)。為使估計(jì)值保持無偏性,修正仿真抽樣計(jì)數(shù)：

(23)

(2)失效偏倚

(24)

(25)

式中:x為失效偏倚系數(shù),用于調(diào)節(jié)傾向故障類轉(zhuǎn)移及修復(fù)類轉(zhuǎn)移之間的相對程度?？梢娛蟹ㄖ桓淖冝D(zhuǎn)移狀態(tài)抽樣空間,不改變轉(zhuǎn)移時(shí)間抽樣空間。x的具體含義是本次抽樣中故障類轉(zhuǎn)移占所有狀態(tài)轉(zhuǎn)移的比例;(1-x)是本次抽樣中維修類轉(zhuǎn)移占所有狀態(tài)轉(zhuǎn)移的比例,通常取值范圍在0.5～0.7[21]。這時(shí)系統(tǒng)轉(zhuǎn)移到故障狀態(tài)的轉(zhuǎn)移速率將遠(yuǎn)大于不偏倚抽樣時(shí)的轉(zhuǎn)移速率。然后根據(jù)式(17)確定轉(zhuǎn)移進(jìn)入的狀態(tài)k。為了使估計(jì)值保持無偏性應(yīng)當(dāng)修正仿真抽樣計(jì)數(shù):

wt=

(26)

(3)抽樣方法的使用條件

每1次隨機(jī)事件抽樣只能從蒙特卡羅間接仿真、受迫轉(zhuǎn)移、失效偏倚3種方法中選用一種,根據(jù)每種方法的特點(diǎn)可在仿真中做以下規(guī)定:① 每次任務(wù)初始化后第1次狀態(tài)轉(zhuǎn)移(k0,t0)→(k1,t1)使用受迫轉(zhuǎn)移法;② 當(dāng)系統(tǒng)中存在維修類轉(zhuǎn)移時(shí)使用失效偏倚法;③ 其余情況使用蒙特卡羅間接仿真法。

3 非指數(shù)分布隨機(jī)抽樣的數(shù)值算法

由于系統(tǒng)中存在服從非指數(shù)分布的狀態(tài)轉(zhuǎn)移,其轉(zhuǎn)移率為時(shí)間的函數(shù),因此系統(tǒng)當(dāng)前狀態(tài)的駐留時(shí)間無“記憶性”,是典型的非齊次馬爾可夫過程。在離散事件仿真中,隨著每次仿真時(shí)間的推進(jìn),一直正常工作或一直處于維修狀態(tài)的設(shè)備的當(dāng)前狀態(tài)駐留時(shí)間應(yīng)按其剩余分布進(jìn)行抽樣。以某設(shè)備的故障轉(zhuǎn)移為例,假設(shè)其故障時(shí)間T的初始分布為F(t),則在累計(jì)工作tc后,其剩余分布函數(shù)[28]為

(27)

作為減小方差的基礎(chǔ)方法,蒙特卡羅間接仿真法首先需要對系統(tǒng)當(dāng)前狀態(tài)駐留時(shí)間進(jìn)行抽樣,則考慮剩余分布后的系統(tǒng)狀態(tài)駐留時(shí)間所服從的分布函數(shù)可表示為

FS-RD(t)=P{min(T1,T2,…,TK)

(28)

步驟 1選取一個(gè)足夠大的數(shù)ts,使得FS-RD(ts)=1。盡管FS-RD(ts)在理論上只能無限接近于1,但在計(jì)算機(jī)上可以實(shí)現(xiàn)等于1。

步驟 2在[0,tu]上均勻選擇若干點(diǎn)組成向量Ts,并根據(jù)式(28)計(jì)算相應(yīng)時(shí)間的概率值Ps=FS-RD(Ts)。由此組成離散樣本點(diǎn)組合(Ts,Ps)=[(t0,p0),(t1,p1), …,(ts,ps)]。

步驟 3對組合(Ts,Ps)進(jìn)行補(bǔ)充和刪選。確定向量Ts中使得FS(tu)>0.99的最小元素tu,并在[0,tu]上選擇若干點(diǎn)以步驟2相同的方式補(bǔ)充離散樣本點(diǎn)組合(Ts,Ps)。Ps中可能存在重復(fù)元素,只保留Ps中不重復(fù)元素和Ts中對應(yīng)的元素,刪除其余元素,為插值法做準(zhǔn)備。

步驟 4利用插值法進(jìn)行隨機(jī)抽樣,以概率Ps為自變量,時(shí)間Ts為因變量。生成0～1之間均勻分布的隨機(jī)數(shù)ξ,若此次抽樣方法為受迫轉(zhuǎn)移法,則生成一個(gè)[0,FS-RD(T0-t′)]間的隨機(jī)數(shù)ξ。在組合(Ts,Ps)中選取自變量與ξ相鄰的兩個(gè)點(diǎn)(ti,pi)、(ti+1,pi+1),再進(jìn)行線性插值計(jì)算得到系統(tǒng)當(dāng)前狀態(tài)駐留時(shí)間Δt：

(29)

該算法不需大量樣本點(diǎn)就可以完整地展示分布函數(shù)全貌,在有限的計(jì)算資源下提高了插值法精度,還可以降低計(jì)算量。并且狀態(tài)轉(zhuǎn)移時(shí)間分布函數(shù)FS-RD(t)是單調(diào)遞增函數(shù),能保證上述算法Δt有唯一值。

以威布爾分布為例,在仿真開始時(shí)令累計(jì)工作時(shí)間tc=0,正常工作及熱儲備設(shè)備的tc應(yīng)累加仿真推進(jìn)時(shí)間,發(fā)生故障需維修的設(shè)備的tc清零,同時(shí)更新該設(shè)備相應(yīng)的狀態(tài)轉(zhuǎn)移率:

γ(tc)=bη-b(tc)b-1

(30)

式中：b為形狀參數(shù);η為尺度參數(shù)。根據(jù)式(27),得到該失效機(jī)理發(fā)生時(shí)間剩余分布為

(31)

系統(tǒng)中存在相互獨(dú)立、相互競爭的多種服從威布爾分布的失效機(jī)理和維修轉(zhuǎn)移,bkk′和ηkk′為相關(guān)轉(zhuǎn)移的形狀參數(shù)和尺度參數(shù)。根據(jù)式(28)可得系統(tǒng)的當(dāng)前狀態(tài)駐留時(shí)間分布函數(shù)。

(32)

圖5展示了該算法的抽樣結(jié)果。在隨機(jī)抽樣生成系統(tǒng)駐留時(shí)間Δt后,根據(jù)式(17)進(jìn)行輪盤賭確定最終發(fā)生狀態(tài)轉(zhuǎn)移的設(shè)備和具體轉(zhuǎn)移狀態(tài)。

圖5 非指數(shù)分布隨機(jī)抽樣數(shù)值算法結(jié)果

4 算例分析

4.1 算例背景

本算例以艦船及其裝備系統(tǒng)為對象進(jìn)行分析,由于篇幅所限,這里僅以艦船動力、電力系統(tǒng)在航渡任務(wù)階段的任務(wù)可靠性和任務(wù)成功率評估為例,其他系統(tǒng)及任務(wù)的相應(yīng)計(jì)算也同樣適用。設(shè)某艦動力、電力系統(tǒng)的系統(tǒng)結(jié)構(gòu)組成如圖6所示。相關(guān)設(shè)備的性能狀態(tài)轉(zhuǎn)移關(guān)系如圖7所示。屬于第I類設(shè)備的有柴油機(jī)、燃?xì)廨啓C(jī),這類設(shè)備有多種性能輸出表現(xiàn),且可能發(fā)生致命故障;屬于第II類設(shè)備的有柴油發(fā)電機(jī)組、軸系、螺旋槳,這類設(shè)備只有完好和失效兩種性能輸出表現(xiàn),且可能發(fā)生致命故障;屬于第Ⅲ類設(shè)備的有監(jiān)控設(shè)備、輔助設(shè)備。在本文應(yīng)用背景中,假設(shè)這類設(shè)備不具有影響任務(wù)成敗的致命故障,可將這類設(shè)備簡化為可修的二態(tài)設(shè)備。系統(tǒng)中的冗余結(jié)構(gòu)包括冷儲備結(jié)構(gòu)(柴油機(jī)1,柴油機(jī)2和燃?xì)廨啓C(jī))和熱儲備結(jié)構(gòu)(軸系和螺旋槳1,螺旋槳2),其中冷儲備結(jié)構(gòu)按照圖2所示規(guī)則啟用備用設(shè)備。

圖6 航渡任務(wù)中艦船系統(tǒng)結(jié)構(gòu)

圖7 相關(guān)設(shè)備性能狀態(tài)轉(zhuǎn)移圖

該艦在航渡階段的主要任務(wù)目標(biāo)[29]為:① 按時(shí)抵達(dá)目標(biāo)海域;② 當(dāng)航渡任務(wù)結(jié)束時(shí)系統(tǒng)仍處于可工作(非故障)狀態(tài)。目標(biāo)①是對完成任務(wù)時(shí)間和空間的要求,指在規(guī)定的任務(wù)時(shí)間T0里完成規(guī)定的任務(wù)航程S0。艦船既不可提前到達(dá),也不可晚于規(guī)定時(shí)間到達(dá)。提前到達(dá)可能導(dǎo)致本艦脫離協(xié)同兵力、提前暴露行蹤;晚于預(yù)定時(shí)間到達(dá)則難以實(shí)現(xiàn)兵力協(xié)同、完成既定作戰(zhàn)任務(wù)。目標(biāo)②是對系統(tǒng)下一階段完成作戰(zhàn)任務(wù)的要求?？傊?當(dāng)且僅當(dāng)以上兩項(xiàng)目標(biāo)均滿足時(shí),方可判定任務(wù)成功(實(shí)際情況下的任務(wù)成功性問題判據(jù)可能更為復(fù)雜,此處略去)。據(jù)此可知,任務(wù)目標(biāo)是否達(dá)成應(yīng)結(jié)合任務(wù)成功性判定規(guī)則,根據(jù)特定的依據(jù)給予評判,如表1所示。在表1中,航速超限是指由于經(jīng)歷系統(tǒng)降級或故障維修,導(dǎo)致該艦在剩余航程中即使以最大航速也不能完成任務(wù)。維修超時(shí)是指任務(wù)結(jié)束時(shí)由于設(shè)備故障正在維修,導(dǎo)致該艦總體處于宕機(jī)、不可用的狀態(tài)。當(dāng)發(fā)生維修超時(shí)時(shí),艦船通常還未到達(dá)目標(biāo)區(qū)域。致命故障是指由于設(shè)備致命故障導(dǎo)致該艦總體宕機(jī)且無法修復(fù),從而無法完成剩余任務(wù)。

表1 航渡任務(wù)成功判據(jù)

4.2 模型假設(shè)條件

根據(jù)上文分析對該案例相關(guān)分系統(tǒng)做進(jìn)一步假設(shè)。

假設(shè) 1在航渡階段,其他分系統(tǒng)不影響該艦總體狀態(tài)。

假設(shè) 2分系統(tǒng)和艦船的可用狀態(tài)包括完好和降級狀態(tài),系統(tǒng)降級時(shí)會影響該艦航速。

假設(shè) 3所有狀態(tài)轉(zhuǎn)移均為獨(dú)立狀態(tài)轉(zhuǎn)移。

假設(shè) 4設(shè)備維修為完全維修(修舊如新),修理后狀態(tài)均返回完好狀態(tài),且維修完成后故障率回歸初始狀態(tài),即λ(t)=λ(0)。

假設(shè) 5不考慮維修資源限制。

4.3 任務(wù)結(jié)果判定

表1中3種任務(wù)失敗類型可通過以下方法進(jìn)行判定。

(1)航速超限

設(shè)S0為任務(wù)航程;St為當(dāng)前累計(jì)航程;T0為任務(wù)總時(shí)間;t為當(dāng)前時(shí)刻;?為航程松弛系數(shù)[30](0

t=t′+Δt

(33)

St=St′+vΔt

(34)

(35)

結(jié)合實(shí)際情況,艦船通常具有特定的經(jīng)濟(jì)巡航航速。航速與動力分系統(tǒng)的功率輸出情況和航速規(guī)則有關(guān),當(dāng)剩余航程計(jì)劃航速vp低于veco時(shí),則以veco航速航行,此次隨機(jī)事件發(fā)生后的實(shí)際航速v為

(36)

式中：θ為功率輸出系數(shù)(0<θ<1),表示艦船在某狀態(tài)下的最大功率占完好狀態(tài)下滿功率的比值,具體取值根據(jù)動力系統(tǒng)狀態(tài)確定。根據(jù)vp與vmax的關(guān)系判定是否航速超限。需要指出的是,在仿真中當(dāng)動力分系統(tǒng)功率輸出降級時(shí),可能暫時(shí)出現(xiàn)vmax>vp>v的情況,此時(shí)暫不判定失敗。當(dāng)修理完成后,功率輸出恢復(fù)完好狀態(tài),使得vpvmax,則認(rèn)定為航速超限,判定任務(wù)失敗。

(2)維修超時(shí)

若本次隨機(jī)事件的發(fā)生時(shí)間超出任務(wù)時(shí)間,即t>T0,并且此時(shí)系統(tǒng)處于宕機(jī)維修狀態(tài),則認(rèn)定為維修超時(shí),判定任務(wù)失敗。

(3)致命故障

若某設(shè)備發(fā)生致命故障,導(dǎo)致系統(tǒng)進(jìn)入吸收態(tài),即系統(tǒng)無當(dāng)場修復(fù)的可能,則認(rèn)定為致命故障,判定任務(wù)失敗。

由于存在多種任務(wù)失敗類型,因此可擴(kuò)展Nd為向量，Nd=[NOS,NUM,NFF],其中元素分別對應(yīng)航速超限、維修超時(shí)和致命故障的頻數(shù)。按照下式計(jì)算不同失敗類型概率,以便于觀察其對任務(wù)成功性的影響。

Dd=Nd/N

(37)

4.4 算例數(shù)據(jù)設(shè)置

(1)系統(tǒng)參數(shù)

通常機(jī)械設(shè)備磨損、電子元件老化等退化失效模式是造成設(shè)備一般故障和降級退化的主要原因,由此可認(rèn)為其失效、退化時(shí)間服從形狀參數(shù)大于1的威布爾分布。而發(fā)生致命故障的隨機(jī)性較大,故此算例將其失效時(shí)間設(shè)為指數(shù)分布(即形狀參數(shù)等于1的威布爾分布)。設(shè)備的維修時(shí)間通常為正態(tài)分布,但正態(tài)分布下可能得到負(fù)數(shù)的抽樣值,這與維修時(shí)間非負(fù)性的要求不符。因此,通常用威布爾分布或?qū)?shù)正態(tài)分布來近似模擬正態(tài)分布。本文假設(shè)維修時(shí)間服從威布爾分布。當(dāng)威布爾分布的形狀參數(shù)大于3時(shí),形態(tài)逼近于正態(tài)分布。例如某隨機(jī)變量X1～N(50,100)可用X2～Wbl(54, 5.8)代替,其概率密度函數(shù)如圖8所示。

圖8 正態(tài)分布與威布爾分布概率密度函數(shù)對比

系統(tǒng)中各設(shè)備狀態(tài)轉(zhuǎn)移時(shí)間所服從威布爾分布的尺度參數(shù)η和形狀參數(shù)b如表2所示。當(dāng)系統(tǒng)處于可用狀態(tài)時(shí),式(36)中功率輸出系數(shù)θ主要由柴油機(jī)、燃?xì)廨啓C(jī)狀態(tài)及功能結(jié)構(gòu)決定,設(shè)備性能狀態(tài)與功率輸出系數(shù)關(guān)系如表3所示。可見系統(tǒng)也會表現(xiàn)出介于完好與故障間的降級退化狀態(tài),這時(shí)最大航速只有完好情況下的60%。

表2 設(shè)備狀態(tài)轉(zhuǎn)移分布的參數(shù)設(shè)置

表3 動力系統(tǒng)設(shè)備性能狀態(tài)與功率輸出系數(shù)關(guān)系

(2)任務(wù)及仿真參數(shù)

設(shè)經(jīng)濟(jì)航速veco=18 kn,滿功率最大航速vmax=25 kn,航程松弛系數(shù)?=0.8,任務(wù)時(shí)間T0=50 h,仿真次數(shù)N=105,失效偏倚系數(shù)x=0.5。

4.5 仿真結(jié)果分析

(1)仿真模型有效性檢驗(yàn)

本文通過解析方法對照驗(yàn)證仿真模型的有效性。由于本文研究的問題具有一定的復(fù)雜性,所以解析法建模和計(jì)算的成本會很高,因此對條件進(jìn)行一定的簡化,使得解析法過程不過于復(fù)雜。簡化假設(shè)條件并不影響仿真模型的有效性。一是令狀態(tài)轉(zhuǎn)移時(shí)間服從指數(shù)分布,即令表2中所有形狀參數(shù)b=1;二是令系統(tǒng)中儲備方式均為熱儲備，其余參數(shù)均不改變。

解析法只能求得系統(tǒng)的可靠性指標(biāo)(可用度),所以只能通過檢驗(yàn)仿真法中系統(tǒng)可用度是否與解析法一致的方式驗(yàn)證其有效性。由于在檢驗(yàn)仿真法中,系統(tǒng)可用度和任務(wù)成功概率為同時(shí)輸出,且任務(wù)成功性實(shí)則是系統(tǒng)可靠性在任務(wù)方面的表現(xiàn)。因此也能夠驗(yàn)證任務(wù)成功性指標(biāo)的有效性。

本文采用馬爾可夫方程與通用生成函數(shù)相結(jié)合的解析方法求解系統(tǒng)可用度。在圖9展示的結(jié)果中,仿真法下可用度曲線與解析法下基本一致,從而驗(yàn)證了仿真方法的有效性。

圖9 仿真法與解析法下的系統(tǒng)可用度對比

(2)抽樣方法效率對比

表4和圖10展示了仿真次數(shù)均為105次時(shí),使用方差減小方法(受迫轉(zhuǎn)移&失效偏倚)和僅使用蒙特卡羅間接仿真法的任務(wù)成功概率仿真結(jié)果。在系統(tǒng)存在冗余設(shè)計(jì)和任務(wù)時(shí)間較短的高可靠性系統(tǒng)背景下,方差減小方法的優(yōu)勢顯著。該方法在一次仿真中能夠抽得更多的故障樣本，在總體仿真中能夠抽得更多的任務(wù)失敗的樣本,因此任務(wù)成功性指標(biāo)收斂速度明顯更快,如圖10所示。方差減小方法所得結(jié)果的精確小數(shù)位數(shù)不受仿真次數(shù)的限制,105次仿真結(jié)果的精確小數(shù)位數(shù)不止5位。

表4 不同抽樣方法下任務(wù)成功性仿真結(jié)果(T0=50 h)

圖10 任務(wù)成功性指標(biāo)收斂速度對比

圖11比較了相同仿真次數(shù)N=105情況下蒙特卡羅間接仿真法與方差減小方法的系統(tǒng)瞬時(shí)可用度,并以樣本方差展現(xiàn)了瞬時(shí)可用度估值的波動程度。與蒙特卡羅間接仿真法相比,圖11中方差減小方法得到的可用度曲線波動幅度和方差范圍減小的效果顯著。兩種方法得到的仿真結(jié)果趨勢基本吻合,而方差減小方法的結(jié)果精度更高。這說明方差減小方法能夠有效結(jié)合本文所設(shè)計(jì)的隨機(jī)抽樣數(shù)值算法,達(dá)到了非齊次馬爾可夫系統(tǒng)加速仿真進(jìn)程和減小方差的目的。

圖11 系統(tǒng)瞬時(shí)可用度(T0=50 h)

(3)任務(wù)過程的仿真結(jié)果

圖12和圖13選取了兩次較為典型的仿真過程,能夠體現(xiàn)系統(tǒng)多狀態(tài)特性。并且能夠明顯地展示出系統(tǒng)狀態(tài)、系統(tǒng)的任務(wù)表現(xiàn)與任務(wù)目標(biāo)之間的實(shí)時(shí)關(guān)系。計(jì)劃航速表示目標(biāo)完成任務(wù)所需要達(dá)到的航速。系統(tǒng)功率輸出系數(shù)表示當(dāng)前艦船狀態(tài)下的性能輸出能力,其代表的艦船狀態(tài)將直接影響實(shí)際航速,從而影響該艦?zāi)芊裢瓿扇蝿?wù)目標(biāo)。需要指出的是,功率輸出系數(shù)表明系統(tǒng)輸出的能力高低,該系數(shù)為1時(shí),艦船可通過動力及傳動裝置在經(jīng)濟(jì)航速至最大航速之間的航速航行。實(shí)際航速將依據(jù)計(jì)劃航速的需要而定,其規(guī)則如式(36)所示。在圖12所示的任務(wù)過程中,該艦先后經(jīng)歷了2次宕機(jī)和1次降級退化,即便在25 h時(shí)已修復(fù)為完好狀態(tài),但是此時(shí)已無法滿足vp>vmax條件,無法按時(shí)到達(dá)指定海域,因此認(rèn)定為航速超限,任務(wù)失敗。圖13所示的過程與圖12類似,但宕機(jī)、降級所延誤的時(shí)間相對較短,且此后一直保持完好的性能狀態(tài),最終任務(wù)成功。

圖12 任務(wù)失敗仿真過程

圖13 任務(wù)成功仿真過程

為了比較計(jì)劃航速、實(shí)際航速和系統(tǒng)功率輸出系數(shù)變化趨勢之間的關(guān)系,圖14展示了上述指標(biāo)值的樣本均值曲線。需要強(qiáng)調(diào)的是,圖14中各指標(biāo)值的每個(gè)數(shù)據(jù)點(diǎn)并非對應(yīng)相應(yīng)時(shí)刻單個(gè)樣本的實(shí)際數(shù)據(jù),而是統(tǒng)計(jì)意義下的平均水平。盡管數(shù)據(jù)的變化幅度較小,但仍能展現(xiàn)各指標(biāo)變化趨勢之間的相互關(guān)系。隨著任務(wù)推進(jìn)、系統(tǒng)性能不斷退化,系統(tǒng)宕機(jī)、降級概率逐漸增加,平均功率輸出系數(shù)呈下降趨勢。因此,需要更高的計(jì)劃航速才能完成航渡任務(wù),表現(xiàn)為平均計(jì)劃航速隨著任務(wù)時(shí)間的推進(jìn)而上升。實(shí)際航速在滿足功率輸出條件的情況下,應(yīng)不低于計(jì)劃航速才能完成航渡任務(wù)。因此,實(shí)際航速隨著計(jì)劃航速的增長而提升。但幅度低于計(jì)劃航速的增幅,這是因?yàn)橄到y(tǒng)狀態(tài)的劣化導(dǎo)致任務(wù)后期實(shí)際航速可能無法達(dá)到計(jì)劃航速的需求。

圖14 實(shí)際航速、計(jì)劃航速、功率輸出系數(shù)的均值對比

(4)狀態(tài)轉(zhuǎn)移分布的影響

按照第4.4節(jié)中的任務(wù)參數(shù)設(shè)置,對非指數(shù)分布與指數(shù)分布狀態(tài)轉(zhuǎn)移下的仿真結(jié)果進(jìn)行對比。為保證仿真具有可比性,應(yīng)使不同狀態(tài)轉(zhuǎn)移分布具有相同期望值。根據(jù)表2中威布爾分布的各項(xiàng)參數(shù),計(jì)算各類轉(zhuǎn)移時(shí)間的期望值作為指數(shù)分布條件下的平均故障、維修時(shí)間,可表示為

E(T)=ηΓ(1+1/b)

(39)

式中:Γ(·)為伽馬函數(shù);E(T)為威布爾分布的期望值。任務(wù)成功性與可靠性結(jié)果如表5和圖15(a)所示,可見指數(shù)分布下的任務(wù)成功概率和可用度相對較低,這是由于在系統(tǒng)運(yùn)行之初指數(shù)分布的失效率更高,這種現(xiàn)象在任務(wù)時(shí)間較短的情況下較為明顯。說明若按照指數(shù)分布進(jìn)行建?？赡軙凸老到y(tǒng)可靠性及任務(wù)成功性。在本算例參數(shù)設(shè)置下,當(dāng)任務(wù)時(shí)間較長時(shí),系統(tǒng)中可修部分對可用度的影響已趨于穩(wěn)定?？捎枚鹊淖兓饕懿豢尚薏糠?致命故障)影響,由于不同狀態(tài)轉(zhuǎn)移分布下致命故障均視為指數(shù)分布,因此可用度的變化趨勢基本相同,但在相同時(shí)刻指數(shù)分布下的估值相對較低,如圖15(b)所示。

表5 不同狀態(tài)轉(zhuǎn)移概率分布下任務(wù)成功性仿真結(jié)果(T0=50 h)

圖15 不同狀態(tài)轉(zhuǎn)移概率分布下的系統(tǒng)瞬時(shí)可用度

5 結(jié) 論

本文研究了考慮時(shí)間相關(guān)狀態(tài)轉(zhuǎn)移的系統(tǒng)任務(wù)可靠性和任務(wù)成功性問題。針對非齊次馬爾可夫過程在離散事件仿真中的剩余分布抽樣,設(shè)計(jì)了服從非指數(shù)分布狀態(tài)轉(zhuǎn)移隨機(jī)抽樣的數(shù)值算法,并結(jié)合受迫轉(zhuǎn)移、失效偏倚等方差減小方法,改善了抽樣效率和估值精度。以艦船航渡任務(wù)為算例,評估了其任務(wù)可靠性和任務(wù)成功率,結(jié)果表明本文建立的仿真模型可以反映系統(tǒng)狀態(tài)、任務(wù)目標(biāo)和系統(tǒng)性能表現(xiàn)間的變化關(guān)系,驗(yàn)證了本文方法對評估復(fù)雜系統(tǒng)行為和故障機(jī)制條件下系統(tǒng)可靠性及任務(wù)成功性的合理性與有效性。本文方法不僅適用于馬爾可夫系統(tǒng)和具有解析形式狀態(tài)轉(zhuǎn)移分布的非齊次馬爾可夫系統(tǒng),還可適用于其他復(fù)雜系統(tǒng)行為及任務(wù)剖面組合的情況,具有一定的通用性。對于競爭性故障、級聯(lián)故障、共因失效等更為復(fù)雜的相關(guān)故障機(jī)制,以及等待維修狀態(tài)、維修資源可用性等維修方面問題,將在進(jìn)一步研究中加以考慮。