高精度相關變量隨機數序列產生方法

馬續波，李耀舟，吳 屈，彭星杰，黃自鋒，朱潤澤，張 斌

(1.華北電力大學 核科學與工程學院，北京 102206；2.中國核動力研究設計院 核反應堆系統設計技術重點實驗室，四川 成都 610213)

蒙特卡羅方法已廣泛應用到了各領域。由于具有已知分布的隨機變量抽樣方法在蒙特卡羅方法中占有重要位置，并已得到廣泛研究[1-4]。采用統計抽樣方法進行不確定度和靈敏度分析也成為了核工程領域研究熱點[5-7]，但采用統計抽樣法進行不確定度分析時，需大量抽樣樣本，而對于有些問題，問題求解本身需消耗大量時間，如求解玻爾茲曼輸運方程、熱傳導方程等，這使采用統計抽樣法進行不確定度分析需消耗大量時間，限制了統計抽樣法的應用和推廣。

采用統計抽樣法進行不確定度分析，首先要產生多維隨機變量序列，而對于多維隨機變量的產生，需根據多維隨機變量之間的協方差矩陣，產生具有相關變量的隨機數序列。文獻[5]給出了相關變量隨機數序列產生方法的理論推導及數值驗證。雖從理論上講，產生的相關變量隨機數序列的協方差矩陣與真實的協方差矩陣應一致，但研究發現，僅在大樣本情況下，產生的相關變量隨機數的協方差矩陣才與真實的協方差矩陣一致，因此，對于實際問題，這將導致耗費大量計算資源。文獻[8]研究利用最小樣本數量進行不確定度分析的方法，并提出了面向協方差矩陣的COST方法，但并未給出COST方法的理論基礎。為采用少樣本也能實現高精度多維隨機變量序列抽樣，本文首先證明COST方法的理論基礎，然后利用兩個變量的例題進行驗證分析，研究變量之間相關性的去除和施加的方法，最后進行工程上的驗證分析。

1 改進的相關隨機變量模擬理論

設Xm×n=(X1,X2,…,Xn)m為需要產生的m組n維相關正態分布隨機變量序列，其均值u=(u1,u2,…,un)，協方差矩陣為A。產生n維獨立標準正態分布隨機變量抽樣序列Ym×n=(Y1,Y2,…,Yn)，每個變量抽樣序列數為m。利用式(1)和獨立標準正態分布隨機變量抽樣Yn即可產生隨機變量X：

X=ΣY+u

(1)

式中，矩陣Σ需要通過對協方差矩陣A分解得到。

V[X]=A=ΣΣT

(2)

由式(2)可知，僅將協方差矩陣A分解成矩陣Σ和ΣT相乘的形式，就可利用式(1)來產生多維相關正態分布隨機變量序列[5]。

上述理論是較嚴格的數學推導，但上述理論中有一個假設，即n維獨立標準正態分布隨機變量抽樣序列Ym×n=(Y1,Y2,…,Yn)之間的協方差矩陣為單位矩陣I，在實際計算過程中，發現上述假設并不完全成立，特別是在抽樣樣本較少的情況下，n維獨立標準正態分布隨機變量抽樣序列之間的相關性非常大。文獻[8]假設這個協方差矩陣為I*，且此協方差矩陣可分解為：

I*=MMT

(3)

由于n維隨機變量正態分布的相關性傳遞到了隨后的抽樣樣本上，導致實際產生的抽樣樣本的協方差矩陣并不完全是A，而是矩陣A與協方差矩陣I*的混合，由此導致了產生的n維相關隨機變量的協方差矩陣與真實的協方差矩陣并不能滿足完全一致。為能使得產生的n維相關隨機變量的協方差矩陣完全一致，首先需產生真正的n維獨立標準正態分布，假設真正的n維獨立標準正態分布為Y′，則Y和Y′的關系為：

Y′=M-1Y

(4)

式(2)就變為：

V[X]=E[(X-E(X))(X-E(X))T]=

E[(ΣY′+u-u)(ΣY′+u-u)T]=

E[(ΣY′)(ΣY′)T]=E[ΣY′Y′TΣT]=

