韋子豪,王 端,王東東,潘翠杰
(1.中國原子能科學研究院,北京 102413;2.核工業研究生部,北京 102413)
壓水堆核電廠投入運行后,在每個堆芯壽期末均要停堆換料。堆芯換料過程中,將產生大量的換料方案,而換料方案的選擇直接關系到核電廠運行的安全性和經濟性。一個優選的堆芯裝載方案,可較好地展平堆芯功率分布、增加卸料燃耗深度、延長堆芯換料周期、提高燃料利用率,從而提高核電經濟性和安全性。
堆芯優化裝載方案的搜索是一項十分費時的工作,因為嚴格說來,壓水堆堆芯換料方案設計是一個多變量、非線性的動態整數規劃問題,現代數學理論已證明這類問題是NP-難問題[1-2]。典型3區裝料的壓水堆堆芯含193個燃料組件,即使在堆芯1/4對稱布置且無可燃毒物的條件下,也有1043量級的堆芯換料方案,而且當堆芯燃料通道增加時,換料方案將會以指數甚至更多的形式增加[3],這使得在眾多的換料方案中快速找到1個全局最優的方案變得十分困難。
目前,國際上已圍繞該問題利用各種優化方法開展了大量的研究,總的來說,主要分為確定性優化算法和隨機優化算法,如線性規劃、非線性規劃、動態規劃、直接搜索、專家系統等屬于確定性優化算法[4-6],而隨機優化算法則主要包括模擬退火算法、遺傳算法、粒子群算法、神經網絡等[7-14]。但迄今為止,尚未有快速搜索出全局最優解的通用方法,其中有兩個需要解決的主要問題。一是堆芯換料方案的優化建模問題,如何根據實際問題的需要建立合適的優化模型。優化模型主要包括目標函數、約束條件和決策變量三要素。在換料優化模型中,決策變量就是換料方案,另外兩個要素的選擇需仔細研究。二是優化算法的選擇和計算效率問題。在換料優化模型中,目標函數與部分約束條件不能用堆芯排布方式的明確數學函數表示,它們的值僅能是已知堆芯排布方式的前提下,通過求解復雜的多群中子擴散方程和燃耗方程獲得,這不僅導致依靠數學表達式的確定性優化算法使用受限,同時也是整個優化計算中最耗時的部分,使得快速、準確搜索出全局最優解成為一個非常困難的問題。
本文主要研究優化算法的選擇和效率問題,自主設計一整套堆芯換料的智能優化程序,首先利用神經網絡算法建立預測模型,再設計遺傳算法快速搜索最優的堆芯燃料組件排布方式。通過神經網絡與遺傳算法結合的復合搜索算法改變傳統依賴物理過程計算堆芯參數的做法所帶來低效和復雜化的不利局面,快速高效地提供數量和質量兼備的優化換料方案,為智能化篩選換料方案提供新的可行方案。
為維持反應堆可控鏈式核反應,堆芯裝載多組燃料組件,這些組件按照一定的規律以陣列形式排布。如圖1所示,壓水堆堆芯上百個燃料組件的排列形狀通常近似八邊形陣列。

圖1 堆芯燃料組件陣列Fig.1 Core fuel assembly matrix
燃料組件排布存在著一定的規律,以秦山二期壓水堆為例,其排布方式呈1/8對稱分布,因此僅需考慮其中1/8部分,如圖2所示。共涉及4大類燃料組件E、F、G、H,除了新燃料組件H,同一燃耗次數的組件燃耗深度也存在差別,因此增加數字編號區分。由于組件僅能在固定區域進行調換,故將這1/8部分分為3個區域,以Ⅰ、Ⅱ、Ⅲ表示(標記為D12的組件是對稱中心,固定不變,不將其納入考慮范圍之內)。

