基于聚類和隨機搜索優化的核反應堆數字孿生參數反演模型

龍家雨,宋美琪,*,柴 翔,劉曉晶,妥艷潔,3

(1.上海交通大學 智慧能源創新學院,上海 200240;2.上海交通大學 機械與動力工程學院,上海 200240;3.國家電力投資集團有限公司,北京 100029)

數字孿生(digital twin)是以數字化方式創建物理實體的虛擬模型,充分利用物理模型、傳感器、運行歷史等數據,集成多學科、多尺度對物理實體在現實環境中的行為進行模擬的仿真過程。作為虛擬空間中對實體產品的鏡像,其反映了相對應物理實體產品的全生命周期過程[1-2]。核電數字孿生有利于核電站以較低成本實現高可靠性、可用性與可維護性[3],可在核電站的預測性運行與維護、自主運行與控制等領域發揮作用,如實現小型模塊化反應堆自動、自主和實時性能要求[4]或結合先進PSA(概率風險評估)方法構建實時決策支持系統[5]等。

核電站為實時監測運行狀態而布置大量的傳感器,產生的數據種類多、總量大,造成存儲空間需求大、傳輸效率不高、數據分析復雜度高等問題,基于核電站構建的數字孿生系統同樣面臨此類問題,如何優化傳感器的布置也需要進一步研究。核電領域常采用數據降維與反向求解的特征工程進行數據處理與分析[6-7]。Peng等[8]利用稀疏自動編碼器(sparse auto encoder)對仿真程序PCTRAN模擬的核電站的幾種不同瞬態參數提取特征,實現對核電站異常狀態的檢測。Yang等[9]用一維卷積神經網絡(CNN)重建核電站傳感器信號,通過逐步減少傳感器進行實驗,在16種傳感器數據確定了5種影響最大的參數。李翔宇等[10]通過改進特征工程算法建立了一種核電站瞬態運行參數數據壓縮和數據復原方法,結合主成分分析法與高斯回歸過程方法對秦山300 MWe全范圍仿真機產生的瞬態運行數據進行壓縮與高精度復原。

目前此類研究在特征提取上面臨如下的問題:通過逐次實驗的方式獲取特征參數,依賴長期數據,運算量大、效率較低;采用主成分分析等方式進行數據壓縮,在壓縮過程中改變了原參數的值,雖然提升了傳輸效率且能達到很高精度,但也因此失去了特征參數的可解釋性。

針對上述問題,本文基于多物理場耦合程序實現空間熱離子反應堆虛擬模型的構建并獲取數字孿生數據,結合K-means聚類算法與ANN(人工神經網絡),提出一種能夠在保留原數據的情況下高效提取特征參數并對其他非特征參數進行反演的方法,對空間熱離子反應堆堆芯的4個區域的溫度參數分別建立參數反演模型并取得較高精度的反演結果。

1 數據集建立

本文的參數反演模型是對空間熱離子反應堆的數字孿生堆芯溫度數據建立的。對堆芯物理系統的數字孿生計算通過上海交通大學搭建的多物理場耦合數值模擬平臺[11]完成。

1.1 空間熱離子堆堆芯系統

與常規壓水堆相比,空間熱離子反應堆體形較小,堆芯三維幾何模型如圖1所示?？臻g熱離子反應堆由37根熱離子燃料元件和ZrH慢化劑芯塊組成。通過旋轉控制鼓引入反應性,從而啟動空間熱離子反應堆。燃料熱量產生與傳遞過程如下:燃料芯塊發生核裂變反應,產生熱量,熱量通過發射極、接收極、氣隙傳遞給冷卻劑,冷卻劑通過對流換熱將會帶走部分熱量,剩余熱量傳遞給慢化劑芯塊,慢化劑通過熱輻射向太空散失部分熱量。

圖1 空間熱離子反應堆堆芯的幾何模型[11]Fig.1 Geometric model of space thermionic reactor core[11]

空間熱離子反應堆堆芯中最核心的部件是熱離子燃料元件,其結構示于圖2。熱離子燃料元件為單節全長多層套筒結構的熱離子燃料元件,采用高濃縮UO2作為燃料,NaK-78作為冷卻劑。發射極材料為單晶鉬和化學氣相沉積制成的鎢表面層,接收極由多晶鉬制成,接收極外包有Al2O3絕緣體薄層。

圖2 單節熱離子燃料元件[11]Fig.2 Single-section thermionic fuel element[11]

1.2 數據獲取

