考慮時間相關(guān)故障的多狀態(tài)系統(tǒng)可靠性與任務(wù)成功性仿真評估方法

楊皓潔, 呂建偉, 徐一帆

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

0 引 言

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

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

為了在考慮時間相關(guān)故障、多狀態(tài)性能輸出、系統(tǒng)冗余儲備等因素的同時,改進抽樣效率和估值精度,本文提出了一種基于蒙特卡羅離散事件仿真,冗余多狀態(tài)系統(tǒng)可靠性及任務(wù)成功性的評估方法,用以針對故障率具有時變特點的時間相關(guān)故障。通過設(shè)計不服從指數(shù)分布狀態(tài)轉(zhuǎn)移隨機抽樣的數(shù)值算法,實現(xiàn)故障事件的隨機抽樣。在此基礎(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ù)過程中的某一時刻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)空間。實際上,系統(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)和運行方式確定的構(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ū)別在于一般故障為可修復故障,致命故障為無法修復或任務(wù)時間內(nèi)不具有維修條件的故障,也標為吸收態(tài)。λi,i′和μi,i′分別為該設(shè)備從狀態(tài)i′→i的故障率和修復率。

為達到提高系統(tǒng)可靠性的目的,系統(tǒng)中存在一些冗余結(jié)構(gòu),包括冷、溫、熱儲備。由冗余結(jié)構(gòu)和系統(tǒng)運行產(chǎn)生設(shè)備的工作狀態(tài)大致可分為3類:工作、儲備待機和維修。設(shè)備工作狀態(tài)根據(jù)系統(tǒng)結(jié)構(gòu)、備用設(shè)備啟動規(guī)則和設(shè)備當前狀態(tài)來確定。假設(shè)x、y分別為一個冗余結(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時刻的性能狀態(tài),數(shù)值為圖1中性能狀態(tài)的編號;Sysstate表示系統(tǒng)狀態(tài);down表示系統(tǒng)宕機,available表示系統(tǒng)可用。當主設(shè)備性能狀態(tài)低于備用設(shè)備或發(fā)生一般故障時應(yīng)當轉(zhuǎn)入維修,同時備用設(shè)備轉(zhuǎn)入工作狀態(tài);主設(shè)備維修完成后立即轉(zhuǎn)入工作狀態(tài),同時備用設(shè)備轉(zhuǎn)入待機狀態(tài);當系統(tǒng)中其他部分性能退化導致系統(tǒng)宕機時,主、備單元均轉(zhuǎn)入待機狀態(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)移速率,也是蒙特卡羅隨機事件抽樣的前提。

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

2.1 仿真流程設(shè)計

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

步驟 1在系統(tǒng)最底層(設(shè)備級)進行系統(tǒng)狀態(tài)轉(zhuǎn)移隨機抽樣,得到故障或維修事件發(fā)生的時間和轉(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)計結(jié)果。設(shè)任務(wù)成功概率為D,系統(tǒng)瞬時可用度為A(t),則有

D=(N-Nd)/N

(5)

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

(6)

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

Nd←Nd+wT0

(7)

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

Nf(t)←Nf(t)+wt

(8)

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

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

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

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

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

(10)

(11)

(12)

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

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

wt=wt′

(13)

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

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

(14)

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

(15)

Ok′=Λk′∪Ψk′

(16)

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

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

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

(17)

式中：ζ為0～1之間均勻分布的隨機數(shù)。當ζ落在式(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)。蒙特卡羅間接仿真法利用兩個隨機數(shù)就可以產(chǎn)生一次狀態(tài)轉(zhuǎn)移,而且不需要比較等步驟,一定程度上提高了計算效率,在此基礎(chǔ)上可使用方差減小方法進行加速抽樣。

2.3 方差減小方法

(18)

(19)

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

(20)

(21)

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

(22)

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

(23)

(2)失效偏倚

(24)

(25)

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

wt=

(26)

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

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

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

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

(27)

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

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

(28)

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

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

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

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

(29)

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

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

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

(30)

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

(31)

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

(32)

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

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

4 算例分析

4.1 算例背景

本算例以艦船及其裝備系統(tǒng)為對象進行分析,由于篇幅所限,這里僅以艦船動力、電力系統(tǒng)在航渡任務(wù)階段的任務(wù)可靠性和任務(wù)成功率評估為例,其他系統(tǒng)及任務(wù)的相應(yīng)計算也同樣適用。設(shè)某艦動力、電力系統(tǒng)的系統(tǒng)結(jié)構(gòu)組成如圖6所示。相關(guān)設(shè)備的性能狀態(tài)轉(zhuǎn)移關(guān)系如圖7所示。屬于第I類設(shè)備的有柴油機、燃氣輪機,這類設(shè)備有多種性能輸出表現(xiàn),且可能發(fā)生致命故障;屬于第II類設(shè)備的有柴油發(fā)電機組、軸系、螺旋槳,這類設(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)(柴油機1,柴油機2和燃氣輪機)和熱儲備結(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ù)目標[29]為:① 按時抵達目標海域;② 當航渡任務(wù)結(jié)束時系統(tǒng)仍處于可工作(非故障)狀態(tài)。目標①是對完成任務(wù)時間和空間的要求,指在規(guī)定的任務(wù)時間T0里完成規(guī)定的任務(wù)航程S0。艦船既不可提前到達,也不可晚于規(guī)定時間到達。提前到達可能導致本艦脫離協(xié)同兵力、提前暴露行蹤;晚于預(yù)定時間到達則難以實現(xiàn)兵力協(xié)同、完成既定作戰(zhàn)任務(wù)。目標②是對系統(tǒng)下一階段完成作戰(zhàn)任務(wù)的要求。總之,當且僅當以上兩項目標均滿足時,方可判定任務(wù)成功(實際情況下的任務(wù)成功性問題判據(jù)可能更為復雜,此處略去)。據(jù)此可知,任務(wù)目標是否達成應(yīng)結(jié)合任務(wù)成功性判定規(guī)則,根據(jù)特定的依據(jù)給予評判,如表1所示。在表1中,航速超限是指由于經(jīng)歷系統(tǒng)降級或故障維修,導致該艦在剩余航程中即使以最大航速也不能完成任務(wù)。維修超時是指任務(wù)結(jié)束時由于設(shè)備故障正在維修,導致該艦總體處于宕機、不可用的狀態(tài)。當發(fā)生維修超時時,艦船通常還未到達目標區(qū)域。致命故障是指由于設(shè)備致命故障導致該艦總體宕機且無法修復,從而無法完成剩余任務(wù)。

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

4.2 模型假設(shè)條件

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

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

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

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

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

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

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

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

(1)航速超限

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

t=t′+Δt

(33)

St=St′+vΔt

(34)

(35)

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

(36)

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

(2)維修超時

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

(3)致命故障

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

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

Dd=Nd/N

(37)

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

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

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

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

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

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

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

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

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

4.5 仿真結(jié)果分析

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

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

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

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

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

(2)抽樣方法效率對比

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(39)

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

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

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

5 結(jié) 論

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