基于自適應MCMC和子集模擬的非能動系統熱工水力可靠性評估

蔣立志，蔡 琦，張永發

(海軍工程大學 核科學技術學院，湖北 武漢 430033)

與能動系統相比，非能動系統驅動力較小，設計參數、初始條件等因素波動帶來的不確定性以及對相關物理現象認識的不充分，均可能造成熱工水力過程失效、無法實現設計功能。對非能動系統進行熱工水力可靠性(TH-R)評估能從可靠性角度增強對非能動系統熱工水力過程不確定性的認識，相關分析結果可作為系統性能評估與優化的參考。最初非能動系統TH-R評估主要采用以蒙特卡羅(MC)為基礎的模擬方法[1-2]，需大量調用系統熱工水力模型進行計算，以獲取足夠的輸入、輸出數據，但以RELAP5為代表的熱工水力程序運行耗時較長，帶來極高的運算代價和較低的評估效率。重要抽樣(IS)方法的引入[3-4]一定程度緩解了運算代價的問題，但在高維問題中重要抽樣密度函數的構造十分困難且存在維數災難問題[5]。子集模擬(SS)方法[6]適用于高維小失效率模型，它的引入[7-8]在評估效率和穩健性之間帶來了一種均衡，其基本思想是將較小失效概率分解為較大失效概率的乘積并依次估計這些條件失效概率。SS中條件失效概率的估計需對各中間失效事件的條件概率密度函數進行精確采樣，通常采用修正Metropolis-Hastings(MMH)[9]算法實現。MMH算法中建議分布的選擇對算法的效率和失效率估計精度影響很大，最優建議分布的選擇較為困難，且隨失效率水平的降低算法性能會出現退化。文獻[10-11]討論了Metropolis-Hastings(MH)算法建議分布的最優問題和自適應馬爾可夫鏈蒙特卡羅(MCMC)的合理框架，但相關優化結論針對特定目標分布假設，并不完全適用于子集模擬的中間條件分布。文獻[12]基于標準正態空間內的條件采樣(CS)算法提出一種適用于SS的自適應MCMC方法，稱為自適應條件采樣(aCS)方法。該方法在當前鏈狀態和候選樣本間引入相關性，避免了建議分布選擇的復雜工作[13]，能在條件采樣過程中直接和靈活地調整算法參數，使候選樣本接受率在目標值附近保持穩定，以獲取更高的算法效率和失效率估計精度。

本文引入aCS-SS方法，提出一種基于自適應MCMC和SS的非能動系統TH-R評估方法，以某型核動力裝置二次側非能動余熱排出試驗系統為例進行驗證。

1 基本理論

1.1 SS的基本思想

假設系統輸入變量的聯合概率密度函數為fX(x)，則系統的失效率PF為：

(1)

式中：x=(x1,x2,…,xn)∈Rn為輸入參數向量，Rn為n維輸入參數空間；F為失效域；IF(x)為示性函數，x∈F時IF(x)=1，否則IF(x)=0。采用MC方法時系統失效率的估計值為：

(2)

(3)

式(3)中條件失效概率的估計值為：

(4)

式中，fX(x|Fj-1)=fX(x)IFj-1(x)/P(Fj-1)為x落入Fj-1內的條件概率密度函數。由于F0是一個必然事件，P1的計算可采用標準的MC模擬進行。當j≥2時Pj的估計需根據條件分布fX(x|Fj-1)進行采樣，高效的條件采樣是SS的核心。

1.2 基于MMH算法的SS

3) 重復上述步驟，直至樣本數量滿足要求。

1) 設定各層的樣本數N及條件失效概率p0，應該盡量保證p0N為正整數。

1.3 CS的基本原理

