李曉琪，鄭友琦，杜夏楠，王永平

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

與壓水堆相比，快堆的運行經(jīng)驗相對匱乏。與此同時快堆的設(shè)計目標(biāo)除壽期更長、功率峰因子更小以外還包括獲得更大的輻照量、更深的燃料增殖或焚燒等。在工程實踐中，還需要考慮諸多約束條件，如熱工、瞬態(tài)安全準(zhǔn)則等?？於褍?yōu)化設(shè)計歸根結(jié)底是一個多變量、多目標(biāo)、有約束的尋優(yōu)問題?？紤]到龐大的搜索空間，枚舉法并不可行；考慮到問題的多峰特性，梯度下降法等數(shù)學(xué)模型亦不可取，因此需要引入智能優(yōu)化算法[1-8]。在本文中使用的遺傳算法作為經(jīng)典的智能優(yōu)化算法，適合于解決這類多輸入、多目標(biāo)、帶約束的問題。

反應(yīng)堆燃料管理包括兩種類型的問題：1) 關(guān)注組件的制造參數(shù)以及燃料裝量等的“Out-of-core”問題；2) 對現(xiàn)成的組件關(guān)注布置方案的“In-core”問題。對這兩種問題進(jìn)行優(yōu)化的方法分為兩類，分別是確定性方法和隨機(jī)性方法。早年關(guān)于壓水堆的燃料管理優(yōu)化方法中，確定性方法發(fā)展相對較早，如對燃料和毒物進(jìn)行脫耦的哈林脫耦方法[9]，通過在中子擴(kuò)散方程中引入0～1變量，將換料方案優(yōu)化轉(zhuǎn)化為數(shù)學(xué)規(guī)劃問題的直接優(yōu)化方法[1,10]等。隨機(jī)性方法多為啟發(fā)式，包括模擬微觀粒子熱運動特性實現(xiàn)優(yōu)化的模擬退火算法、對搜索方向進(jìn)行標(biāo)記進(jìn)而提高搜索效率的禁忌搜索算法[11]、模擬動物群落通過標(biāo)記信息尋找目標(biāo)的粒子群算法[12]，以及模擬生物種群進(jìn)化規(guī)律的遺傳算法等。上述隨機(jī)性算法通過從某種自然或社會現(xiàn)象中獲得啟發(fā)以實現(xiàn)對方案的定向搜索，具有較強的通用性，基本不依賴于物理背景本身，適用于本文要解決的快堆優(yōu)化設(shè)計問題。近年來蓬勃發(fā)展的人工智能推動了神經(jīng)網(wǎng)絡(luò)在堆芯換料優(yōu)化及設(shè)計優(yōu)化領(lǐng)域中的應(yīng)用，通過樣本訓(xùn)練代理模型進(jìn)行堆芯參數(shù)預(yù)測以實現(xiàn)優(yōu)化成為堆芯優(yōu)化研究的一重要方向[13-14]。

遺傳算法最早由從事自適應(yīng)系統(tǒng)研究的Professor Holland團(tuán)隊提出。1967年Holland教授的學(xué)生Bagley首次在論文中應(yīng)用了遺傳算法，并使用了包括選擇、交叉、變異、占優(yōu)、倒置在內(nèi)的遺傳算子。1992年，Poon與Parks最早將遺傳算法用于反應(yīng)堆燃料管理優(yōu)化問題。在早年用遺傳算法進(jìn)行反應(yīng)堆優(yōu)化設(shè)計的研究中，優(yōu)化目標(biāo)基本為keff和功率峰因子[15]。

