應(yīng)用分層模型進(jìn)行核電站設(shè)備可靠性參數(shù)估計(jì)

陳 妍，何 亮，余少青，陳海英，張春明

(環(huán)境保護(hù)部核與輻射安全中心，北京100082)

?

應(yīng)用分層模型進(jìn)行核電站設(shè)備可靠性參數(shù)估計(jì)

陳 妍，何 亮，余少青，陳海英，張春明

(環(huán)境保護(hù)部核與輻射安全中心，北京100082)

目前國(guó)內(nèi)對(duì)于從多個(gè)核電站統(tǒng)計(jì)的設(shè)備失效數(shù)據(jù)進(jìn)行各特定核電站設(shè)備可靠性參數(shù)估計(jì)的方法研究尚少。本文研究了用于可靠性參數(shù)估計(jì)的分層模型以及實(shí)現(xiàn)分層模型的兩種方法：帶Kass-Steffey修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法和馬氏鏈蒙特卡洛方法。以設(shè)備需求失效的稀少失效數(shù)據(jù)樣本為例，推導(dǎo)了帶Kass-Steffey修正的Beta-Binomial模型原理并編程求解，研究了馬氏鏈蒙特卡羅方法及軟件計(jì)算，對(duì)比了核電站后驗(yàn)失效概率的計(jì)算結(jié)果。計(jì)算表明：兩種方法得到的部分失效概率后驗(yàn)估計(jì)的均值相差0～25%；95分位值相差5%～15%，兩種方法都可用于稀少數(shù)據(jù)的樣本估計(jì)。

參數(shù)經(jīng)驗(yàn)貝葉斯方法； Kass-Steffey修正；MCMC；WinBUGS

設(shè)備可靠性參數(shù)是核電廠設(shè)備可靠性管理的重要內(nèi)容，其在設(shè)備管理、老化管理、備品備件、維修規(guī)則等領(lǐng)域廣泛應(yīng)用。目前國(guó)內(nèi)對(duì)于如何從多個(gè)核電站統(tǒng)計(jì)的失效數(shù)據(jù)直接得到各特定核電站設(shè)備可靠性參數(shù)的計(jì)算方法研究尚少[1]。隨著我國(guó)《中國(guó)核電廠設(shè)備可靠性數(shù)據(jù)報(bào)告》(2015年)的發(fā)布[2]，已經(jīng)具備采集并統(tǒng)計(jì)我國(guó)各核電站設(shè)備失效的數(shù)據(jù)樣本的能力，亟須進(jìn)一步完善有關(guān)利用統(tǒng)計(jì)的失效數(shù)據(jù)計(jì)算可靠性參數(shù)的貝葉斯方法，以得到我國(guó)的特定電站設(shè)備的可靠性參數(shù)。

本文研究了可靠性參數(shù)估計(jì)通常采用的分層模型原理[3，4]，以及實(shí)現(xiàn)分層模型的兩種方法：帶Kass-Steffey修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法和馬氏鏈蒙特卡羅方法。這兩種方法都可以給出核行業(yè)級(jí)的可靠性參數(shù)，也可以計(jì)算電站級(jí)的可靠性參數(shù)。鑒于核電站的失效數(shù)據(jù)多是稀少失效數(shù)據(jù)樣本，本文以設(shè)備需求失效的稀少失效數(shù)據(jù)樣本參數(shù)估計(jì)為例，詳細(xì)推導(dǎo)了帶Kass-Steffey 修正的Beta-Binomial模型的數(shù)學(xué)原理并編程計(jì)算，同時(shí)研究了馬氏鏈蒙特卡羅方法及軟件計(jì)算，最后對(duì)比了兩種方法計(jì)算的特定電廠失效概率的后驗(yàn)估計(jì)值。

1 分層模型(Hierarchical Model)

m個(gè)核電站的某類失效數(shù)據(jù)若在部件類型、邊界定義、運(yùn)行狀態(tài)、運(yùn)行環(huán)境等方面相似，則失效數(shù)據(jù)xi可組成樣本x=(x1，…，xm)，各個(gè)核電站失效數(shù)據(jù)的待估參數(shù)通常采用2層的分層模型進(jìn)行估計(jì)[4]。

第1層：由于失效數(shù)據(jù)的相似性，將m個(gè)核電站的失效數(shù)據(jù)看成一個(gè)整體。盡管各個(gè)電廠的待估參數(shù)不同，但各個(gè)待估參數(shù)服從總體分布g。通過調(diào)用g分布m次，獨(dú)立產(chǎn)生各個(gè)核電站參數(shù)。將g分布作為先驗(yàn)分布，g分布中的參數(shù)λ稱為超參數(shù)。