圖2 1/8堆芯組件排布示意圖Fig.2 Schematic diagram of 1/8 core fuel assembly arrangement
圖2中共有16種不同燃耗深度的組件,為簡便起見,若假設5個H組件的燃耗深度也不相同,則根據排列組合關系,可計算排布方式總數為6!×10!×4!=6.27×1010,其中,6!、10!、4!分別是Ⅰ、Ⅱ、Ⅲ區組件的排布方式總數。對于627億種組合的排布方式,如果全部計算選擇最優解,則會失去時效性與經濟性。按如下3個步驟進行改進:1) 在所有排布方式中均勻抽樣,選擇10 000種組合作為樣本;2) 建立高精度BP神經網絡模型,根據樣本數據進行訓練,模擬堆芯排布-關鍵參數對應關系;3) 使用遺傳算法快速搜索最優換料方案。
樣本抽樣過程為:從3個區所有排布方式中各隨機抽取1種方式并進行組合,獲得1個樣本,重復抽樣10 000次,即可獲得10 000個一定程度上反映整個組合空間的樣本。從這個過程可看到,作為樣本數據僅提供輸入輸出的對應關系,不會涉及其中包含的復雜物理過程。
利用CASMO5程序計算每種排布方式下的堆芯關鍵參數,包括有效增殖因數(keff)、組件功率峰因子(Rad)和棒功率峰因子(FΔH),這3種參數與堆芯排布方式一起,生成學習樣本。
1) 樣本總體特征
從樣本數據讀取其中獨立的1/8排布信息作為輸入,keff、Rad和FΔH作為標簽,考察樣本標簽的分布規律。圖3示出原始數據分布直方圖。
2) 樣本平衡處理
對于keff數據進行細致考察發現,其最大值為1.202 02,最小值為1.123 65,精度為0.000 01,意味著樣本最多7 838個不同的數據,而實際是僅有4 304個,并且根據圖3可知,值集中,一些值僅有1、2個。對樣本進行如下平衡處理:(1) 該值在樣本中僅有1個,進入訓練集;(2) 該值在樣本中存在N個,進入訓練集的概率為(1/N)-1/2。如果未進入訓練集,則進入測試集。
最終獲得的keff數據訓練集和測試集樣本分布如圖4所示。可看出,訓練樣本分布更均勻,方差更大。

圖3 原始數據分布直方圖Fig.3 Raw data distribution histogram

圖4 處理后的樣本分布直方圖Fig.4 Processed data distribution histogram
同樣對Rad、FΔH數據進行平衡處理,但由于它們更加集中,為不破壞原始數據的分布,采用簡化的平衡樣本的方式:從標簽在均值±1左右的樣本中選取1 000組作為測試集,其余作為訓練集。
3) 輸入向量選擇
在文獻[15]中,將E、F、G、H 4種組件表示-1到1的等差數列,即-1、-1/3、1/3、1。這種方法過于簡化,對預測精度有較大影響,因此本文使用4位二進制數表示16種不同燃耗深度的組件,對應情況列于表1。這樣,每個包含20個組件的排布方式x=(1,1,1,1,0,0,0,0,…,1,0,1,1,1,0,0,1)T均是1個80維的輸入向量。