使用遺傳算法進(jìn)行換料優(yōu)化與堆芯設(shè)計的研究中，通常從隨機(jī)或人為給定的一組初始方案出發(fā)，通過不斷迭代改變搜索方向進(jìn)行尋優(yōu)。如對每個組件進(jìn)行獨立編號，以編號的排列作為方案的基因[2]，并配合使用改進(jìn)的算子獲得新基因，從而避免無效方案的產(chǎn)生[16]。早期的研究大多對裝載方案的十進(jìn)制排列進(jìn)行二進(jìn)制編碼，以模擬遺傳現(xiàn)象的本質(zhì)[17]。對于規(guī)模較大并涉及到換料的問題，文獻(xiàn)[18]提出一種對裝載方案基因分段，分別表示組件類別與循環(huán)次數(shù)的方式。由于搜索空間對于整個尋優(yōu)過程是未知的，初始方案的偏移可能導(dǎo)致搜素過程陷入局部最優(yōu)，因此需要慎重確定初始方案。本文提供一種新的編碼方案，對搜索空間中含有的所有個體進(jìn)行編碼，在優(yōu)化的起始在多個自變量維度上均勻抽取初始方案，從而一定程度避免陷入局部最優(yōu)。

1 理論與優(yōu)化框架開發(fā)

1.1 DAKOTA/SARAX耦合

在本文的優(yōu)化框架中，使用了開源工具箱DAKOTA中的多目標(biāo)遺傳算法模塊MOGA。DAKOTA工具箱由桑迪亞國家實驗室開發(fā)，包含優(yōu)化、敏感性分析等功能。優(yōu)化模塊MOGA可作為黑箱使用，由中子學(xué)計算程序SARAX[19-20]對個體(方案)進(jìn)行分析(計算)并反饋給DAKOTA作為遺傳操作的憑據(jù)。SARAX是由西安交通大學(xué)開發(fā)的快堆中子學(xué)計算程序，適用于快堆的臨界和次臨界計算。使用Python腳本將兩個外部程序進(jìn)行耦合，功能是將DAKOTA提供的個體(方案)翻譯為堆芯計算的參數(shù)或建模作為輸入卡片傳遞給SARAX，并調(diào)用SARAX進(jìn)行穩(wěn)態(tài)計算、燃耗計算、反應(yīng)性計算等，在計算結(jié)束后將計算結(jié)果處理成遺傳算法可以讀取的“適應(yīng)度”即衡量方案優(yōu)劣的指標(biāo)返還給DAKOTA進(jìn)行下一步操作。

SARAX是一套中子學(xué)計算程序，包括用于截面產(chǎn)生的模塊Tulip、堆芯穩(wěn)態(tài)計算模塊Lavender及瞬態(tài)模塊Daisy。本文的工作集中在使用預(yù)先準(zhǔn)備的截面實現(xiàn)堆芯層面的參數(shù)或裝料的優(yōu)化方案搜索，使用堆芯穩(wěn)態(tài)計算模塊Lavender。

在燃料管理與堆芯設(shè)計中使用遺傳算法，流程如圖1所示，分為以下步驟：1) 遺傳算法模塊產(chǎn)生第1代隨機(jī)個體(方案)的基因文件，翻譯成堆芯程序輸入后傳遞給堆芯計算模塊；2) 堆芯計算的結(jié)果處理成適應(yīng)度函數(shù)返回給遺傳算法模塊；3) 遺傳算法模塊進(jìn)行選擇、交叉、變異操作，形成包含新個體(方案)的種群；4) 重復(fù)步驟3和4直到種群收斂。

圖1 堆芯設(shè)計和燃料管理中的遺傳算法流程Fig.1 Genetic algorithm flowchart in design of reactor core and fuel management

該流程在本框架中的數(shù)據(jù)流動如圖2所示。

圖2 優(yōu)化框架的流程圖Fig.2 Flowchart of optimal frame

1.2 編碼策略

遺傳算法模仿生物進(jìn)化的過程，其中包含一個重要的步驟即基因與設(shè)計參數(shù)或裝載方案之間的編碼/解碼過程。將二進(jìn)制或十進(jìn)制的編碼翻譯成連續(xù)的堆芯參數(shù)相對簡單，而翻譯成裝載方案則需要一定的策略。

