γ輻照裝置多目標(biāo)優(yōu)化排源算法研究及其應(yīng)用

楊 磊，劉 燕，羅志平，李文革，周毅吉，張美琴

(1.煙臺(tái)大學(xué)，山東 煙臺(tái) 264005；2.湖南省農(nóng)業(yè)科學(xué)院，湖南 長沙 410300)

60Co γ輻照加工技術(shù)在農(nóng)業(yè)、工業(yè)、醫(yī)療衛(wèi)生、食品、環(huán)境保護(hù)等多個(gè)領(lǐng)域應(yīng)用十分廣泛[1],目前總體發(fā)展平穩(wěn)，新建裝置呈現(xiàn)大型化、自動(dòng)化及泛用化等特點(diǎn)[2]。因60Co源半衰期較短，為維持和擴(kuò)大加工能力，需定期補(bǔ)充和更換60Co源棒，而源價(jià)格昂貴，為此提高源的整體利用率十分重要。生產(chǎn)實(shí)踐表明，優(yōu)化60Co輻照裝置排源設(shè)計(jì)是提高源利用率的核心技術(shù)因素，對(duì)優(yōu)化裝置劑量學(xué)性能、改善貨物加工劑量、降低單位成本、提高裝置加工量和經(jīng)濟(jì)效益等均有顯著影響。

排源方法主要包括人工經(jīng)驗(yàn)排源和計(jì)算機(jī)自動(dòng)排源兩大類，行業(yè)公司仍以使用前者為多，此類方法對(duì)工作人員專業(yè)素質(zhì)要求高，人工和時(shí)間花費(fèi)多[3]，研究快速、可靠、結(jié)果方案足夠優(yōu)化的自動(dòng)排源方法成為必然需要。目前國內(nèi)外有基于貪心算法、神經(jīng)網(wǎng)絡(luò)算法(深度學(xué)習(xí)算法)等的排源方法，貪心算法快速有效但全局搜索能力不夠[4]。李磊等[5]結(jié)合經(jīng)驗(yàn)方法提出了一種混合求解方法，提高了全局搜索水平且排源效果較好。神經(jīng)網(wǎng)絡(luò)算法潛力巨大，但其難點(diǎn)在于需搜集足夠多的優(yōu)秀排源方案對(duì)算法進(jìn)行預(yù)訓(xùn)練[6]。文獻(xiàn)[7-8]提出了一種全局隨機(jī)搜索型方法，針對(duì)輻照裝置劑量不均勻度(DUR)進(jìn)行優(yōu)化排源設(shè)計(jì)并證實(shí)其可行有效，文獻(xiàn)[9-10]在此基礎(chǔ)上對(duì)相關(guān)算法進(jìn)行并行化改進(jìn)并利用GPU加速，可滿足任意裝機(jī)規(guī)模輻照裝置的排源工作。近年來通過與多家單位的交流合作發(fā)現(xiàn)，企業(yè)不僅對(duì)DUR，而且對(duì)射線能量利用率、源棒排布工作量等更多目標(biāo)均有需求，因此本文對(duì)原方法的評(píng)價(jià)數(shù)學(xué)模型、求解算法及代碼均進(jìn)行改進(jìn)，形成一種新型多目標(biāo)優(yōu)化排源方法。

1 排源評(píng)價(jià)數(shù)學(xué)模型

影響60Co γ輻照裝置性能的主要因素及相互關(guān)系如圖1所示[11-13]，優(yōu)化選源、排源是保證和提高60Co γ輻照裝置加工質(zhì)量和加工效率的關(guān)鍵，因目前相關(guān)單位多直接購買成品源棒，源的活度分布選擇余地小，為此排源成為工作核心。相關(guān)單位在評(píng)判排源設(shè)計(jì)方案時(shí)會(huì)綜合考慮DUR、射線能量利用率及排布工作量等，而主流評(píng)價(jià)方法——參考平面法[14-15]難以直接評(píng)判貨箱中劑量分布DUR和射線能量利用率。本文從真實(shí)輻照情景出發(fā)，通過計(jì)算貨箱在輻照室中完成動(dòng)態(tài)輻照后的累積劑量分布，可直接得到上述參量值。