上述基于CS算法的MCMC策略其穩態分布漸進收斂于條件分布fX(x|Fj-1)，文獻[12]給出了詳細的證明過程。CS算法的效率和性能僅依賴于相關系數ρi的選擇，與MMH算法相比避免了建議分布形式、參數的選擇，避免了預候選樣本的接受/拒絕過程，算法形式更加簡潔、靈活，運算代價更低。但針對不同的系統類型和失效率水平，如何合理選擇相關系數ρi依然是一個困難的問題，較高和較低的相關系數均會導致條件樣本的相關性過高，降低算法效率和估計精度。

1.4 基于aCS的子集模擬

文獻[10-11]的相關研究表明，目標分布為一維正態分布時MH算法存在一近似的最優接受率a*≈0.44，文獻[15]的算例表明SS中MMH算法的最優接受率在0.3～0.5之間。CS算法中候選樣本的接受/拒絕僅依賴于系統極限狀態函數的響應，因此從條件分布fX(x|Fj-1)采樣可等效為一維條件采樣問題，文獻[12]提出一種aCS方法，能使候選樣本接受率較穩定地保持在目標值附近。利用aCS算法從條件分布fX(x|Fj-1)采樣的詳細流程如下。

2) 自適應迭代過程：令迭代次數r=1,2,…,Ns/Na，重復執行以下步驟。

(1) 根據標準差初值σ0i和當前尺度參數λr計算相關系數：

σi=min(λrσ0i,1)

(5)

(6)

(7)

(8)

2 算例分析及結果討論

本文以某型核動力裝置二次側非能動余熱排出試驗系統[16-18]作為算例驗證aCS-SS算法的有效性，系統原理如圖1所示。采用文獻[16]建立的熱工水力模型并選擇啟動工況5為運行工況進行分析，該模型使用RELAP5/Mod3.2作為建模工具。由于缺少充足的試驗數據和運行經驗，通過專家分析及工程判斷方法識別影響系統熱工水力行為的32個不確定性輸入參數，參數分布類型列于表1。由于試驗系統未模擬反應堆一回路，因此選擇關注的系統輸出為系統運行期間(20 000 s)冷凝器管側壓力峰值pMAX：若運行期間pMAX不大于管道所能承受的壓力上限值pm，則認為蒸汽發生器模擬體中的熱量能被成功排出至冷卻水箱，即g(x)=pm-pMAX<0時系統功能失效，否則系統功能正常。

圖1 非能動余熱排出系統原理簡圖Fig.1 Simplified figure of passive residual heat removal system

表1 不確定性參數分布類型Table 1 Distribution type of uncertainty parameter

注：N表示正態分布，U表示均勻分布

為比較不同失效率水平下的算法性能，設置4組壓力上限值pm1

1) MMH-SS：采用MMH算法產生條件樣本，建議分布選取為正態分布，其標準差根據種子樣本計算得到。

2) CS-SS：采用CS算法產生條件樣本，為分析不同相關系數對失效率估計精度的影響，各分量相關系數取3種配置：ρi=0.7、ρi=0.8、ρi=0.9，分別記為CS-SS-a、CS-SS-b和CS-SS-c。

3) aCS-SS-a：采用aCS算法產生條件樣本，標準差初值取σ0i=1，尺度參數初值取λ1=0.8。

4) aCS-SS-b：采用aCS算法產生條件樣本，標準差初值根據種子樣本計算得到，尺度參數初值取λ1=0.8。

由于a*≈0.44是在目標分布為一維正態分布的假設下獲得的結論，并不完全適用于本文所評估的系統。經過測試，系統失效率為1.94×10-3時MMH-SS算法的候選樣本接受率約為0.5，且失效率估計精度較高，因此將aCS-SS的目標接受率設定為a*=0.5。設定N=1 000、p0=0.1，在不同失效率水平下，利用4種策略各進行100次失效率估計，各策略下采用相同的100個隨機數種子。最終得到整個SS過程中候選樣本的平均接受率隨失效率水平的變化，如圖2所示。失效率估計平均誤差隨失效率水平的變化如圖3所示，失效率估計平均變異系數(SS變異系數的近似估計方法參考文獻[7,13])隨失效率水平的變化如圖4所示。