本文提出了一種直接實現(xiàn)十進(jìn)制基因文件與裝載排列或組合方案之間的一一映射的編碼/解碼思路，并用于此框架，具體實現(xiàn)如下。

1) 從十進(jìn)制到全排列

在這種裝載策略下，每一個組件具有獨一無二且連續(xù)的編號，組件數(shù)量與裝載位數(shù)量一致，如圖3所示。使用一種映射算法“康托展開”，實現(xiàn)十進(jìn)制整數(shù)與全排列之間的一一映射。表1列出對3個數(shù)字1/2/3進(jìn)行全排列下排列方式與序號之間的對應(yīng)關(guān)系。此過程可以通過式(1)實現(xiàn)。

圖3 全排列問題Fig.3 Arrangement problem

X=a[n]×(n-1)!+a[n-1]×

(n-2)!+…+a[i]×(i-1)!+

…+a[1]×0!

(1)

式中：a[1]，…，a[n]為全排列時第n位上的組件編號(對組件從1開始編號)；X為該全排列的十進(jìn)制序列號(從0開始)，對于本框架相當(dāng)于基因從十進(jìn)制基因到全排列的過程，需要將表1右側(cè)的十進(jìn)制整數(shù)(即基因)翻譯成左側(cè)的排列方案，進(jìn)而作為組件的布置，此過程可通過康托展開的逆過程實現(xiàn)。

表1 十進(jìn)制整數(shù)與全排列的對應(yīng)Table 1 Correspondence between decimal integer and full permutation

2) 從十進(jìn)制整數(shù)到組合數(shù)

對于大多數(shù)裝載策略，同一類型的組件有多個，都需要參入到排列中。如果直接套用上述的全排列思路，則會產(chǎn)生多個基因?qū)?yīng)本質(zhì)上同一種排列方式的問題，導(dǎo)致大量重復(fù)方案產(chǎn)生，進(jìn)而影響優(yōu)化效果。對于這種組合問題，本文提供一種編碼方案：(1) 先針對一種類型的組件，挑選相應(yīng)數(shù)量的裝載位，此過程是1個組合問題，其中需要進(jìn)行1次十進(jìn)制整數(shù)與組合數(shù)之間的映射，即將十進(jìn)制整數(shù)轉(zhuǎn)換為一系列組合數(shù)，而這個組合數(shù)就是裝載該種類型組件的裝載位，此十進(jìn)制數(shù)是基因的一段；(2) 對另一種類型的組件，在剩下的裝載位上重復(fù)過程1，以此類推；(3) 當(dāng)只剩下一種類型的組件，則直接填充剩下的裝載位。

組合問題的編碼過程如圖4所示。

圖4 組合問題的編碼Fig.4 Coding of combinatorial problem

考慮到如圖5所示的情況，待裝載的組件總數(shù)比裝載位多。解決方案是在原基因前增加一段，用于決定哪些組件不參與此次排列。

圖5 考慮的后備組件Fig.5 Considering of reserved assembly

由于快堆不需要考慮1/4對稱或1/8對稱的情況，上述編碼方案基本可以應(yīng)對快堆中的編碼問題，并且具有以下優(yōu)勢：(1) 有效方案與指定范圍內(nèi)連續(xù)變化的基因一一對應(yīng)，搜索空間中每個個體都有1個獨一無二的編號，利用這些編號可在優(yōu)化初始進(jìn)行均勻抽樣，無需提供初始方案，如圖6所示；(2)基因型的微小變化對應(yīng)表現(xiàn)型的微小變化，該特點對于隨機(jī)性優(yōu)化算法是重要的評估標(biāo)準(zhǔn)。

a——傳統(tǒng)遺傳算法初值選??；b——本文改進(jìn)方法初值選取圖6 初始方案的分布Fig.6 Distribution of initial scheme

1.3 多目標(biāo)處理

