集成學習模型的堆芯物理關鍵參數預測

馬季， 郝琛， 謝曉芹， 生義

(1.哈爾濱工程大學 核安全與仿真技術國防重點學科實驗室，黑龍江 哈爾濱 150001； 2.哈爾濱工程大學 計算機科學與技術學院，黑龍江 哈爾濱 150001)

針對某些堆芯關鍵參數(如功率分布)的不確定性分析需要采用抽樣統計的方法，不可避免地需要進行大量重復計算，使得分析過程面臨計算量巨大、耗時很長的問題，需要得到精確不確定性分析結果的同時減少計算量。作為人工智能的核心，機器學習在近幾年來得到了迅速的發展，眾多成果已用到反應堆工程計算領域，其中以機器學習預測模型替代數值求解過程[1-3]及經驗模型預測[4-5]的應用最廣泛。其優勢在于通過分析輸入數據與輸出參數的關系訓練機器學習模型，訓練后的模型在給定輸入參數后可快速得到目標輸出參數。當實際計算過程所需的計算代價較大時，采用訓練完備的機器學習預測模型替代傳統計算過程，可以獲得精度和效率的雙重收益。采用高保真技術的反應堆精細化中子輸運計算由于能群、計算網格劃分更加精細，需要較長的計算時間，基于抽樣統計方法使用此計算工具開展針對堆芯關鍵參數的不確定性分析必然會耗費不可接受的計算時間。因此采用機器學習預測模型替代精細化中子輸運計算以實現進一步不確定性分析是解決上述問題的一種可行方案。

本文開展了通過訓練機器學習模型替代精細化中子輸運計算的方法研究及分析，以實現對堆芯關鍵參數的準確、快速預測。構建2D、3D C5G7基準算例預測模型，實現以宏觀核截面為輸入的堆芯有效增殖因子預測；以C5G7單組件瞬態問題宏觀核截面及瞬態參數為輸入的全局功率、全局反應性、中子代時間和緩發中子份額預測。同時，選擇平均絕對誤差、平均標準誤差及方差偏離矩陣為評價指標檢驗預測模型的預測效果。從計算精度、計算時間2個角度對預測結果進行評價，同時從訓練樣本數目及機器學習模型2個角度分析預測效果的影響因素。

1 集成學習模型及評價指標

本文建立的預測模型為集成學習模型，即通過構建、結合多個學習器來完成學習任務，以達到預測能力和魯棒性明顯優于單一學習器的學習性能。

1.1 集成學習模型-堆疊法

堆疊法[6]是一種分層模型集成學習框架。以2層堆疊為例，第1層由多個基學習器組成，它的輸入為原始訓練集，每個基學習器進行K折訓練；第2層模型以第1層基學習器的輸出作為訓練集進行再訓練，最終得到完整的堆疊法模型。圖1為單個基學習器采用5折交叉驗證的堆疊法模型示意圖，簡要過程為：

圖1 堆疊法示意Fig.1 Sketch map of stacking framework

1) 使用一個基礎模型作為Model進行5折交叉驗證，5折交叉驗證的含義為取4折作為訓練數據，另外1折作為測試數據。交叉驗證過程為：基于訓練數據模型，隨后基于訓練數據訓練生成的模型對測試數據進行預測。在第1次交叉驗證完成后會得到關于當前測試數據的預測值，記為a1；

2) 對數據集自身的測試數據進行預測得到預測值，這部分預測值將作為第2層模型測試數據的一部分，記為b1。5折交叉驗證即為將以上過程實施5次，最終生成針對測試數據a1、a2、a3、a4、a5，對測試數據的預測結果為數據b1、b2、b3、b4、b5。

3) 在完成對Model1的整個步驟之后，得到的a1、a2、a3、a4、a5，即為對整個訓練數據的預測值，將它們的組合記為矩陣A1；將b1、b2、b3、b4、b5相加取平均值，得到矩陣B1。則A1即為第2層模型的訓練數據，B1為第2層模型的測試數據。

以上為單個基學習器的stacking模型實施流程，通常情況下stacking模型的第1層包含多個基學習器，即可能存在Model2至ModelN，對于這些模型可重復以上步驟，在整個流程結束之后將得到A2,A3，…,AN,B2,B3，…,BN矩陣，隨后將A2,A3，…,AN合并為N列的矩陣作為訓練數據，B2,B3，…,BN合并為N列的矩陣作為測試數據，并讓第2層模型基于新的訓練數據及測試數據進行訓練和測試。

