陳法國,李國棟,楊明明,韓 毅,梁潤成
(中國輻射防護研究院,山西 太原 030006)
屏蔽設計本質上是一個帶有約束條件的多目標尋優過程。一方面,屏蔽效果與射線類型和能譜分布密切相關;另一方面,在實際工程中,屏蔽設計還需要考慮重量、體積、成本、結構性能等多種設計限定條件。傳統的屏蔽設計通常以設計者的個人經驗為基礎,優化設計效率相對較低,設計結果的全局最優性受人因等不確定因素影響較大。為此,遺傳算法(GA)、人工神經網絡算法(ANN)、特征統計算法(CSA)等多種優化計算方法[1-5],被應用于屏蔽結構設計、復合屏蔽材料組分優化等問題求解中,以實現屏蔽設計方案的自動優化。
中子屏蔽設計需要考慮中子的慢化、吸收以及次級光子的屏蔽,需要同時采用快中子非彈性散射占優的重核材料、中能中子彈性散射截面大的輕核材料以及熱中子俘獲截面大的中子吸收材料,是一個典型的優化求解問題[6]。以中子復合屏蔽材料組分優化設計為例,以屏蔽劑量和材料密度最小化為目標,開展了遺傳算法在多目標屏蔽優化設計中的應用研究。
利用優化算法自動設計中子復合屏蔽材料的組分比例,使其以最小的密度實現對特定能譜中子的最好屏蔽效果。中子源項設為單向面源,能譜分布為Watt裂變譜;中子平均能量為2.348 MeV,能譜分布如圖1所示;復合材料由聚乙烯、鉛、硼、鋰、鐵、鋁6種常見材料組成。假設各組分均勻混合且混合前后總體積不變,用單位注量的中子穿過厚度為30 cm平板結構后的中子和次級光子劑量之和(定義為屏蔽劑量D)來表征屏蔽效果,屏蔽結構示意圖如圖2所示。

圖1 入射中子的能譜分布Fig.1 Energy spectrum of incident neutron

圖2 屏蔽結構示意圖Fig.2 Schematic diagram of shielding structure
中子復合屏蔽材料組分優化的數學函數模型可表示為:
minF(X)=[D(X),ρ(X)]T
s.t. 0≤xi≤1,∑xi=1,1/ρ(X)=
∑xi/ρi,i=1,2, …6
(1)
式中,F(X)為兩維設計目標,D(X)為組分為X屏蔽材料對應的屏蔽劑量,X=[x1,x2, …x6]分別為6種組分的質量分數,ρ為材料密度。
傳統的屏蔽優化問題求解,采用加權法將多目標降維到單目標進行處理。加權法簡化了問題求解,但子目標權重因子的設置具有一定的主觀性,并且優化結果在理論上只是Pareto最優解集中的一個解,在權重因子變化后需要再次運行優化過程進行求解;而多目標優化算法可以在單次運行過程中,求解Pareto最優解集[7]。為此,采用遺傳算法直接對式(1)中的雙目標優化問題進行求解。
遺傳算法是由Holland提出的模擬生物種群遺傳和自然選擇規律而演化的一種隨機搜索方法。為提高優化效率、避免過早收斂,在基本遺傳算法的基礎上,發展了向量評估遺傳算法(VEGA)、多目標遺傳算法(MOGA)、非支配排序遺傳算法(NSGA)、基于距離的Pareto遺傳算法(DPGA)等[7]。其中,NSGA是學者Srinivas和Deb基于非支配排序思想提出的,通過分類排序可提高適應度較高的個體保留至下一代的概率,并采用適應度共享策略使個體分布更均勻,以保持種群的多樣性,避免過早收斂;為優化NSGA的時間復雜度較高、未采用精英策略、優化過程和尋優效率對共享參數設置較為敏感的問題,Deb等人又提出了快速非支配排序遺傳算法(NSGA-II)[8]。
利用NSGA-II建立式(1)數學模型所表示的中子復合屏蔽材料組分優化算法。優化算法主要包括屏蔽計算模塊和遺傳算法模塊兩部分,算法流程圖如圖3所示。遺傳算法的主要構成要素有基本運行參數、染色體編碼方式、個體適應度評價、遺傳算子。