1.1 輻射場(chǎng)分布和貨箱累積劑量分布的計(jì)算方法

(1)

圖1 影響60Co γ輻照裝置性能的主要因素及相互關(guān)系Fig.1 Main factor affecting 60Co gamma irradiation facility and interrelation

常數(shù)，對(duì)60Co，Γ=2.503×10-18C·m2·kg-1·Bq-1·s-1；A為源棒活度，Bq；L為源棒長度，m；a為待計(jì)算空間點(diǎn)到線源垂直距離的絕對(duì)值，m；b為垂足到線源最近點(diǎn)的距離，m；C為校準(zhǔn)因子(測(cè)量值/理論值)，以修訂源棒包殼、劑量計(jì)、輻照裝置附屬機(jī)構(gòu)(如導(dǎo)源管、護(hù)源板(罩))帶來的射線散射和吸收影響，對(duì)已定裝置而言C近似為常數(shù)且歸一化處理后實(shí)際對(duì)當(dāng)前研究無影響。

圖2 源架、貨箱及劑量計(jì)的空間布置Fig.2 Spatial arrangement of source rack, container and dosimeter

(2)

其中：oi為流水線層編號(hào)；si為流水線上貨箱工位編號(hào)。

進(jìn)一步考慮貨箱換面照射，則貨箱中參考點(diǎn)的最終累積劑量Di,j,k可表示為：

(3)

1.2 基于多目標(biāo)的評(píng)價(jià)函數(shù)定義

(4)

2 排源模型的局部復(fù)合求解算法

從評(píng)價(jià)模型可知，排源問題為典型的非線性規(guī)劃問題，相關(guān)研究表明模擬植物生長算法(PGSA)求解排源問題可行有效，該算法通用性強(qiáng)，原則上只需替換評(píng)價(jià)函數(shù)即可[16]。為減少現(xiàn)場(chǎng)布源工作量，本文改進(jìn)了PGSA而形成局部復(fù)合求解算法。

2.1 局部復(fù)合求解的基本原理

PGSA潛在需至少1個(gè)在所有變量維數(shù)上盡量均勻分布的初始方案，文獻(xiàn)[16]表明強(qiáng)魯棒性雖使得初始方案對(duì)最終結(jié)果的優(yōu)化影響不大，但計(jì)算量變化劇烈。為保證初始方案的合理性及降低計(jì)算量，在求解過程中將源架分割為若干區(qū)域，源棒按總活度進(jìn)行比例分配，分別排源后組成1個(gè)次優(yōu)化方案，再進(jìn)行整體優(yōu)化計(jì)算。

2.2 算法的求解流程

算法的求解流程及功能如下。

1) 讀入計(jì)算所需的各種初始數(shù)據(jù)，包括源架、棒位、源棒、貨箱及填充材料、相應(yīng)空間位置及控制算法本身等的系列參數(shù)。

2) 源架可按子源板自然分為N0個(gè)子區(qū)域，或按源架對(duì)稱軸劃分；將源棒按總活度對(duì)應(yīng)分為N0組，接近退役的源棒優(yōu)先分配到一起。

3) 依次將源棒分組排布到各源板子區(qū)域，若全部排布完，則進(jìn)行步驟11；若沒有，則進(jìn)行步驟4。

5) 以xB為起點(diǎn)，以λ(≥1的正整數(shù))為步長，在各維變量上沿坐標(biāo)軸的正負(fù)方向?qū)ふ覞M足約束條件[a,b]的生長點(diǎn)xTemp，若f(xTemp)

6) 求出X內(nèi)目標(biāo)函數(shù)值最小點(diǎn)xmin，若f(xmin)≤f(Xmin)且不重復(fù)，則把該點(diǎn)放入最小生長點(diǎn)集合Xmin。

7) 若生長點(diǎn)集合X容量太大，則以一定的修剪比例對(duì)X內(nèi)的劣勢(shì)生長點(diǎn)進(jìn)行裁減，經(jīng)驗(yàn)值取0.4～0.6為佳；本輪產(chǎn)生的xTemp應(yīng)全部保留從而保持生長點(diǎn)集合X的生物多樣性以防算法早熟而導(dǎo)致優(yōu)化停滯。