使用基于OpenFOAM開發的多物理場耦合程序對空間熱離子反應堆構建數值模擬平臺。網格劃分、工況設置、相關方程及計算可見文獻[11]?；緟档某绦蛴嬎阒蹬c設計值吻合良好,在設計范圍內。

20%熱管失效工況下的仿真結果中包含燃料、發射極、接收極、冷卻劑4個區域的溫度計算結果,采樣時間間隔為0.1 s,瞬態計算從750 s開始到2 425 s結束,共16 750個時間節點。沿燃料元件軸向(z軸)均勻分布的100個測點,高度從0到0.375 m。剔除最初原始數據發生突變的時間點,最終生成了16 748×1 601維的數據集,對應1 600個參數的16 748個時間節點。

用相同的方式獲取15%熱管失效工況下的堆芯溫度數據,包含從750.2 s到1 893.9 s共11 438個時間節點,參數種類與20%熱管失效工況下的相同,將用于參數反演模型的效果驗證。

2 參數反演模型構建

本文構建的參數反演模型結構如圖3所示。不同于自動編碼器與數據壓縮復原方法,本文構建的模型在特征提取環節通過機器學習方法快速選取出特征參數并且保留了其原始數據。

圖3 本文的參數反演模型Fig.3 Parameter inversion model proposed in this paper

依據如圖4所示流程構建參數反演模型,包括特征提取與神經網絡搭建兩部分,并在最后進行模型效果驗證。

圖4 參數反演模型搭建流程Fig.4 Process of constructing parameter inversion model

2.1 基于K-means的特征選取環節

特征選取是一個數據降維環節,將原始數據轉換為低維度的特征向量,從而提高數據傳輸與分析效率,也節約了存儲空間。K均值聚類算法(K-means算法)是一種無監督的分類方法,在設定聚類數后可自動通過算法完成不同參數類別的劃分。不同類簇中的中心參數表征了該類參數中的典型值,選擇到聚類中心有最小平方歐氏距離的參數作為特征參數。由此,通過K-means聚類可以高效選取特征參數。聚類數的選擇影響著聚類最終效果的好壞,選擇輪廓系數(silhouette coefficient)作為評估聚類效果的指標,并比較不同聚類數的輪廓系數以選擇最佳聚類數。

特征選取環節包括如下步驟:將空間熱離子反應堆的數字孿生數據歸一化、使用K-means算法分類、使用輪廓系數指標評估聚類效果、提取各類最接近聚類中心的參數作為特征參數,流程如圖5所示。

圖5 特征提取流程Fig.5 Process of feature extraction

首先,分別對4個位置的溫度參數進行歸一化處理,將數據映射到[0,1]之間。使用如下公式:

(1)

K-means算法的核心思想是找出K個聚類中心c1,c2,…,cK,使得每個數據點xi和與其最近的聚類中心cv的平方距離和被最小化(該平方距離和被稱為偏差D)[12],主要步驟[13]如下。

1) 聚類中心初始化:隨機指定K個聚類中心c1,c2,…,cK。

2) 對樣本歸類:每個樣本xi僅選擇1個聚類中心S(t)作為歸屬,找到與它具有最小平方歐氏距離的聚類中心cv,并將其分配到所標明的類cv:

?n,1≤n≤K}

(2)

3) 更新聚類中心cn:將每個cn移動到其標明的類的中心:

(3)

4) 計算偏差D:

(4)

5) 判斷偏差D是否收斂:如果D收斂,則返回c1,c2,…,cK并終止本算法;否則,返回步驟2。

輪廓系數是用來評估聚類效果好壞的一種指標,衡量聚類中樣本的內聚度和分離度。內聚度ai為第i個樣本到同簇中其他點的距離的平均值,分離度bi為第i個樣本到其他簇中點的距離的平均值,第i個樣本的輪廓系數[14]為:

(5)

式中:Si越接近1,說明樣本i的聚類合理;Si為0時,說明樣本i處于兩簇的邊界;Si為-1時,說明樣本i應該被分到其他的簇。對所有樣本的輪廓系數取平均值,得到聚類結果的總輪廓系數。

以接收極為例,輪廓系數隨聚類數的變化如圖6所示。在聚類數為31時,輪廓系數達到最大,表明類數選為31類時的聚類效果最好。

圖6 接收極的輪廓系數隨聚類數的變化Fig.6 Variation of silhouette coefficient of receiving pole with number of clusters

確定聚類數后,選取每類中到聚類中心平均距離最小的參數曲線作為特征參數,保存特征索引。參數索引與幾何位置一一對應,其中軸向方向序號z-1～z-100分別代表0.003 75 m的整數倍,如z-3代表軸向0.003 75 m×3處位置;徑向方向序號r-1代表徑向中心處,r-2、r-3、r-4分別代表3個不同徑向位點。仍以接收極為例,提取的31個特征參數的對應位點列于表1。