ΣV(Y′)ΣT=ΣΣT

(5)

改進后的方法與之前的方法不同之處在于，首先需先求出n維獨立標準正態分布隨機變量抽樣序列的協方差矩陣，然后利用式(3)和(4)將其轉換成真正的n維獨立標準正態分布隨機抽樣序列，在此基礎上把實際的相關性加上去，才能產生滿足協方差矩陣守恒的n維相關隨機變量序列。

本文提出的高精度相關變量隨機數序列產生方法的具體實現方法為：根據樣本總數N將概率區間等分為N個，使用RANDOM_NUMBER函數抽取隨機數ξ1，根據標準正態分布的分布函數與誤差函數的對應關系計算ξ1對應的概率，舍選抽樣直到保證N個等概率區間中每個區間均能抽取到有且僅有1個隨機數，則抽樣結束，即單變量拉丁超立方正態分布抽樣，重復上述抽樣可得到多序列抽樣樣本，然后利用式(4)進行去除相關變量之間的相關性，即可得到高精度相關變量隨機數序列。

2 測試和驗證

2.1 簡單例題測試和驗證

采用文獻[5]中的例題，要求產生隨機變量X1、X2、X3，其中X1、X2∈[-1,1]，X3∈[0,4]。3個隨機變量具有如下相關系數矩陣R：

(6)

利用Fortran2003語言分別開發了Cholesky和Jacobi分解法對矩陣R進行計算，假設X1、X2、X3的標準偏差分別為0.1、0.2、0.3，均值為1、2、3。利用上述方法產生的20個樣本列于表1。采用高精度相關變量隨機數序列方法(改進方法)與傳統的Jacobi方法和Cholesky方法的計算結果[5]亦列于表1，由表1可見，即使對于3個抽樣序列，傳統的Jacobi方法和Cholesky方法均需產生大量樣本(大于2萬)才能保證相關系數一致。而采用高精度相關變量隨機數序列方法，僅采用20個樣本，3個參數的相關系數即能保持與原相關系數一致(表2)。

表1 不同方法的相關系數對比Table 1 Comparison of correlation coefficient with different calculation methods

表2 利用改進方法產生的20個樣本Table 2 Twenty samples using high-precision correlated variable sampling method

2.2 矩陣分解

為研究該方法對多個隨機數序列的適用性，本文對238U的輻射俘獲反應道33群相關系數矩陣進行了研究。基于ENDF/B-Ⅶ.1核評價數據庫，采用國際著名核截面處理程序NJOY2016，制作了238U的輻射俘獲反應道33群相關系數矩陣，如圖1所示。由圖1可見，主對角線表示本能群與本能群的相關系數，所有能群的自相關系數均為1，并且不同能群之間具有較強的相關性。采用本文提出的高精度相關變量隨機數序列方法產生了33個抽樣序列，即n=33，每個抽樣序列樣本數為34，即m=34。研究發現，若m≤n，則式中的矩陣M無法求逆矩陣，因此，m的最小值為n+1。所以本文抽樣的每個抽樣序列的樣本數為34。圖2示出了利用抽樣得到的樣本序列計算得到的33個能群之間的相關系數矩陣，與圖1給出的原矩陣符合很好。圖3示出了抽樣樣本得到的相關系數矩陣(圖2)與原矩陣的相對誤差分布，最大相對誤差僅0.061%，絕大部分的相對誤差為0，由此說明了高精度相關變量隨機數序列方法具有較強的處理多變量的能力，且具有高精度的特點。

圖1 238U 33群俘獲截面的相關系數矩陣Fig.1 Correlation coefficient matrix of 238U capture with 33 neutron groups

圖2 利用34個相關變量樣本序列計算得到的238U俘獲反應道相關系數矩陣Fig.2 Calculated correlation coefficient matrix of 238U capture with 34 high-precision correlated variable pseudo-random number sequences

圖3 238U俘獲相關系數矩陣的相對誤差Fig.3 Relative error of correlation coefficient matrix of 238U capture

2.3 鉛基快堆LFR有效增殖因數不確定度