表1 組件與編碼對應關系Table 1 Correspondence between component and coding
4) 輸出向量的處理
每個樣本標簽為3個關鍵參數,為提高精度,采用3套不同網絡訓練,分別設置輸出層為1個節點,輸出向量長度為1。
理論研究證明,多層前饋網絡即能以任意精度逼近任意復雜度的連續函數(萬能逼近定理)[16]。但要使神經網絡能較精準地擬合1個復雜函數,設置其網絡結構如隱藏層數、節點數量等,并進行有效訓練是比較困難的,目前尚未有通用辦法,僅能根據經驗反復嘗試。為更快速獲得較好的模型,對訓練的模型提出評價指標。本文主要使用的誤差是絕對誤差:
(1)
式中,ypredict、ylabel分別為預測值、真實值。
1) 敏感性與適應性分析
利用神經網絡算法解決問題時,要求輸入輸出數據之間強的相關性,一般需進行特征分析和選擇,本問題中輸入數據為燃料組件的類型與排布規律,對輸出數據有決定性的影響,因此是強相關的特征,比較可靠。輸出數據的3個堆芯關鍵參數在精度上不同,對噪音干擾的敏感性也不同。keff精度要求高,因此模型預測產生的誤差也同樣要求嚴苛,而對于其余兩個參數數據則可降低標準。本文認為,模型對keff預測的平均誤差不高于0.05,對FΔH預測的平均誤差不高于0.5,而對Rad預測的平均誤差不高于0.5即是可靠的。
另外,對于與訓練樣本相似或相鄰的數據,擬合結果是相對精確的,但對于整個解空間的極值(這些值在訓練樣本中沒有或很少),擬合的結果會差一些,而且受模型優劣的影響較大。
2) 評價指標
用以下6個指標評價模型的優劣。(1) 訓練誤差(MAEtrain)與測試誤差(MAEtest),即10 000組樣本在訓練模型中最終對應的訓練誤差和測試誤差,它們能反映訓練情況,之間的差距還能反映擬合情況,即過擬合或欠擬合。(2) 模型最大誤差(MAEmax),即模型在所有樣本中的最大絕對誤差。(3) 模型平均誤差(MAEavg),即模型在所有樣本中產生的平均誤差。(4) 驗證數據(td),即一些未參與訓練的關鍵數據,如CASMO5計算出的當前堆芯最優數據為:keff=1.127 21,Rad=1.20,FΔH=1.388。另外一些表征樣本邊界的數據,即樣本中標簽極值,如keff樣本最大值為1.202 02,Rad樣本最小值為1.28,即為樣本邊界數據(bd)。(5) 超閾值誤差數量(ot),即超過一定絕對誤差閾值(如0.1)的樣本數量。(6) 模型預測極限(lp),即模型能預測到的整個解空間中最值。由于實際全局最優解難以確認,以遺傳算法搜索到的最值作為替代。
泛化誤差為模型在一些未知數據上預測的誤差,是模型優劣最核心的評價標準。但需預測的數據僅有輸入已知,標簽未知,無法計算誤差。因此,1個模型優劣的評判也是神經網絡算法所面臨的現實問題之一。除了通常使用的訓練誤差和測試誤差之外,實驗發現以上幾個指標也可從不同方位給出一個模型的訓練情況和泛化能力。
因此,在使用遺傳算法搜索最優結構時,使用的適應度應綜合考慮這些指標,為此,可擬定適應度函數如下:
(2)
式中:MAEtrain為訓練誤差;MAEtest為測試誤差;MAEmax為樣本總體的最大誤差;MAEavg為樣本總體的平均誤差;|xp1-bd|為對邊界值的誤差;|xp2-td|為對驗證數據的誤差;ot為誤差超過0.05的樣本數量;lp為模型預測的極限值;δi(i=1,2,…,8)為各項系數,代表各項的影響權重。
對于keff數據,由于精度較高,基于奧卡姆剃刀原則,單隱層BP神經網絡即可訓練出一個比較滿意的模型。但Rad和FΔH的樣本精度低、方差小,為更好學習到數據特征,使用深度神經網絡(DNN)模型,其拓撲結構如圖5所示。

圖5 包含3個隱藏層的BP神經網絡結構示意圖Fig.5 Schematic diagram of BP neural network structure with three hidden layers
經實驗確定的keff、Rad和FΔH預測模型的結構列于表2。
圖6為keff、Rad和FΔH訓練過程誤差和學習率的變化,使用dropout減少過擬合,使用持續減小的學習率使模型盡可能向全局最優收斂,可見訓練的收斂情況表現比較理想,得到的最佳模型是較為可靠的,其中模型的評價指標列于表3。模型對樣本的平均誤差均小于0.05,與敏感性分析中預期的值相符合,對邊界值和特殊的驗證數據的預測誤差也很小,說明模型擁有一定的泛化能力。

表2 3個預測模型網絡結構信息Table 2 Structure information of three predictive models
注:隱藏層后括號內為dropout節點重置率,即不進行參數調整的神經元比例