從表1可看出,特征參數沿軸向分布較為均勻,在徑向分布上,位于中心與4號徑向位點的特征參數占所有特征參數的74.2%,表明這兩個徑向位點的參數值得關注。由于上述特征參數均可在堆芯中確定其所屬位置,因此也可確認每個特征參數的含義。

對各區域分別進行聚類,4個區域的聚類數最終確定如表2所列。重復上述流程,得到4個區域的特征參數及其對應的位點信息。將特征參數與非特征參數重新編排為新的數據集,用于下一步的參數反演。

表2 不同區域的聚類數Table 2 Clusteringnumbers of different regions

2.2 搭建神經網絡反演模型

完成對數據的降維后,需要進行反向求解,通過使用降維后的數據求出對應的物理場分布與相應探測器的理論測量值,從而完成數據的反演。本文為4個區域分別搭建全連接神經網絡進行參數反演,并使用隨機配置法優化神經網絡結構和學習率,根據各區域的反演結果評估參數反演模型。

神經網絡的計算單元完成輸入xi與權重wi的加權求和,添加偏置bi并經過激活函數Φ得到輸出yi,構成下一層的輸入或者最終的輸出[15],如下式:

(6)

全連接代表每個神經元與相鄰層的所有神經元相連,且同層中的神經元互不連接。選擇MSE(均方差)作為損失函數,以歐氏距離表達反演的預測值yi與真實值y′i的距離:

(7)

使用SGD(隨機梯度下降)函數,通過反向傳播更新神經網絡權值,每個連接的權值通過自身減去學習率與損失函數對該權值的一階偏導數進行更新。使用原始數據輸入神經網絡,溫度值的范圍在[750 ℃,2 500 ℃],由此選擇神經網絡各隱藏層的激活函數為Leaky-ReLU。將16 748個時間節點數據按12 000∶2 000∶2 748的比例劃分為訓練集、驗證集、測試集,在前14 000個時間節點中按照12 000∶2 000的比例隨機劃分為訓練集與驗證集。收斂條件設置為兩輪MSE之差小于1×10-6且MSE<100。以接收極為例,初步搭建神經網絡。經嘗試,當學習率大于1×10-7時訓練易發散。采用1×10-7的學習率、網絡結構為3層、每層各20個節點的全連接神經網絡,設置1 500輪訓練。損失隨迭代輪數的變化如圖7所示,最終訓練集與驗證集損失分別為0.716 4和0.844 8。然后圍繞上述參數,對學習率、神經網絡層數與節點數進行優化。

圖7 訓練損失與驗證損失隨迭代次數的變化Fig.7 Variation of training loss and validation loss with number of iterations

對神經網絡結構的優化方法包括試驗法、窮舉法、啟發式算法、修剪與構造算法等[16],其中能夠自動搜索來完成的有隨機配置法、系統化搜索、窮舉法、啟發式方法。本文使用隨機配置法選取神經網絡的隱藏層層數及各層節點數、學習率。4個區域各自的神經網絡在初步訓練后實現收斂,圍繞初步訓練的神經網絡參數進行配置,為4個區域的隱藏層數、各層節點數、學習率設定搜索區間進行隨機配置。

根據Robert Hecht Nielson的證明,1個3層的BP網絡可以完成任意M維到m維的映射[17]?？紤]到堆芯溫度數據的輸入維度遠小于輸出維度,將隱藏層數的選擇范圍設定為2～5。

隱藏層節點數的選取具有高度任意性,節點數的范圍參照如下經驗公式[18]并適當擴大范圍:

(8)

式中:n1為隱藏層節點數;n為該層的輸入神經元數;m為該層的輸出神經元數;a為1～10之間的常數。節點搜索范圍確定為[5,25]或[5,30]。

根據初步訓練的配置,將學習率的搜索范圍設為[1×10-9,1×10-7]。各區域的神經網絡參數搜索范圍列于表3。最大輪數設置為20 000,滿足收斂條件時結束。

表3 不同區域的神經網絡參數配置隨機搜索區間Table 3 Random search interval of neural network parameter configuration for different regions

神經網絡訓練與優化流程如圖8所示。

圖8 神經網絡訓練與優化流程Fig.8 Neural network training and optimization process