遺傳算法廣泛用于多目標(biāo)優(yōu)化問題，包括反應(yīng)堆設(shè)計和燃料管理。由于優(yōu)化目標(biāo)不止1個，如何讓算法綜合考慮每個優(yōu)化目標(biāo)以評估方案的優(yōu)劣，直接影響到優(yōu)化效果。處理方式分為兩種類型：將多目標(biāo)與約束條件按權(quán)重組合成1個適應(yīng)度函數(shù)作為評估準(zhǔn)則，或直接使用多目標(biāo)。本文搭建的框架中，使用了DAKOTA中內(nèi)置的多目標(biāo)適應(yīng)度處理功能domination_count和layer_rank。這兩種功能對被評估的個體在每個優(yōu)化目標(biāo)維度上進(jìn)行比較，將“占優(yōu)”數(shù)量作為評估個體優(yōu)劣程度的指標(biāo)。這與傳統(tǒng)的將多目標(biāo)用1組權(quán)重線性相加合成1個適應(yīng)度函數(shù)相比，不需要人為給定權(quán)重，使用更簡單且不易因為權(quán)重選取不合適導(dǎo)致對多目標(biāo)評估的偏差，是遺傳算法適應(yīng)度函數(shù)處理當(dāng)下發(fā)展的一個主流做法。

1.4 約束條件處理

除優(yōu)化目標(biāo)，方案優(yōu)化問題中可能需要設(shè)置約束條件，對不滿足約束條件的方案直接淘汰。由于遺傳算法通過適應(yīng)度函數(shù)對方案進(jìn)行評估，因此約束條件也將通過作用于方案的適應(yīng)度函數(shù)來實現(xiàn)，稱為罰函數(shù)，即對不滿足某個約束條件的方案在其適應(yīng)度函數(shù)后加上1個數(shù)值，降低該方案在后續(xù)的選擇操作中被選中的概率。對于多目標(biāo)問題，當(dāng)采用將多目標(biāo)與約束條件按權(quán)重組合成1個適應(yīng)度函數(shù)的權(quán)重法處理方式時，將罰函數(shù)按式(2)所示添加在適應(yīng)度函數(shù)中即可。本文中直接使用多目標(biāo)優(yōu)化，則將罰函數(shù)依次添加在不符合某一條件的方案的每一個優(yōu)化目標(biāo)函數(shù)值內(nèi)。

fit=ω1f1+ω2f2+…+P

(2)

式中：ωi為優(yōu)化目標(biāo)i的權(quán)重；fi為優(yōu)化目標(biāo)i在當(dāng)前方案下的值；P為罰函數(shù)，當(dāng)前方案不滿足某個約束條件時P不為0，都滿足則取為0。

2 優(yōu)化框架驗證

2.1 旅行商問題

在對反應(yīng)堆優(yōu)化設(shè)計問題進(jìn)行框架驗證之前，先通過1個旅行商問題(TSP)驗證編碼方案是否合理以及DAKOTA內(nèi)置選項使用是否正確。

此處設(shè)計了1個簡單的TSP：若干城市以1個單位長度為間隔排列成1圈，則最短路徑是依次環(huán)繞這些城市1圈，如圖7所示。該模型的優(yōu)點是容易確定最短路徑與對應(yīng)的距離，并且易于增加或減小問題的規(guī)模。

圖7 簡化旅行商問題Fig.7 Simplified TSP problem

本文考慮了12個城市，即進(jìn)行12個整數(shù)的全排列問題搜索，解空間的大小為479 001 600個方案，最優(yōu)方案數(shù)目為24個，隨機(jī)生成初始方案。適應(yīng)度即為當(dāng)前方案下遍歷12個城市的路徑長度。遺傳算法參數(shù)設(shè)置列于表2。方案搜索過程中，每代最短路徑變化如圖8所示。DAKOTA內(nèi)置小生境選項，即根據(jù)當(dāng)前種群個體的分散或聚集程度，通過調(diào)整個體適應(yīng)度，避免個體分布過于集中，以提高種群的基因多樣性從而避免陷入局部最優(yōu)。測試表明，本文給出的編碼方案可在無初值的條件下，可接受計算時間內(nèi)，配合小生境的使用，對12個城市的TSP收斂到全局最優(yōu)，從而驗證了編碼方案的合理性，以及DAKOTA遺傳算法參數(shù)的選擇的正確性。