本文采用的集成學習中的堆疊模型分為2層：第1層采用的基學習器為彈性網絡模型和梯度推進模型，每個基學習器采用5折交叉驗證：彈性網絡是一種同時使用L1范數和L2范數作為先驗正則項訓練的線性回歸模型，能實現輸入參數數量遠小于輸出參數數量時的精確預測；梯度推進模型屬于集成學習中的boosting算法，在大樣本、復雜非線性系統中具有較大的優勢[7]；第2層采用的元學習器為最小絕對值選擇與收縮算子模型，它是一種有偏估計方法，由于在模型中加入了懲罰項，使變量集合簡化、稀疏，從而使模型更加精煉[8]。本文構建的堆疊融合模型，可以集中各基學習器的優勢，在實現精確預測的同時保持計算的穩定性。

本文采用開源Python發行版本Anaconda 4.9.2中的科學包Jupyter 6.2.0進行以上算法建模及運行。

1.2 評價指標

機器學習模型的預測結果可能受到多種因素影響而導致學習率不高、預測結果不準確等問題，通常情況下，隨著訓練數據的增加，模型預測結果的準確率會接近于1，最終在一定范圍內將不再變化[9]。為了更好地對學習效果進行評價，應當選擇恰當的評價指標，本文使用的評價指標為均方誤差MSE、平均絕對值誤差MAE及回歸平方和與總平方和之比Q2分別為：

(1)

變量含義同均方誤差，MAE值越小表明訓練效果越好。

方差偏離矩陣Q2為：

(2)

對于通量分布、堆芯棒功率等分布量，以誤差分布作為評價指標，即計算堆芯每個棒位置處的均方根誤差，最后給出均方根誤差的分布。均方根誤差為：

(3)

式中r為棒在堆芯中的位置。

2 C5G7問題預測驗證與分析

本文選擇的驗證問題為由經濟合作與發展組織核能機構發布的確定論中子輸運計算驗證基準題C5G7問題。采用建立的訓練模型針對上述計算問題開展輸入參數和輸出響應的特征捕捉，可以學習得到代替精細化物理計算過程的預測模型。本節將就C5G7的3個計算問題，對預測模型的計算精確性、計算效率進行驗證，并開展訓練效果影響因素分析。

2.1 C5G7問題

2.1.1 C5G7穩態問題

C5G7穩態問題[10]包括2D、3D問題。基準問題是由16個燃料組件構成的小型水堆全堆芯問題，因此它可以簡化為由4個燃料組件構成的1/4堆芯。2D模型的規模為64.26 cm×64.26 cm，組件寬度為21.42 cm，下、右邊界采用真空邊界條件，上、左邊界采用全反射邊界條件以實現堆芯的對稱結構。單個燃料組件采用17×17方形柵元布置，其中包括264個燃料柵元，24個導向管柵元(用于插拔控制棒)和1個位于中心的測量通道。對于UO2燃料組件，僅有一種燃料類型，而MOX燃料組件則包括4.3%、7.0%和8.7% 3種富集度的燃料。2D C5G7模型的幾何布置如圖2所示。

圖2 C5G7柵元排布及2D堆芯布置Fig.2 Fuel pin configuration and 2D layout of C5G7

3D C5G7模型則是在2D模型的基礎上將燃料區延伸至軸向高度為128.52 cm，同時堆芯的上部和下部增加厚度均為21.42 cm軸向反射材料，最終構成規模為64.26 cm×64.26 cm×171.36 cm的三維計算問題。3D C5G7模型的幾何布置如圖3所示。

圖3 3D C5G7布置Fig.3 3D configuration of C5G7 problem

本文對2D/3D C5G7問題進行訓練所采用的輸入數據為堆芯材料宏觀截面，輸出數據為全堆芯有效增殖因子。

2.1.2 3D C5G7單組件瞬態問題

C5G7單組件瞬態問題描述了單個UO2組件在0.1 s內，由控制棒插入深度為20 cm變為堆芯活性區不插棒狀態的過程。堆芯四周采用全反射邊界條件，堆芯上下由于布置了反射材料因此采用真空邊界條件，3D C5G7單組件瞬態問題的幾何描述如圖4所示。為了實現動態過程的精確模擬，需要對事件發生時間進行時間步的精確劃分，本文選擇的時間步長為5 ms。

圖4 3D C5G7單組件幾何布置Fig.4 Geometry layout of 3D C5G7 single assembly

不同于穩態問題，本文對3D C5G7單組件瞬態問題進行訓練所采用的輸入數據除了堆芯材料宏觀截面外，還包括中子速度、緩發中子份額、緩發中子能譜及緩發中子衰變常數；輸出數據包括堆芯相對功率、反應性、有效緩發中子份額、中子代時間。

2.2 訓練規模及參數設置