圖2 平均接受率隨失效率水平的變化Fig.2 Average acceptance rate plotted against varying failure probability level

圖3 失效率估計平均誤差隨失效率水平的變化Fig.3 Average error of failure probability estimate plotted against varying failure probability level

圖4 失效率估計平均變異系數隨失效率水平的變化Fig.4 Average coefficient of variation of failure probability estimate plotted against varying failure probability level

由圖2可知：隨系統失效率水平的不斷降低，4種估計方法的平均接受率均有所下降，但aCS-SS-a和aCS-SS-b能使接受率在目標值附近保持穩定水平；CS-SS的平均接受率與相關系數的取值有直接關系，相關系數越大平均接受率越高，反之亦然。結合圖3、4的結果可知：隨失效率水平的降低，4種方法的平均估計精度均降低、平均變異系數均增大，這是由于較低失效率水平下失效區域非常狹小，有效的條件采樣變得更加困難；失效率處于較高水平時4種方法精度差別不大，但失效率處于較低水平時aCS-SS-b的平均精度最高、平均變異系數最小；CS-SS算法中相關系數ρi取不同值時失效率平均估計精度和平均變異系數變化劇烈，對于本文算例ρi=0.8時失效率平均估計精度較高、變異系數較大，但在實際工程應用中對于不同的研究對象很難準確地找到一個比較合理的相關系數取值；與aCS-SS-a相比，aCS-SS-b的平均精度更高、變異系數更小，這是由于aCS-SS-b中標準差初值σ0i由種子樣本估計得到，這一過程本質上是對失效區域的一種自適應探索，當各輸入參數對系統性能的影響不相同時aCS-SS-b的性能必然高于aCS-SS-a，但當各輸入參數對系統性能的影響基本相同時二者的性能無明顯差別。

SS的失效率估計精度取決于各子集邊界的準確劃分和條件失效域的高效采樣。采樣過程中候選樣本接受率過低或過高，均會導致條件失效樣本間相關性較高、重復樣本較多，最終導致各子集邊界劃分出現偏差以及最后一層子集的失效率估計誤差增大。在失效率為2.90×10-6的水平下，分別利用MMH-SS、CS-SS-b、aCS-SS-a和aCS-SS-b進行1次失效率估計，第2～6層子集對應的候選樣本接受率變化列于表2。由表2的變化規律可看出：隨各子集的逐漸劃分，候選樣本接受率逐漸降低，這是由于后續子集對應的區域越來越狹小、采樣也越來越困難；MMH-SS和CS-SS算法中各子集接受率的降低較為劇烈，這是因為各子集的條件采樣使用的均是相同的策略或相關系數；aCS-SS中接受率的下降較為緩慢，即使在最后一層子集也能保持接近0.5的水平，這是由于aCS-SS的自適應調節過程可令各子集對應的相關系數逐漸增大，特別是aCS-SS-b中對失效區域的自適應探索過程還可使各輸入參數對應的相關系數隨參數重要度的不同自動調節。

表2 各子集候選樣本的接受率Table 2 Acceptance rate of candidate sample for each subset

3 小結

本文引入一種基于aCS的自適應MCMC算法，提出一種基于自適應MCMC和SS的非能動系統熱工水力可靠性評估方法，以某型核動力裝置二次側非能動余熱排出試驗系統為算例進行驗證，對4種估計方法的性能進行了比較。計算結果表明：aCS-SS避免了MMH-SS中建議分布形式和參數的選擇，算法形式更加簡潔、靈活，算法中的參數可自適應調整以保持穩健性；aCS-SS能更好地使候選樣本接受率在目標值附近保持穩定水平，在較低失效率水平(10-5、10-6)下aCS-SS的估計精度和穩健性均高于目前常用的MMH-SS算法。隨非能動技術的不斷發展、運行經驗的不斷積累以及各類非能動安全系統的不斷改進，反應堆非能動系統的物理過程失效會進一步降低，本文方法對于較低失效率水平下非能動系統TH-R評估問題是非常有前景的。