表2 TSP的遺傳算法參數(shù)設(shè)置Table 2 Parameter setting in TSP

圖8 城市數(shù)為12的TSP收斂Fig.8 Convergence of 12 cities TSP

2.2 基于ABTR的優(yōu)化與枚舉法

本文基于先進(jìn)燃燒實驗堆ABTR的燃料裝載優(yōu)化，通過與同一問題的枚舉結(jié)果對比，驗證本文提出的編碼方式對于十進(jìn)制與組合數(shù)轉(zhuǎn)換的處理。ABTR是一個用于驗證嬗變技術(shù)的概念堆，堆中有兩種不同富集度的燃料組件，其中高富集度燃料分布在外圈，低富集度燃料分布在內(nèi)圈，周期為4個月，批數(shù)分別為12和15。此外還有一些用于進(jìn)行材料輻照實驗的測試組件分布在活性區(qū)中。堆芯活性區(qū)的1/3示意圖如圖9所示。

圖9 ABTR原始活性區(qū)布置Fig.9 Original active zone layout of ABTR

對給定數(shù)量的3種燃料組件，重新布置其排列方式，目標(biāo)是對keff和功率峰因子擾動最小的前提下在測試組件中獲得最大的中子通量。以1/6堆為對稱條件，解空間有1 260個方案，通過枚舉獲得每個可能方案的3個參數(shù)，作為判斷優(yōu)化搜索是否尋找到最優(yōu)解的依據(jù)。在該問題下，基因文件具有兩段信息，分別用于決定外圈高富集度燃料(未必仍然排在外圈)和內(nèi)圈低富集度燃料在堆芯中的位置。當(dāng)兩種燃料的位置確定后，剩余位置即為測試組件的位置。對3個優(yōu)化目標(biāo)的函數(shù)值采用1∶1∶1的權(quán)重加和為1個適應(yīng)度函數(shù)，其中keff和功率峰因子擾動量的絕對值直接取為優(yōu)化目標(biāo)函數(shù)值，測試組件的中子通量取負(fù)作為優(yōu)化目標(biāo)函數(shù)值，如式(3)所示。

fit=ω1keff+ω2|ppf-ppf0|+ω3(-φtest)

(3)

式中：ω1、ω2、ω3為權(quán)重，ω1=ω2=ω3=1；ppf和ppf0分別為當(dāng)前方案下的功率峰因子和原始方案的功率峰因子；φtest為測試組件處的中子通量。

經(jīng)過對353個個體的計算，獲得了一最終優(yōu)化結(jié)果，堆芯活性區(qū)布置示意圖如圖10所示。此外對整個優(yōu)化空間的1 260個方案依次進(jìn)行中子學(xué)計算，并記錄下keff、功率峰因子、測試組件處的中子通量密度，與最終優(yōu)化結(jié)果一起繪制成三維示意圖，如圖11所示。用kh表示功率峰因子，T_flux表示測試組件處的中子通量密度，紅色點表示搜索出的最優(yōu)方案，在二維坐標(biāo)上的投影如圖12所示。

圖10 優(yōu)化框架提供的最優(yōu)方案Fig.10 The best scheme provided by optimal frame

圖11 最優(yōu)方案在完整解空間中的分布Fig.11 Distribution of optimal schemein whole solution space

圖12 最優(yōu)方案在二維坐標(biāo)的投影Fig.12 Mapping of optimal scheme in two-dimensional coordinate

3 優(yōu)化框架在實際問題中的使用