8) 若當(dāng)前生長周期內(nèi)的f(xmin)相比f(Xmin)沒有進(jìn)步即視為生長停滯，且停滯次數(shù)超過設(shè)定值，或總的生長周期數(shù)大于設(shè)定值，那么本分組計(jì)算結(jié)束，保存Xmin、f(Xmin)等結(jié)果待合并，進(jìn)入步驟3開始下一分組排源，若不滿足則進(jìn)入步驟9。

9) 計(jì)算X內(nèi)各生長點(diǎn)的形態(tài)素濃度，即各自的狀態(tài)空間概率[16]。

10) 產(chǎn)生隨機(jī)數(shù)β并計(jì)算落入的狀態(tài)空間，將所落區(qū)間對(duì)應(yīng)的生長點(diǎn)記為新生長基點(diǎn)xB，轉(zhuǎn)回步驟5開始新的生長周期。

11) 將所有分組的Xmin合并，完成所有源棒在整體源架上的排放。

12) 若不需進(jìn)一步整體排源，則輸出Xmin形成系列最終排源方案，進(jìn)而精細(xì)計(jì)算得到這些方案在貨箱中的DUR、能量利用率、平均累積劑量等評(píng)價(jià)用結(jié)果。若進(jìn)一步需要整體計(jì)算，則將Xmin作為整體排源用初始生長點(diǎn)X0，依次調(diào)用步驟4～10完成整體優(yōu)化計(jì)算，最終得到系列優(yōu)化排源方案供客戶選用。

2.3 基于GPU的并行加速計(jì)算

本文提出的新評(píng)價(jià)模型相比原模型，其理論復(fù)雜度更高而計(jì)算復(fù)雜度偏低，原模型所耗計(jì)算時(shí)間占總計(jì)算時(shí)間的90%以上[7]，新模型經(jīng)初步測(cè)試計(jì)算占45%左右，需進(jìn)一步提升其計(jì)算效率。文獻(xiàn)[17]中基于GPU通用計(jì)算技術(shù)實(shí)現(xiàn)了PGSA的并行化，在費(fèi)用增加不明顯的前提下使得原方法的計(jì)算效率顯著提高。基于同樣思路，改進(jìn)后的PGSA也進(jìn)行了基于GPU的并行化計(jì)算改進(jìn)。

限于篇幅，本文僅給出當(dāng)前算法中最重要的改進(jìn)部分，即目標(biāo)函數(shù)f的GPU并行化實(shí)現(xiàn)原理，如圖3所示，其中步驟1、2、4、5計(jì)算壓力小，以GPU最大Block和Thread數(shù)執(zhí)行即可，步驟3計(jì)算壓力大，執(zhí)行時(shí)設(shè)計(jì)為每個(gè)Block負(fù)責(zé)計(jì)算生長點(diǎn)集合1個(gè)生長點(diǎn)的目標(biāo)函數(shù)值，Block中每個(gè)Thread計(jì)算1根源棒造成的劑量。

新評(píng)價(jià)模型對(duì)每個(gè)Thread的任務(wù)壓力明顯下降，求解代碼不僅計(jì)算效率顯著提升，且在GPU單精度計(jì)算模式下即可獲得滿意結(jié)果，可避免購買昂貴的專業(yè)雙精度計(jì)算卡。

圖3 PGSA中目標(biāo)函數(shù)f基于GPU的并行化原理Fig.3 GPU parallelization principle of target function f by PGSA

3 實(shí)際應(yīng)用及結(jié)果分析

3.1 相關(guān)裝置參數(shù)

以不同大型輻照裝置的實(shí)踐排源設(shè)計(jì)結(jié)果為例，驗(yàn)證算法的可靠性和有效性，表1列出了輻照裝置的有關(guān)參數(shù)。

表1 輻照裝置參數(shù)Table 1 Parameter of irradiation facility

算法代碼的測(cè)試平臺(tái)主要參數(shù)為：CPU，Xeon E3 1230 V2(數(shù)量1)；GPU，GTX285(數(shù)量1，硬件規(guī)格1.3)；內(nèi)存，DDR3 1333(8 GB)；軟件平臺(tái)，Windows 7 64bit、VS2008、CUDA4.2、C++。