穩態及瞬態問題的輸入參數和輸出參數均不一致，本文對2種問題進行訓練和驗證時采用的訓練樣本和測試樣本數目也不相同。表1展示了不同計算問題的輸入、輸出參數規模及訓練、測試數目設置情況。對于訓練樣本和測試樣本的選擇目前尚無定論，一般按照8∶2或7∶3的比例進行設置，本文根據7∶3并結合實際產生的總樣本數設置了不同問題的訓練參數。

表1 訓練規模及參數設置Table 1 Lists of parameters used in training and testing for cases

對于2D及3D C5G7穩態問題，其輸入-輸出關系為：

1) 材料宏觀截面(輸入)：包括7個能群結構的UO2、富集度分別為4.3%、7.0%、8.7%的MOX燃料、導向管、測量通道、慢化劑及控制棒等棒束的材料總截面Σt、吸收截面Σa、裂變截面Σf、中子產出截面υΣf、群-群的散射截面Σg→g′、裂變譜χ，共672個參數(值為0也包含在內)；

2) HNET-穩態模塊(計算模型)：自主開發的高保真精細化穩態中子輸運計算程序[11]；

3) 有效增殖因子keff(輸出)：有限大增殖介質的中子增殖過程中，某代中子數與相鄰前一代中子數之比，用于描述堆芯鏈式反應的狀態。

對于3D C5G7單組件瞬態問題，其輸入-輸出關系為：

1) 材料宏觀截面及瞬態參數(輸入)：包括7個能群結構的UO2、富集度分別為4.3%、7.0%、8.7%的MOX燃料、導向管、測量通道、慢化劑及控制棒等棒束的材料總截面Σt、吸收截面Σa、裂變截面Σf、中子產出截面υΣf、群-群的散射截面Σg→g′、裂變譜χ、中子速度v；8個能群結構的緩發中子份額、緩發中子能譜及緩發中子衰變常數，共1 248個參數(值為0也包含在內)；

2) HNET-瞬態模塊(計算模型)：自主開發的高保真精細化瞬態中子輸運計算程序[12]；

3) 輸出數據為各時間步下堆芯相對功率、反應性、有效緩發中子份額、中子代時間，共128個參數。

2.3 計算精度驗證

2.3.1 2D C5G7穩態問題

2D C5G7穩態問題的輸入參數為7個宏觀截面，輸出參數為有效增殖因子。訓練集樣本數目為108個，測試集樣本數目為60個。圖5直觀地展示了由預測模型得到的預測值與精細化物理計算結果對比，可以看到，2種計算方法得到的結果基本一致，這在一定程度上表明了建立的機器學習模型實現了對2D C5G7穩態問題輸入、輸出關系的特征捕捉，具有實現替代精細化計算的可能性。

圖5 2D C5G7堆芯有效增殖因子預測驗證Fig.5 Prediction test for keff of 2D C5G7 problem

通過引入評價指標可以定量判斷模型的訓練效果。針對2D C5G7穩態問題堆芯有效增殖因子預測結果的MAE、MSE及Q23種評價指標計算結果分別為1.537 60×10-5及1.104 93×10-9，Q2值為0.999 59。文獻[3]中開展的反應堆物理計算替代模型研究針對65個少群宏觀截面，采用200組訓練樣本及100組驗證樣本對組件無限增殖因子進行預測，評價指標MAE及Q2分別為4.11×10-4及0.993。

與文獻[3]預測效果對比可以發現，本文建立的訓練模型針對2D C5G7穩態問題實現了多輸入參數、少訓練樣本下對計算結果的精確預測。

2.3.2 3D C5G7穩態問題

3D C5G7穩態問題的輸入參數同樣為7群宏觀截面，輸出參數為有效增殖因子。訓練集樣本數目及測試集樣本數目與2D問題一直，分別為108和60。圖6展示了3D問題下由預測模型得到的有效增殖因子預測值與精細化物理計算結果的對比，2種計算方法得到的結果同樣基本一致。

圖6 3D C5G7堆芯有效增殖因子預測驗證Fig.6 Prediction test for keff of 3D C5G7 problem

3D C5G7穩態問題有效增殖因子預測結果的MAE、MSE及Q23種評價指標計算結果分別為1.546 68×10-5、1.103 19×10-9、0.999 56。結果表明本文所建立的訓練模型針對3D C5G7穩態問題同樣能實現堆芯有效增殖因子的精確預測。

2.3.3 3D C5G7單組件瞬態問題

對于3D C5G7單組件瞬態問題，利用預測模型實現了堆芯在不同時間下的相對功率、反應性、有效緩發中子份額及中子代時間的預測。圖7展示了預測結果與高保真瞬態計算程序的相對誤差，可以看到，堆芯功率的相對誤差小于1.2%，反應性的相對誤差小于0.014%，有效緩發中子份額相對誤差小于0.065%，中子代時間相對誤差小于0.59%，誤差均處于可接受范圍內。這也表明了本文所建立的訓練模型同樣實現了3D C5G7瞬態問題堆芯關鍵參數的精確預測。