第2層：依據(jù)各個(gè)核電站參數(shù)及失效數(shù)據(jù)分布，產(chǎn)生各個(gè)核電站的失效數(shù)據(jù)xi。

分層模型通過反復(fù)迭代，最后給出超參數(shù)λ和待估參數(shù)的后驗(yàn)估計(jì)值。

2 帶Kass-Steffey修正的Beta-Binomial模型

Beta-Binomial模型用于描述設(shè)備的需求失效，例如閥門的開啟/關(guān)閉失效、泵的啟動(dòng)失效、開關(guān)的打開/關(guān)閉失效、斷路器的斷開失效等。設(shè)第i個(gè)核電站某設(shè)備的失效次數(shù)為xi，需求次數(shù)為ni，失效概率pi，i=1, …,m，由于設(shè)備在每次需求時(shí)只有成功或失效2個(gè)結(jié)果，則xi服從參數(shù)pi的二項(xiàng)分布binomial(ni，pi)。

帶Kass-Steffey修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法中，第一層：首先依據(jù)m個(gè)核電站的失效數(shù)據(jù)xi，采用最大似然法估計(jì)g分布的超參數(shù)。第二層：將g分布作為先驗(yàn)分布，使用貝葉斯方法，結(jié)合單個(gè)電廠數(shù)據(jù)，得到特定電廠參數(shù)pi的后驗(yàn)估計(jì)值[5, 6]。

2.1 Beta-Binomial模型的超參數(shù)求解

鑒于失效概率pi服從二項(xiàng)分布的共軛分布，因此總體分布g(p|α,β)服從beta(α,β)分布，即g的概率密度函數(shù)為：

(1)

因此，x的無條件分布服從Beta-Binomial分布，即：

(2)

m個(gè)核電站的似然函數(shù)為xi的無條件分布函數(shù)之積：

(3)

因此，似然函數(shù)的對(duì)數(shù)為：

(4)

其中，(4)式利用了伽馬函數(shù)性質(zhì)Γ(a)=(a-1)Γ(a-1)。為了便于似然函數(shù)的求導(dǎo)計(jì)算，需進(jìn)行參數(shù)變換如下：

(5)

則(4)式變換為：

(6)

似然函數(shù)對(duì)數(shù)的偏導(dǎo)數(shù)方程組為：

(7)

超參數(shù)的估計(jì)值通過求解似然函數(shù)對(duì)數(shù)的偏導(dǎo)數(shù)方程組得到，即：

(8)

由于(8)式?jīng)]有解析解，通常的算法[7]：先設(shè)置δ初始值，利用False position求根算法，求解方程組第一式中的μ值。然后將μ值帶入方程組第二式，利用二分法再調(diào)整δ值。調(diào)整后的δ值再次作為初始值，帶入第一式求解μ。如此反復(fù)迭代，直到方程組都收斂，最后給出μ和δ。

2.2 Beta-Binomial模型的特定電站參數(shù)pi估計(jì)

一階修正下特定電廠參數(shù)pi的后驗(yàn)期望值不變[4, 5]，因此pi的后驗(yàn)期望值為：

(9)

帶Kass-Steffey 修正的pi后驗(yàn)方差估計(jì)為：

(10)

其中，

(11)

同時(shí)，由于二階協(xié)方差矩陣等于Fisher矩陣的逆，則:

var(μ)=1/J11

var(α)=1/J22

cov(μ,δ)=-J12/D

(12)

Fisher矩陣元為：

(13)

可以看到，帶Kass-Steffey修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法是通過最大似然函數(shù)求導(dǎo)數(shù)解得超參數(shù)，然后利用超參數(shù)和特定電廠失效樣本xi，計(jì)算特定電廠參數(shù)pi。值得注意的是，帶Kass-Steffey修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法中的總體分布g選取的是二項(xiàng)分布的共軛分布。

3 馬氏鏈蒙特卡羅方法

3.1 分層貝葉斯方法及MCMC

分層貝葉斯方法是完全的貝葉斯方法，不同于先計(jì)算超參數(shù)，對(duì)于需求失效數(shù)據(jù)，分層貝葉斯方法是將2個(gè)超參數(shù)和各個(gè)電廠待估參數(shù)組成維數(shù)為2+m維的先驗(yàn)參數(shù)矢量θ。θ的先驗(yàn)分布先假定：包括超參數(shù)α、β的超先驗(yàn)分布，m個(gè)核電站的pi條件分布及總體分布g，g分布可以假設(shè)為共軛分布，也可以是非共軛分布。