圖6 keff(a)、Rad(b)與FΔH(c)模型訓練中收斂情況數據圖Fig.6 Figures of convergence in model training in keff (a), Rad (b) and FΔH (c) data
表3 模型的評價指標
Table 3 Evaluation index of models

關鍵參數MAEtrainMAEtestMAEmaxMAEavgbdtdotlpkeff0.001 10.001 10.003 20.001 11.201 02(1.202 02)1.128 13(1.127 21)01.206 89Rad0.0440.0430.110.0441.3(1.28)1.23(1.20)301.143FΔH0.0490.0460.120.0481.403(1.411)1.422(1.388)141.321 7
注:括號內為對應樣本的標簽值
1) 優化目標[17-18]
本文數據針對的是秦山二期平衡循環(換料周期12月),堆內無可燃毒物,因此以參數keff、Rad和FΔH作為優化目標。具體目標是,在可行域范圍內,搜索使堆芯參數最優的組件排布方式,使keff最大化,Rad和FΔH最小化。以f(x)表示目標函數,則表達式為:
f(x)=?kefffkeff(x)+?RadfRad(x)+
?FΔHfFΔH(x)+δ
(3)
式中:δ為1個常量,用以保證f(x)為正值;?keff、?Rad、?FΔH為3個堆芯參數的適應度所占的權重。以f(x)的最大化為目標,其約束條件為?keff≥0,?Rad≤0,?FΔH≤0。
權重的正負決定對應參數的優化方向,權重絕對值的大小決定參數的重要程度。
2) 適應度選擇與尺度變換
單純以目標函數為適應度會發生早熟和停滯現象,因此引入簡單的適應度尺度變換。最終確定適應度計算公式為:
fitness=f4(x)=(2fkeff(x)-
3fRad(x)-5fFΔH(x)+24)4
(4)
遺傳算法認為生物進化的本質是染色體的選擇、交配與突變。借助該思想,遺傳算法主要包含的算子為選擇、交叉和變異算子。而在第1.1節已指出,將不同燃耗次數組件以E、F、G和H表示,不同燃耗深度的組件則使用字母與數字結合標記,因此1個20個組件的排布方式可使用按照一定順序連接形成1個基因序列,即對應遺傳算法的編碼,如“H02F34F26F35-G02E34H01H08G10G04F32H05E07G19F07-H04F01E36F29G27”。
1) 選擇算子
使用輪盤賭選擇與最優保存策略結合的選擇算子。輪盤賭選擇使用隨機采樣方法,即蒙特卡羅方法,從N個個體中有放回的隨機抽取N次,作為下一算子的輸入,每個個體被抽取的概率由其適應度決定,其公式為:
(5)
式中:fi為第i個個體經過尺度變換后的適應度;Pi為其被選中的概率。
進行輪盤賭選擇前,將當前適應度最高的個體保存,作為新一代種群的1個個體,保證適應度最高的個體總參與遺傳、變異,這就是最優保留策略。
2) 交叉算子
使用保留基因片段的方法進行交叉,即在子代中保留父代被選擇保留的基因段。在父代(1)中選擇保留的基因段,在子代(1’)中將保留并保持位置不變,并將父代(2)中與這些保留基因相同的基因去掉,然后將剩下的基因按照一定的順序(如原來的順序)填充到子代(1’)余下的位置形成新的子代(1’),同理可獲得子代(2’)。與一般的交叉算子不同,為使生成的子代依然滿足約束條件,即包含在解空間中,要在保留1個父代的基因后,對另一個父代基因進行篩選。
每次保存的基因均應保持一定的長度以保持父代基因的特征,不同的長度意味著信息的不同的保留率,可在實驗中檢驗不同信息保留率對搜索結果的影響。本文對3個區的基因設置不同的交叉概率參數,以向量Pc=(pc1,pc2,pc3)T表示,其中元素分別表示Ⅰ、Ⅱ、Ⅲ區的交叉概率。在迭代進行過程中,根據一定規則改變交叉概率,有利于算法平衡全局搜索能力與搜索時間的關系,本文在前期以較高的概率進行交叉,而后期則逐步降低交叉的發生概率。當發生交叉時,使用保留上一代1/2的基因信息。
3) 變異算子
設置變異概率向量Pm=(pm1,pm2,pm3)T,分別對染色體的3個區間進行變異操作,即隨機選擇兩個不同位置,交換該位置上的1對基因順序從而生成新的1代。變異前后個別組件發生調換而產生變異,與交叉類似,變異隨著時間的進行,其概率有所變動以改變其局部搜索的能力,在迭代后期可適當增加這一能力以盡可能在局部上達到最優。
遺傳算法的流程如圖7所示。
對于圖7,有如下兩點說明。1) 通過抽樣數據獲取學習數據集,通過數據處理生成訓練集和測試集,使用訓練集對預測模型進行訓練,當在測試集的絕對平均誤差小于1個預設的預定值時獲得最終的預測模型。2) 遺傳算法通過隨機數初始化一定規模的種群,獲得輸入Xm,將Xm作為預測模型的輸入獲得預測值,從而根據適應度計算公式計算出適應度。種群根據適應度進行3個算子的變換,并記錄過程中最優的解。當迭代次數或最優解滿足停止條件時,停止迭代,輸出結果。
設置原始雜交概率Pc=(0.6,0.9,0.5)T,變異概率Pm=(0.05,0.1,0.05)T,種群數量N=100~400,迭代次數為200~1 000,進行多次搜索,搜索時間為5~10 min。其中1個結果的情況如圖8所示。表4列出優選的換料方案對比。圖9為其中1個最優解(預測值1)對應的排布方式,圖9a為基因序列按照從右到左、從上到下的順序填充的1/8堆芯,圖9b為基因序列對應的完整堆芯。
圖8a指出遺傳算法在迭代300次左右、耗時數min內可搜索到最優解,圖8b指出適應度限制了對keff的最大化搜索,因此如果需要可改變適應度函數以獲得更高的keff。表4列出本文提供的集中優選堆芯方案和對應的參數的預測值,在誤差范圍內與目標值相當或更好。