引入中國實驗快堆(CEFR)作為研究目標(biāo)。CEFR是中國第一座鈉冷實驗快堆，運行周期為80 d，活性區(qū)布置如圖13所示，是一個用于各種燃料、材料輻照實驗與放射性同位素生產(chǎn)的優(yōu)良平臺[21]。快堆的用途之一是同位素焚燒和生產(chǎn)，可以用于焚燒隨著壓水堆的運行產(chǎn)生的大量對環(huán)境有危害的237Np。237Np可通過嬗變轉(zhuǎn)換成能制作同位素電源的238Pu。基于CEFR堆芯，裝載含有237Np的組件以實現(xiàn)238Pu在1個運行周期內(nèi)產(chǎn)量最大化為目標(biāo)，同時保證后備反應(yīng)性、功率形狀分布盡可能不偏離原始堆芯。237Np以NpO2的形式裝載到堆芯中，包括在燃料活性區(qū)和增殖區(qū)的裝載，如圖14所示。通過搜索在哪個組件的活性區(qū)或增殖區(qū)中裝載以及活性區(qū)中的裝載量，實現(xiàn)表3所列的優(yōu)化目標(biāo)與約束條件。其中對功率波動的優(yōu)化與約束是由于考慮Np的裝載影響堆芯的中子價值，進(jìn)而影響控制棒組的價值，而7%的功率波動限值是考慮了CEFR的反應(yīng)性控制系統(tǒng)設(shè)計誤差限給出的經(jīng)驗值。

圖13 CEFR的活性區(qū)布置Fig.13 Layout of active zone in CEFR

表3 優(yōu)化目標(biāo)與約束Table 3 Optimization objective and constraint

圖14 裝載NpO2的組件Fig.14 Subassembly with NpO2

此處采用直接的多目標(biāo)優(yōu)化方法，故不構(gòu)造適應(yīng)度函數(shù)，而是由DAKOTA直接對每個方案的每個優(yōu)化目標(biāo)之間進(jìn)行比較。

使用DAKOTA內(nèi)部提供的rank_method作為選擇算子，個體在keff和功率波動兩個維度上的分布如圖15所示，其中紅框內(nèi)是符合約束條件的個體。顯而易見，符合條件的個體在這種搜索條件下并沒有得到很好的保存，原因是在rank_method內(nèi)，以罰函數(shù)形式出現(xiàn)的約束條件并不能發(fā)揮作用。隨后將選擇算子改為傳統(tǒng)的權(quán)重函數(shù)加和的形式，獲得的個體在keff和功率波動兩個維度上的分布如圖16所示。其中最優(yōu)方案的堆芯布置如圖17所示，對應(yīng)的功率波動分布如圖18所示。在該方案下，237Np的轉(zhuǎn)換率為0.616%，80 d238Pu產(chǎn)量為0.549 kg。

圖15 rank_method個體的keff和功率波動Fig.15 keff and power fluctuation of individual in rank_method

圖16 權(quán)重法個體的keff和功率波動Fig.16 keff and power fluctuation of individual in weighting method

圖17 最終方案的堆芯布置Fig.17 Core layout of final scheme

圖18 最終方案的功率波動Fig.18 Power fluctuation of final scheme

4 總結(jié)

本文搭建了一個用于快堆設(shè)計的優(yōu)化框架，可實現(xiàn)讀入一定范圍的決策變量、優(yōu)化目標(biāo)以及約束條件。搜索過程基于初始均勻抽樣，對應(yīng)提出一套對整個解空間進(jìn)行編碼的編碼方式，可以對排列問題、組合問題進(jìn)行連續(xù)編碼。基于該框架，驗證了簡化的TSP、ABTR裝料問題，從而證明編碼方式的合理性和框架的可行性。最后通過一個基于CEFR實現(xiàn)Np-Pu生產(chǎn)優(yōu)化的算例示范了優(yōu)化框架在工程實際中的具體應(yīng)用。