馬氏鏈蒙特卡羅方法(MarkovChainMonteCarlo，MCMC)提供了一種計(jì)算分層貝葉斯模型的方法，其利用計(jì)算機(jī)模擬從后驗(yàn)分布抽樣然后用樣本估計(jì)值推斷參數(shù)。Gibbs抽樣是一種特殊的馬氏鏈蒙特卡羅算法，具體是從某個(gè)初始點(diǎn)出發(fā)，通過全條件分布的循環(huán)抽樣產(chǎn)生馬氏鏈，算法如下[4]：

3.2 BUGS的Doodle模型的建立

BUGS(Bayesian inference Using Gibbs Sampling)是采用MCMC算法(利用Gibbs抽樣)解決貝葉斯估計(jì)的軟件。在WinBUGS軟件中根據(jù)(1)式，在Doodle中建立相應(yīng)的節(jié)點(diǎn)模型，如圖1所示[8]。

圖1 需求失效的Doodle模型Fig.1 The Doodle model of failures on demand

4 實(shí)例計(jì)算及結(jié)果比較

本文以68個(gè)核電站輔助給水系統(tǒng)電動(dòng)部分啟動(dòng)失效的參數(shù)計(jì)算為例[4](累計(jì)失效個(gè)數(shù)為6個(gè))，樣本數(shù)據(jù)見表1，分別用Kass-Steffey 修正的Beta-Binomial模型和WinBUGS的方法進(jìn)行計(jì)算，比較計(jì)算結(jié)果及比較見表2。

表1 AFW系統(tǒng)電動(dòng)部分啟動(dòng)失效數(shù)據(jù)

表2 部分核電站后驗(yàn)需求失效概率的兩種方法計(jì)算結(jié)果比較

在WinBUGS的Gibbs抽樣模型中，假設(shè)超參數(shù)α服從均值和標(biāo)準(zhǔn)差都為1的指數(shù)超先驗(yàn)分布，β服從gamma(1.0, 0.035)超先驗(yàn)分布，p(i)為第i個(gè)核電站輔助給水系統(tǒng)電動(dòng)部分啟動(dòng)失效概率，n(i)描述需求次數(shù)，x(i)描述失效次數(shù)，i=1,, 68。在文本中輸入失效數(shù)據(jù)及需求次數(shù)數(shù)據(jù)，編譯后加載超參數(shù)的初始值。在Sample Monitor Tool的節(jié)點(diǎn)中設(shè)置監(jiān)測(cè)變量alpha、beta及p[i]，在Update中設(shè)置迭代次數(shù)。迭代完成后，在node中輸入alpha、beta，查看stats獲得參數(shù)估計(jì)結(jié)果。計(jì)算先驗(yàn)分布超參數(shù)計(jì)算的部分核電站后驗(yàn)需求失效概率見表2第3列。需要注意的是，WinBUGS軟件使用時(shí)需要判斷馬氏鏈的收斂性、設(shè)置合理的迭代次數(shù)以及合適的超參數(shù)的超先驗(yàn)分布。

選取具有代表性的6個(gè)核電站(有失效數(shù)據(jù)次數(shù)的5個(gè)核電站和1個(gè)有0失效次數(shù)的核電站)，比較兩種方法得到的后驗(yàn)估計(jì)值相對(duì)變化量見表2第4列。可以看到，兩種方法得到的失效概率后驗(yàn)估計(jì)的均值和95分位值相差較小，均值相差0～25%，95分位值相差5%～15%；5分位值相差大，達(dá)32%～348%。

可以看到，對(duì)于后驗(yàn)參數(shù)的分位值估計(jì)，分位值越低，對(duì)模型假設(shè)的依賴越強(qiáng)，這一結(jié)論對(duì)于稀少樣本的參數(shù)估計(jì)尤為明顯，不同模型的計(jì)算結(jié)果差異越大。然而，在風(fēng)險(xiǎn)評(píng)價(jià)中，更多關(guān)注均值和95分位值(均值常作為點(diǎn)估計(jì)值，95分位值常作為風(fēng)險(xiǎn)估計(jì)的上限)，不同模型計(jì)算的低風(fēng)險(xiǎn)值(5分位值)的差異一定程度上可以不考慮[4]，這也是通常的可靠性參數(shù)采用均值和誤差因子，而誤差因子(Error Factor)通常采用95分位值除以均值計(jì)算的原因。因此，帶Kass-Steffey修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法，也通常可以作為完全貝葉斯方法——馬氏鏈蒙特卡羅方法的一階近似，廣泛應(yīng)用于美國(guó)NRC的核電站的可靠性參數(shù)估計(jì)中[3]。此外，采用帶Kass-Steffey修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法可以避免馬氏鏈蒙特卡羅方法在收斂性判斷、迭代次數(shù)、超先驗(yàn)分布設(shè)置合理性等方面的問題。