圖7 神經網絡-遺傳復合算法的總流程圖Fig.7 General flowchart of neural network-genetic composite algorithm

圖8 復合算法的適應度變化曲線Fig.8 Fitness curve of composite algorithm
表4 優選的換料方案對比
Table 4 Comparison of preferred refueling scheme

項目keffRadFΔH目標真實值1.127 211.201.388目標預測值1.128 19±0.001 101.23±0.041.422±0.048最優預測值11.136 76±0.001 101.166±0.0401.353±0.048最優預測值21.134 89±0.001 101.187±0.0401.360±0.048最優預測值31.135 03±0.001 101.184±0.0401.335±0.048

圖9 其中1個最優解(預測值1)對應的排布方式Fig.9 One of optimal solutions (prediction 1) corresponding arrangement
實驗表明,采用神經網絡-遺傳復合算法可預測不同堆芯排布方式下的關鍵參數,并利用生物的遺傳進化機制,在大量排布方式下快速搜索更優的堆芯排布方式。1個訓練很好的神經網絡模型對超出樣本數據所在區間外的預測精度會下降,導致搜索的范圍不能太大。解決這一問題的方法是,根據已獲得的最優結果再次訓練網絡、再次搜索,不斷迭代,向更優的結果靠近,直到搜索到全局最優解,這是下一步的工作。
建立自適應神經網絡提供堆芯參數,再結合改進的遺傳算法進行堆芯換料方案的篩選,快速、準確選擇較優的堆芯換料方案。其結果可推廣至任意堆芯堆型、任意關鍵參數預測,并且可很容易地優化修改目標函數,篩選符合實際需要的換料方案。該方法不需特定領域的專業知識,僅要有足夠多的數據,即可實現精度較高的快速預測,但需對模型作進一步的敏感性分析和適應性研究。