圖7 瞬態參數預測相對誤差Fig.7 Relative error of transient parameters prediction results

2.4 計算效率驗證

基于機器學習的替代模型代替傳統高保真計算需要保證精度和計算時間的雙重收益，本節對計算時間進行測試。表4列出了不同計算問題使用傳統計算方式與預測模型的計算時間對比，可以看到對于2D、3D C5G7問題訓練模型進行60組結果預測所需的計算時間在70 s以內，遠低于采用高保真計算使用的時間；對于瞬態問題而言，由于輸出參數增多，因此時間收益相較于穩態問題有所減少，但相比于傳統計算仍有33倍的計算時間收益。這表明采用基于機器學習的預測模型代替高保真計算所帶來的計算時間收益十分可觀，也表明了預測模型達到了精度與計算時間的雙重收益，滿足了開展進一步不確定性分析的2項基本前提。

表2 計算時間對比Table 2 Comparison of computational time

值得一提的是，訓練時間應當包含訓練樣本產生的時間，即進行多次高保真計算的時間，這也帶來了一定的計算代價。然而，以訓練模型開展后續分析相比采用高保真計算仍能減少大量計算時間。綜合考慮，采用訓練模型替代高保真計算的方案仍是可以采納的。

2.5 訓練集數目對訓練效果的影響

訓練樣本數是模型訓練過程中非常重要的參數，它對模型的訓練效果具有很大影響。通常情況下，訓練樣本數越大，模型訓練效果越好，但無限制的增加訓練樣本數勢必會增加很多不必要的訓練時間；同時，采用數量不足的訓練樣本數則達不到好的訓練效果，因此有必要進行訓練樣本數目對訓練效果影響的分析。

表3列出了在不同訓練樣本數目下，對2D C5G7問題keff的預測模型評價指標的計算值，測試樣本數均為60。隨著樣本數目增加，MAE、MSE值逐漸減小，Q2值逐漸增大，結果符合預期。從圖8可以更直觀的看到這樣的變化趨勢，但在訓練樣本數目為20時，Q2值大于采用30、40個訓練樣本。事實上，這并不能說明樣本數目為20時即可得到合理的預測結果。雖然Q2在訓練樣本數目為20時較大，但預測結果與實際計算結果的方差相對誤差為11%，這對后續的不確定性分析結果會帶來負面的影響。同時，從圖8中可以發現當訓練樣本數大于100時，Q2值已經很接近于1，此時預測結果和實際計算結果的方差相對誤差為1%，即采用大于等于100個訓練樣本時，模型的訓練效果是符合要求的。因此，在評價訓練效果時，應綜合考慮預測精度與方差分布，避免由訓練模型引入新的不確定度。

表3 訓練集數目對訓練效果的影響Table 3 Effect of training data number on prediction accuracy

圖8 訓練樣本數目對預測效果的影響Fig.8 Prediction accuracy change with number of training sets

2.6 不同模型對訓練效果的影響

本文采用2層堆疊的集成學習模型開展針對堆芯關鍵參數的預測，第1層計算模型分別為ENet、GBoost及KRR，第2層計算模型為lasso。事實上，不同模型可以單獨完成訓練、預測工作，因此本節將分別單獨使用以上4種模型對2D C5G7問題的堆芯keff進行訓練和測試，訓練樣本數均為100，測試樣本數均為60。表4展示了采用以上模型單獨預測的評價指標計算值。可以發現，KRR及GBoost方法建立的預測模型預測效果較差，但基于ENet、lasso方法的預測模型可以得到很好的預測結果。

表4 訓練模型對訓練效果的影響Table 4 Effect of training model on prediction accuracy

單就結果而言，采用單一方法建立的預測模型也可以得到精確的預測結果，但堆疊模型相比單一模型具有更顯著的優勢，如不需要太多的調參和特征選擇，模型魯棒性更強等。因此，堆疊模型仍是推薦的訓練模型。

3 結論

1) 本文開展的樣本數目對預測效果的分析結果表明，當訓練樣本數目為100時可以得到滿足要求的預測結果，堆疊模型在實現精確預測的同時能保證其魯棒性。

2) 以基于機器學習的預測模型替代抽樣統計過程中大量重復的精細化物理計算開展不確定性分析是可行的研究方案。

本文建立的堆疊模型為進一步實現以核截面為不確定性來源的堆芯關鍵參數不確定性量化與分析提供了良好的工具基礎。該模型的應用范圍并不限于此，進一步可開展共振計算、燃耗計算、物理熱工耦合計算的替代計算。同時，在開展規模更大、耗時更長的計算結果預測時，本文建立的替代模型將獲得更大的時間收益。