5 結(jié)論

本文以設(shè)備需求失效的稀少失效數(shù)據(jù)樣本為例，分別用帶Kass-Steffey修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法和馬氏鏈蒙特卡羅方法計(jì)算了特定核電廠失效概率的后驗(yàn)估計(jì)值。計(jì)算表明：兩種方法得到的部分核電廠失效概率后驗(yàn)估計(jì)的均值相差0～25%，95分位值相差5%～15%。對(duì)于大量失效數(shù)據(jù)的可靠性參數(shù)計(jì)算，兩種方法的計(jì)算結(jié)果更為接近。一定程度上，帶Kass-Steffey修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法可以作為分層貝葉斯馬氏鏈蒙特卡洛方法的一階近似。需要注意的是，帶Kass-Steffey 修正的參數(shù)經(jīng)驗(yàn)貝葉斯方法由于沒有解析解，需要進(jìn)行數(shù)值方法求解，有時(shí)可能因?yàn)閿?shù)據(jù)的特殊性，存在無法求解的可能，且總體分布選取二項(xiàng)分布的共軛分布。馬氏鏈蒙特卡羅方法在可靠性參數(shù)計(jì)算時(shí)則需要討論收斂性判斷、超先驗(yàn)分布設(shè)置的合理性、迭代次數(shù)等問題。這些問題都需要進(jìn)一步深入研究。

[1] 茆定遠(yuǎn), 薛大知.核電站PSA分析中可靠性數(shù)據(jù)處理的貝葉斯方法[J].核動(dòng)力工程, 2000, 21(5):451-455.[2] 國(guó)家核安全局.中國(guó)核電廠設(shè)備可靠性數(shù)據(jù)報(bào)告，2015.

[3] Eide S A, Wierman T E, Gentillon C D, et al.Industry-Average Performance for Components and Initiating Events at U.S.Commercial Nuclear Power Plants [R].NUREG/CR-6928.2007.

[4] Atwood C L, LaChance J L, Martz H F, et al.Handbook of Parameter Estimation for Probabilistic Risk Assessment[R].NUREG/CR-6823.2003.

[5] Robert E.Kass, Duane Steffey, Approximate Bayesian Inference in Conditionally Independent Hierarchical Models (Parametric Empircal Bayes Models) [J].Journal of the American Statistical Association, 1989, 84(407)717-726.

[6] 陳妍, 鄭鵬, 李朝君, 等.PSA 分析中可靠性參數(shù)的Kass-Steffey修正原理及應(yīng)用研究[J].核動(dòng)力工程, 2015, 36(3):70-74.

[7] Atwood, C.L. .1994, Hits per Trial: Basic Analysis of Binomial Data, EGG-RAAM-1 1041, Idaho National Engineering and Environmental Laboratory, Idaho Falls, ID.

[8] 陳妍, 余少青.WinBUGS軟件在核電站可靠性參數(shù)估計(jì)中的應(yīng)用[J], 科學(xué)技術(shù)與工程.2016, 16(1):151-154.

Calculation of Reliability Parameters of Nuclear Power Plant with Hierarchical Model

CHEN Yan，HE Liang，YU Shao-qing，CHEN Hai-ying，ZHANG Chun-ming

(Nuclear and Radiation Safety Center, Ministry of Environmental Protection, Beijing 100082)

There is little research on the calculation methods of evaluation reliability parameters directly from our own statistics failure data.The parametric empirical Bayes model with Kass-Steffey adjustment and the MCMC method are two Bayesian methods to estimate the reliability parameters.Take the failures on demand with sparse data as an example, this paper derives the principle of Beta-Binomial with Kass-Steffey adjustment, adopts two methods to calculate the hyperparameter and plant-specific posterior parameter.It shows that posterior mean of NPPs change 0～25%, posterior 95thchange 5%～15%.Two methods are all suitable for estimating the reliability parameters from sparse failure data.

Parametric empirical Bayes models; Kass-Steffey adjustment; MCMC; WinBUGS

2017-04-11

國(guó)家科技重大專項(xiàng)(2013ZX06002001-008)項(xiàng)目資助

陳 妍(1982—)，女，江蘇徐州人，高級(jí)工程師，博士，理論物理，現(xiàn)從事核電站概率安全評(píng)價(jià)方面研究

何 亮：heliang@chinansc.cn

TL364

A

0258-0918(2017)03-0458-06