對各區域分別設置30種隨機組合,綜合考慮訓練集收斂速度最快、收斂時MSE最小的組合作為神經網絡配置。當收斂輪數與最終訓練損失不能同時達到最小時,選擇相近的訓練損失中收斂輪數最小的作為配置。以發射極為例,從30種配置結果中選出損失最小的5組,如表4所列。從中可以看出,組別為4號的配置具有最小的最終訓練損失與最快的收斂速度,故作為最后的神經網絡配置,訓練損失變化如圖9所示。

表4 發射極訓練損失最小的幾組配置Table 4 Several configurations of minimal training loss in emitter

圖9 發射極神經網絡最終配置下的訓練損失變化Fig.9 Training loss variation of emitter neural network with final configuration

通過上述流程獲得4個區域的不同神經網絡配置,結果列于表5。

3 反演結果

對反演值與實際數值進行比較,采用如下方式:對角線圖能直觀地比較反演值與實際值;相對誤差頻數分布直方圖展示不同誤差區間的數據數目;相對誤差的RMSE(均方根誤差)反映了反演的總體誤差,計算方法為:

(9)

將上述訓練好的神經網絡用于15%熱管失效工況下的堆芯溫度數據,各部分的反演結果如下。

3.1 燃料溫度反演結果

燃料溫度范圍在[1 370 ℃, 2 450 ℃]。燃料溫度的反演值與真實值貼近直線y=x(圖10),相對誤差分布如圖11所示。最大相對誤差在±1.4%內,94.34%的數據反演值相對誤差在1%以內,相對誤差均方根不超過0.6%。

圖10 燃料溫度的反演值與真實值對角線圖Fig.10 Diagonals of inverse and actual values of fuel temperature

圖11 燃料溫度相對誤差頻數分布直方圖Fig.11 Relative error frequency distribution histogram of fuel temperature

3.2 發射極溫度反演結果

發射極溫度范圍在[1 350 ℃, 2 220 ℃]。發射極溫度的反演值與真實值的對角線圖如圖12所示,相對誤差分布如圖13所示。最大相對誤差不超過1.2%,相對誤差均方根在0.4%左右。整體誤差較燃料更小。

圖12 發射極溫度的反演值與真實值對角線圖Fig.12 Diagonals of inverse and actual values of emitter temperature

圖13 發射極溫度相對誤差頻數分布直方圖Fig.13 Relative error frequency distribution histogram of emitter temperature

3.3 接收極溫度反演結果

接收極溫度范圍在[750 ℃, 935 ℃]。接收極溫度的反演值與真實值的對角線圖如圖14所示,誤差分布如圖15所示。99.75%的反演值相對誤差不超過0.6%,相對誤差均方根不超過0.2%。

圖14 接收極溫度的反演值與真實值對角線圖Fig.14 Diagonals of inverse and actual values of receiving pole temperature

圖15 接收極溫度相對誤差頻數分布直方圖Fig.15 Relative error frequency distribution histogram of receiving pole temperature

3.4 冷卻劑溫度反演結果

冷卻劑溫度范圍在[730 ℃, 920 ℃]。冷卻劑溫度的反演值與真實值的對角線圖如圖16所示,相對誤差分布如圖17所示。92.09%的反演值相對誤差在±0.3%以內,98.67%的反演值相對誤差在[-0.3%,0.9%],相對誤差均方根不超過0.2%。

圖16 冷卻劑溫度的反演值與真實值對角線圖Fig.16 Diagonals of inverse and actual values of coolant temperature

4 結論

本文在當前相關研究無法同時兼顧特征參數的自動、高效選取與保留特征參數的可解釋性的情況下,基于多物理場耦合程序構建空間熱離子反應堆的數字孿生,結合聚類算法與隨機搜索優化的人工神經網絡對20%熱管失效工況下的堆芯溫度建立了溫度參數反演模型,并通過15%熱管失效工況數據進行反演效果驗證。結果表明,運用本文流程搭建的參數反演模型,對于空間熱離子反應堆堆芯各區域的溫度參數,反演最大相對誤差在1%上下,相對誤差均方根在0.001 8～0.005 5之間;能實現參數比例小于1∶10的反演,對空間熱離子反應堆堆芯溫度場的反演取得良好的精度,提高傳輸、儲存與分析效率;且獲取的特征參數可通過索引定位到不同構件/區域的軸向與徑向位置,從而得到特征參數的物理意義。本研究的模型可以拓展應用到空間熱離子反應堆其他區域以及其他堆型,高效獲取特定區域的特征參數并保留特征參數的位置信息,為未來小型堆探測器的安裝提供依據,也為研究大規模的核電數據提供了一種可行思路。