算法的魯棒性較好，參數(shù)限制少。計(jì)算時(shí)源板按每層板塊數(shù)劃分，源棒按總活度等比例劃分并由程序隨機(jī)抽取后均勻布放到棒位上，貨箱中參考點(diǎn)的空間間隔為5 cm，算法修剪比例為50%，編寫參數(shù)輸入文件，利用不同版本的代碼程序計(jì)算，其耗費(fèi)時(shí)間情況列于表2。從表2可知，計(jì)算同等規(guī)模裝置，本文提出的新方法對(duì)應(yīng)代碼的計(jì)算時(shí)間明顯縮短，可顯著縮短排源方案設(shè)計(jì)周期。

表2 不同版本排源代碼的計(jì)算時(shí)間Table 2 Computation time for source pencil deployment of various codes

3.2 結(jié)果及分析

現(xiàn)場(chǎng)部署排源方案后，按國標(biāo)采用重鉻酸銀劑量計(jì)布置測(cè)試用貨箱，表3為測(cè)試貨箱參數(shù)、裝置劑量場(chǎng)的計(jì)算結(jié)果及實(shí)測(cè)結(jié)果。

從表3可知，計(jì)算采取均衡策略，在測(cè)試貨箱充滿模擬貨物前提下，對(duì)于裝置1和裝置2，新模型的理論DUR相比原模型分別增大了7.0%和16.5%，實(shí)測(cè)DUR相比計(jì)算結(jié)果分別增加了4.195%和0.70%，均顯著低于國標(biāo)要求限值而實(shí)際影響不大，同時(shí)能量利用率顯著提高5.72倍和3.02倍。實(shí)測(cè)DUR相比理論值升高，其主要原因是計(jì)算模型相比實(shí)際輻照裝置做了合理簡化，未考慮護(hù)源板、貨箱壁，另外計(jì)算模型以空氣為受照介質(zhì)，而實(shí)測(cè)采用了密度更大的材料，射線的衰減和散射差別較大，IAEA證實(shí)這類問題可采用ρ-DUR關(guān)系曲線修訂而不影響排源[18]，且排源工作核心點(diǎn)是尋找優(yōu)化方案，DUR和能量利用率的值僅用于評(píng)判方案的優(yōu)化程度。

優(yōu)化的排源方案在裝置1、2參考面(圖2)上的等劑量線分布如圖4所示。

表3 測(cè)試貨箱參數(shù)、裝置劑量場(chǎng)的計(jì)算結(jié)果及實(shí)測(cè)結(jié)果Table 3 Test container parameter and result of facility dose field calculation and measurement

注：( )內(nèi)為模型改進(jìn)后的數(shù)據(jù)；[ ]內(nèi)為相對(duì)偏差

圖4 裝置1(a)、2(b)參考面上的劑量分布Fig.4 Dose distribution along reference plane in facility 1 (a) and facility 2 (b)

從圖4可知，排源方案在裝置1、2參考面上的劑量分布很均勻，基本以裝置中心為圓心呈同心圓分布且梯度均勻，這對(duì)降低DUR有利，證實(shí)了排源方案設(shè)計(jì)的可靠性。

4 結(jié)論

本文提出了一種新型多目標(biāo)優(yōu)化排源方法，建立了綜合考慮貨物DUR、射線能量利用率、現(xiàn)場(chǎng)布源工作量等多目標(biāo)參數(shù)的新數(shù)學(xué)模型，設(shè)計(jì)了新型的GPU并行PGSA算法用于模型求解，并對(duì)兩座大型γ輻照裝置進(jìn)行了實(shí)踐排源。計(jì)算和測(cè)量結(jié)果表明：新方法計(jì)算精確度得到保證，對(duì)昂貴專業(yè)計(jì)算卡的依賴度下降；計(jì)算效率提升顯著，計(jì)算時(shí)間從原有方法的h量級(jí)下降到min量級(jí)；新方法的理論和實(shí)測(cè)DUR比原方法增大但仍顯著低于國標(biāo)控制標(biāo)準(zhǔn)，且能量利用率顯著提升，同時(shí)分區(qū)域排源降低了現(xiàn)場(chǎng)布源工作量。檢測(cè)結(jié)果驗(yàn)證了新方法的可行性和有效性。