圖3 中子屏蔽材料組分優化算法流程圖Fig.3 Flow chart of optimization algorithm for neutron shielding material component
基本運行參數包括種群大小M(通常設為20~100)、終止進化代數T(通常設為50~200)、交叉因子Pc、變異因子Pm。交叉因子和變異因子決定了父代和子代個體的差異性,直接影響算法的收斂性和尋優性能;增加差異性可擴大搜索空間,但不利于收斂,反之亦然。
編碼和解碼是應用遺傳算法的關鍵,用于建立優化問題的設計變量與遺傳基因的映射關系。屏蔽材料組分優化是一個定和約束問題,即∑xi=1,0≤xi≤1;若直接用質量分數xi進行交叉或變異操作,子代種群中會出現不滿足約束條件的個體。為此,采用球坐標轉換方法,引入一個采用浮點數編碼方式的5維向量θ=[θ1,θ2…θ5]T,0≤θi≤π/2,來表征設計變量X=[x1,x2…x6]T,其中:
xi=(sinθ1·sinθ2·…sinθ5)2,i=1
xi=(cosθi-1·sinθi·…sinθ5)2,i=2,3,4,5
(2)
xi=(cosθ5)2,i=6
在執行交叉或變異算子時,直接對θi進行操作,然后按照式(2)將其解碼為質量分數xi,用于計算復合材料的密度和屏蔽劑量。通過球坐標轉換,既能滿足定和約束條件,又可保證交叉和變異操作的隨機性。
遺傳算法通過解碼得到特定個體的材料組分后,調用屏蔽計算模塊來計算對應的屏蔽劑量。在整個優化算法流程中,對每代種群中的新個體都需要調用一次屏蔽計算模塊;若直接采用蒙特卡羅(MC)程序作為屏蔽計算模塊,多次調用會增加時長。為了兼顧屏蔽計算精度與運行時間,屏蔽計算模塊采用了基于BP(誤差逆向傳播算法)神經網絡算法的屏蔽快速計算方法。
針對屏蔽優化問題的具體模型,建立可根據材料組分比例預測屏蔽劑量的BP神經網絡,以便在優化過程中快速調用。首先,在整個搜索空間內對6種材料的組分比例進行均勻抽樣,并根據抽樣參數建立MC模型計算屏蔽后劑量D,形成初始訓練樣本。其次,利用訓練樣本對BP神經網絡進行訓練和學習,尋找組分比例與劑量之間的計算規律。最后,利用訓練后的網絡來預測任意組分比例所對應的屏蔽計算結果。
利用297個訓練樣本訓練后,BP神經網絡預測結果和MC計算結果的相對偏差在-5.2%~8.1%范圍內;隨機抽取33種新組分比例作為測試樣本,神經網絡預測結果與MC計算結果的相對偏差在-3.9%~5.1%范圍內。在執行優化算法時,將個體編碼進行解碼后,可直接調用屏蔽快速計算模塊,并將計算結果用于個體適應度評價。
適應度函數是用來衡量種群中的個體接近最優解的優良程度。NSGA-II算法中,利用非支配序rank和擁擠度兩個屬性來評價個體適應度。
對于兩個變量個體X1和X2,若所有子目標函數滿足fi(X1)≤fi(X2),且存在子目標函數使fj(X1) 在對兩個個體進行選擇時,非支配序小的個體占優,非支配序相同則擁擠距離大的個體占優;非支配序小可保證被選擇的個體較優,擁擠距離大可保證個體多樣性,避免過早收斂于局部最優。 選擇、交叉、變異是遺傳算法中的基本算子。 NSGA-II中,選擇算子采用錦標賽方法。從父代種群M(種群大小為N)中隨機選擇s1個個體,進行適應度評價,再從s1個個體中選擇占優的s2個個體進入子代種群Ms(N>s1>s2);重復上述過程,直到子代種群的個體數量與父代種群一致,然后按照設定的交叉因子和變異因子對子代種群中的個體進行交叉、變異操作。 交叉算子對兩個父代個體中的編碼向量進行部分重組,將其中的父代基因按照特定算法處理,生成兩個新的子代個體;變異算子用一個基于特定概率分布的隨機數替代父代個體中的某個基因位,從而產生新的子代個體。NSGA-II中采用模擬二進制的交叉算子(SBX)和多項式變異算子。為了防止種群中的較優個體因交叉或變異操作而丟失,NSGA-II引入了精英策略。在精英策略中,將父代種群與子代種群混合(種群大小為2N),重新計算混合種群中所有個體的非支配序和擁擠度,并選擇占優的N個個體作為新一代種群。 對每一代種群重復執行以上操作,當滿足迭代條件后,將當前種群中的個體進行解碼,就可輸出屏蔽材料組分比例的優化設計結果。 針對Watt譜分布的入射中子,在未屏蔽的條件下,單位注量對應的周圍劑量當量為388.6 pSv·cm2。6種組分材料的密度取值列于表1,屏蔽計算時各組成元素的同位素設置為天然豐度。 利用Matlab軟件編譯了優化算法程序,并利用CPU主頻為3.7 GHz的計算機運行。種群大小設為100,交叉因子設為20,變異因子設為20,以材料密度最小和屏蔽后劑量最小為目標,經過100代后中子復合屏蔽材料組分的優化方案分布如圖4所示;圖中顯示了初始方案以及優化1代、10代、50代、100代后的方案(a包含了所有數據,b為局部放大圖)。由圖可知,初始方案在搜索區域內隨機分布;隨著迭代次數的增加,優化方案逐漸接近Pareto前沿區域,并且在該區域內較為均勻地分布。 表1 組分材料的密度(g/cm3)Tab.1 Densities of component materials 圖4 中子屏蔽材料組分的優化結果Fig.4 Optimization results of neutron shielding material component 對于每一代種群,用該代種群中所有個體對應的密度和屏蔽后劑量的乘積之和,來表示種群個體與Pareto前沿的接近程度,即種群優良指數。如圖5所示,隨著進化代數的不斷增加,種群優良指數逐漸降低;當進化代數超過40代時,優良指數最終趨于穩定。 圖5 種群優良指數隨進化代數的變化Fig.5 Figure of merit for population against evolutional generation number 在最終優化方案的100個個體中,劑量最小方案、密度最小方案、典型中間最優方案的組分比列于表2。由表可知,密度最小方案中,復合材料的主要組分幾乎全部為鋰;劑量最小方案中,復合材料的主要組分是86.5%的聚乙烯、10.4%的天然硼以及少量的鋁和鉛;對于一系列中間方案,隨著復合材料密度的增大,鋰的含量逐漸減少,聚乙烯的含量逐漸增加,同時含有少量的硼。 Watt譜分布中子能量主要在2 MeV以下,屏蔽材料的有效原子量越小,慢化作用越好。6種材料中鋰的密度最低、有效原子量最小且6Li與熱中子的相互作用截面大,單位質量的屏蔽效果最好;對于密度最小方案的理想情況,優化結果幾乎全部為鋰,符合理論預期。劑量最小方案更符合實際情況,并且實際應用經驗中含硼聚乙烯是對該能譜中子屏蔽效果最好的材料,與優化結果一致。上述定性結果說明了該屏蔽材料組分優化算法的可行性。 在得出優化方案后,根據優化的材料組分進行MC建模計算,并將屏蔽計算結果與優化結果進行對比。在最終的100個可行解中,遺傳算法優化結果與MC計算結果的相對偏差在-14.9%~4.8%范圍內。 1)引入基于BP神經網絡的快速屏蔽計算模塊,避免了優化過程中對MC計算程序的大量調用,節省了計算時間。計算時間包括297個初始樣本的MC計算時間、神經網絡訓練時間和遺傳算法程序運行時間,分別為127.1 min、10.6 s和71.6 s。若直接調用MC程序,100個個體的種群進化100代,需要運行10 100次MC程序,則整個優化過程至少耗時72.4 h。但神經網絡訓練樣本不足會導致預測結果與MC計算值的相對偏差增大,甚至誤導遺傳算法的尋優方向而陷入局部最優。最終的優化解集和訓練樣本相比,屏蔽后劑量的BP神經網絡預測結果與MC計算結果的相對偏差從-5.2%~8.1%變到-14.9%~4.8%,整體呈現變大趨勢。BP神經網絡算法的初始訓練樣本在整個搜索空間內隨機抽樣得到,在整個變量區域內的分布是相對均勻的;隨著種群的不斷進化,優化解集逐漸向Pareto前沿區域靠近,初始訓練樣本在該區域的代表性不足,基于初始訓練樣本的神經網絡預測結果的不確定度就會增大。為了降低神經網絡預測的不確定度對優化結果的影響,可采用分階段優化或部分調用MC計算程序的處理辦法。 2)在適應度評價中引入不確定度因子,可提高種群的多樣性,避免因屏蔽計算誤差導致的過早收斂。對于屏蔽計算,不論是MC方法、點核積分方法,還是基于BP神經網絡算法的快速預測,都存在不可忽視的計算誤差。對個體適應度進行評價時,若直接以屏蔽后劑量的大小為判據,可能造成對個體優劣的誤判,從而誤導尋優過程;在算法調試時,多次獨立執行優化程序,就出現了劑量最小方案中的組分比例隨機陷入不同局部最優解的情況。為此,在適應度評價中引入了一個劑量不確定因子δD,若兩個個體X1和X2對應的劑量子目標存在如下關系: |D(X1)-D(X2)|/(D(X1)+D(X2))≤δD 則認為X1和X2劑量子目標的適應度相等;根據屏蔽計算誤差設置不確定因子δD,在本工作中為計算誤差的二分之一。通過引入不確定因子,在選擇算子中對個體基于錦標賽進行選擇時,可避免因屏蔽計算誤差造成的個體優劣誤判,從而保證種群的多樣性。 3)NSGA-II算法在求解雙目標和三目標問題時,具有很好的優化效率。但實際應用中,屏蔽設計通常還需要考慮成本、體積、材料性能等多種因素,優化子目標可能超過三個。NSGA-II算法在較多目標求解中,存在以下問題:①隨著子目標數量增加,一個隨機種群中非支配解的數量占比將成指數增長,經過幾代進化后會很快占據整個種群,使搜索過程變慢;②隨著設計目標維度增加,計算時間成本將大大增加;③對于三維以下的設計目標,可以借助于優化解集分布前沿的圖形化表示,輔助對優化結果的評價和決策;但設計目標維度超過三個后,不容易用圖形直觀表示。對于三目標以上的優化問題,可以考慮參考自適應加權法,將多維設計目標降維后進行求解;也可以采用Deb等人最新提出的基于參考點的非支配遺傳算法[9]。 利用NSGA-II遺傳算法,開展了屏蔽劑量和材料密度最小化的雙目標中子復合材料組分優化研究,以Watt裂變中子譜以及聚乙烯、鉛、硼、鋰、鐵、鋁6種常見材料組成30 cm厚平板結構為例,進行了自動優化設計;優化結果定性說明了基于算法的屏蔽材料組分優化方法的可行性。本工作目前只關注遺傳算法在多目標屏蔽設計中的應用,暫未關注屏蔽結構長期使用后的中子活化對屏蔽效果的影響及深穿透等問題,這些問題可在后續工作中繼續關注。以設計者的經驗為先驗條件,結合該優化算法,可更好地進行輻射屏蔽設計,為更高效確定優化設計方案提供現實指導。2.5 遺傳算子
3 優化結果及對比分析
3.1 優化結果



3.2 結果對比分析
4 討論和結語
4.1 討論
4.2 結語