鉛基快堆作為世界上第四代先進核反應堆代表堆型之一，與傳統核反應堆相比有著獨特的優勢及特點，且有望最早實現商業化的第四代核反應堆[9]。本文以歐共體核能局(NEA)組織的國際評價合作工作組(簡稱NEAWPEC)提出的鉛基快堆LFR為研究對象[10]，研究主要核素不同反應道對有效增殖因數keff的不確定度的貢獻，分析了高精度相關變量隨機數序列產生方法工程應用的可行性。LFR的二維rz模型如圖4所示，該反應堆主要由堆芯、反射層、屏蔽層、軸向屏蔽層、控制棒等組成。采用鉛(Pb)作冷卻劑、U-TRU-Zr合金作燃料，堆芯中次錒系核素(MA)的含量為2.0%，反應堆滿功率為900 MWt，反應堆循環壽期為310 d。

圖4 鉛基快堆LFR的rz模型Fig.4 rz model of lead cooled fast reactor LFR

選取了對LFR的keff具有重要貢獻的6個核素多個反應道進行研究，對比高精度相關變量隨機數序列產生方法與傳統抽樣方法計算得到keff的不確定度的相對誤差。239Pu和238U均為LFR燃料中的重要核素，206Pb、207Pb和208Pb是本基準例題LFR冷卻劑的成分，56Fe為結構材料的重要組成部分，而燃料、冷卻劑和結構材料幾乎構成了LFR堆芯的全部材料，因此，這幾個核素具有很好的代表性。計算結果列于表3。表3中的參考解是采用直接擾動方法得到每個核素每個反應道的敏感性系數，然后再利用Sandwich公式計算得到。每個反應道的抽樣樣本總數為50，傳統正態分布抽樣方法的抽樣樣本為100。由表3可見，雖然傳統正態分布抽樣樣本總數較高精度相關變量隨機數序列產生方法的樣本總數高1倍，但其相對誤差最大可達12.5%，而高精度相關變量隨機數序列產生方法的最大相對誤差僅1.7%，計算精度有明顯提高。

表3 不同計算方法計算各核素相對于keff的不確定度對比Table 3 Comparison of uncertainty of keff with respect to different isotopes by different sampling methods

3 結論

當采用蒙特卡羅方法研究許多問題時，需要產生多維相關隨機數序列。但之前的研究表明：需特別大的抽樣樣本序列才能基本保證樣本的協方差矩陣與真實的協方差矩陣一致。由于采用大量抽樣樣本，導致很大的時間消耗。本文從理論上分析了以前的多維相關變量隨機數序列的協方差矩陣與真實的協方差矩陣有差別的原因是之前產生的多維相關變量隨機數序列，在考慮真實協方差矩陣前，相關變量隨機數序列之間并不能完全滿足不相關假設條件，通過對產生后的樣本的進一步處理，即可得到，在考慮真實協方差矩陣之前，即能滿足各變量之間完全滿足不相關條件。然后，利用數值方法對上述理論進行了驗證，通過3個變量的例題可見，原來的抽樣方法需上萬個抽樣樣本才能計算得到與原相關協方差矩陣一致的矩陣，而采用高精度相關變量隨機數序列產生方法，20個抽樣樣本即可完全重構出原協方差矩陣。樣本數量減小了5個量級，而計算精度有了提高。在此基礎上，通過對33群的238U輻射俘獲反應的相關系數矩陣進行分析，得到了產生的樣本數一定要大于樣本變量數，利用高精度相關變量隨機數序列產生方法，產生的樣本可高度還原原相關系數矩陣，最大相對誤差小于0.061%。最后，利用鉛基快堆LFR對該方法的工程應用性進行了研究，研究結果表明，雖然傳統正態分布抽樣樣本總數較高精度相關變量隨機數序列產生方法的樣本總數高1倍，但其相對誤差仍較大，最大達12.5%，而高精度相關變量隨機數序列產生方法的最大相對誤差僅1.7%，計算精度有明顯提高。采用此方法時，抽樣序列樣本數要大于變量的個數，即m>n，否則，因為矩陣M無法求逆矩陣，而無法采用該方法，m的最小值為n+1。通過鉛基快堆LFR驗證表明，該方法具有工程應用前景。