李晓琪，郑友琦，杜夏楠，王永平

(西安交通大学 核科学与技术学院，陕西 西安 710049)

与压水堆相比，快堆的运行经验相对匮乏。与此同时快堆的设计目标除寿期更长、功率峰因子更小以外还包括获得更大的辐照量、更深的燃料增殖或焚烧等。在工程实践中，还需要考虑诸多约束条件，如热工、瞬态安全准则等。快堆优化设计归根结底是一个多变量、多目标、有约束的寻优问题。考虑到庞大的搜索空间，枚举法并不可行；考虑到问题的多峰特性，梯度下降法等数学模型亦不可取，因此需要引入智能优化算法[1-8]。在本文中使用的遗传算法作为经典的智能优化算法，适合于解决这类多输入、多目标、带约束的问题。

反应堆燃料管理包括两种类型的问题：1) 关注组件的制造参数以及燃料装量等的“Out-of-core”问题；2) 对现成的组件关注布置方案的“In-core”问题。对这两种问题进行优化的方法分为两类，分别是确定性方法和随机性方法。早年关于压水堆的燃料管理优化方法中，确定性方法发展相对较早，如对燃料和毒物进行脱耦的哈林脱耦方法[9]，通过在中子扩散方程中引入0～1变量，将换料方案优化转化为数学规划问题的直接优化方法[1,10]等。随机性方法多为启发式，包括模拟微观粒子热运动特性实现优化的模拟退火算法、对搜索方向进行标记进而提高搜索效率的禁忌搜索算法[11]、模拟动物群落通过标记信息寻找目标的粒子群算法[12]，以及模拟生物种群进化规律的遗传算法等。上述随机性算法通过从某种自然或社会现象中获得启发以实现对方案的定向搜索，具有较强的通用性，基本不依赖于物理背景本身，适用于本文要解决的快堆优化设计问题。近年来蓬勃发展的人工智能推动了神经网络在堆芯换料优化及设计优化领域中的应用，通过样本训练代理模型进行堆芯参数预测以实现优化成为堆芯优化研究的一重要方向[13-14]。

遗传算法最早由从事自适应系统研究的Professor Holland团队提出。1967年Holland教授的学生Bagley首次在论文中应用了遗传算法，并使用了包括选择、交叉、变异、占优、倒置在内的遗传算子。1992年，Poon与Parks最早将遗传算法用于反应堆燃料管理优化问题。在早年用遗传算法进行反应堆优化设计的研究中，优化目标基本为keff和功率峰因子[15]。

使用遗传算法进行换料优化与堆芯设计的研究中，通常从随机或人为给定的一组初始方案出发，通过不断迭代改变搜索方向进行寻优。如对每个组件进行独立编号，以编号的排列作为方案的基因[2]，并配合使用改进的算子获得新基因，从而避免无效方案的产生[16]。早期的研究大多对装载方案的十进制排列进行二进制编码，以模拟遗传现象的本质[17]。对于规模较大并涉及到换料的问题，文献[18]提出一种对装载方案基因分段，分别表示组件类别与循环次数的方式。由于搜索空间对于整个寻优过程是未知的，初始方案的偏移可能导致搜素过程陷入局部最优，因此需要慎重确定初始方案。本文提供一种新的编码方案，对搜索空间中含有的所有个体进行编码，在优化的起始在多个自变量维度上均匀抽取初始方案，从而一定程度避免陷入局部最优。

1 理论与优化框架开发

1.1 DAKOTA/SARAX耦合

在本文的优化框架中，使用了开源工具箱DAKOTA中的多目标遗传算法模块MOGA。DAKOTA工具箱由桑迪亚国家实验室开发，包含优化、敏感性分析等功能。优化模块MOGA可作为黑箱使用，由中子学计算程序SARAX[19-20]对个体(方案)进行分析(计算)并反馈给DAKOTA作为遗传操作的凭据。SARAX是由西安交通大学开发的快堆中子学计算程序，适用于快堆的临界和次临界计算。使用Python脚本将两个外部程序进行耦合，功能是将DAKOTA提供的个体(方案)翻译为堆芯计算的参数或建模作为输入卡片传递给SARAX，并调用SARAX进行稳态计算、燃耗计算、反应性计算等，在计算结束后将计算结果处理成遗传算法可以读取的“适应度”即衡量方案优劣的指标返还给DAKOTA进行下一步操作。

SARAX是一套中子学计算程序，包括用于截面产生的模块Tulip、堆芯稳态计算模块Lavender及瞬态模块Daisy。本文的工作集中在使用预先准备的截面实现堆芯层面的参数或装料的优化方案搜索，使用堆芯稳态计算模块Lavender。

在燃料管理与堆芯设计中使用遗传算法，流程如图1所示，分为以下步骤：1) 遗传算法模块产生第1代随机个体(方案)的基因文件，翻译成堆芯程序输入后传递给堆芯计算模块；2) 堆芯计算的结果处理成适应度函数返回给遗传算法模块；3) 遗传算法模块进行选择、交叉、变异操作，形成包含新个体(方案)的种群；4) 重复步骤3和4直到种群收敛。

图1 堆芯设计和燃料管理中的遗传算法流程Fig.1 Genetic algorithm flowchart in design of reactor core and fuel management

该流程在本框架中的数据流动如图2所示。

图2 优化框架的流程图Fig.2 Flowchart of optimal frame

1.2 编码策略

遗传算法模仿生物进化的过程，其中包含一个重要的步骤即基因与设计参数或装载方案之间的编码/解码过程。将二进制或十进制的编码翻译成连续的堆芯参数相对简单，而翻译成装载方案则需要一定的策略。

本文提出了一种直接实现十进制基因文件与装载排列或组合方案之间的一一映射的编码/解码思路，并用于此框架，具体实现如下。

1) 从十进制到全排列

在这种装载策略下，每一个组件具有独一无二且连续的编号，组件数量与装载位数量一致，如图3所示。使用一种映射算法“康托展开”，实现十进制整数与全排列之间的一一映射。表1列出对3个数字1/2/3进行全排列下排列方式与序号之间的对应关系。此过程可以通过式(1)实现。

图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为该全排列的十进制序列号(从0开始)，对于本框架相当于基因从十进制基因到全排列的过程，需要将表1右侧的十进制整数(即基因)翻译成左侧的排列方案，进而作为组件的布置，此过程可通过康托展开的逆过程实现。

表1 十进制整数与全排列的对应Table 1 Correspondence between decimal integer and full permutation

2) 从十进制整数到组合数

对于大多数装载策略，同一类型的组件有多个，都需要参入到排列中。如果直接套用上述的全排列思路，则会产生多个基因对应本质上同一种排列方式的问题，导致大量重复方案产生，进而影响优化效果。对于这种组合问题，本文提供一种编码方案：(1) 先针对一种类型的组件，挑选相应数量的装载位，此过程是1个组合问题，其中需要进行1次十进制整数与组合数之间的映射，即将十进制整数转换为一系列组合数，而这个组合数就是装载该种类型组件的装载位，此十进制数是基因的一段；(2) 对另一种类型的组件，在剩下的装载位上重复过程1，以此类推；(3) 当只剩下一种类型的组件，则直接填充剩下的装载位。

组合问题的编码过程如图4所示。

图4 组合问题的编码Fig.4 Coding of combinatorial problem

考虑到如图5所示的情况，待装载的组件总数比装载位多。解决方案是在原基因前增加一段，用于决定哪些组件不参与此次排列。

图5 考虑的后备组件Fig.5 Considering of reserved assembly

由于快堆不需要考虑1/4对称或1/8对称的情况，上述编码方案基本可以应对快堆中的编码问题，并且具有以下优势：(1) 有效方案与指定范围内连续变化的基因一一对应，搜索空间中每个个体都有1个独一无二的编号，利用这些编号可在优化初始进行均匀抽样，无需提供初始方案，如图6所示；(2)基因型的微小变化对应表现型的微小变化，该特点对于随机性优化算法是重要的评估标准。

a——传统遗传算法初值选取；b——本文改进方法初值选取图6 初始方案的分布Fig.6 Distribution of initial scheme

1.3 多目标处理

遗传算法广泛用于多目标优化问题，包括反应堆设计和燃料管理。由于优化目标不止1个，如何让算法综合考虑每个优化目标以评估方案的优劣，直接影响到优化效果。处理方式分为两种类型：将多目标与约束条件按权重组合成1个适应度函数作为评估准则，或直接使用多目标。本文搭建的框架中，使用了DAKOTA中内置的多目标适应度处理功能domination_count和layer_rank。这两种功能对被评估的个体在每个优化目标维度上进行比较，将“占优”数量作为评估个体优劣程度的指标。这与传统的将多目标用1组权重线性相加合成1个适应度函数相比，不需要人为给定权重，使用更简单且不易因为权重选取不合适导致对多目标评估的偏差，是遗传算法适应度函数处理当下发展的一个主流做法。

1.4 约束条件处理

除优化目标，方案优化问题中可能需要设置约束条件，对不满足约束条件的方案直接淘汰。由于遗传算法通过适应度函数对方案进行评估，因此约束条件也将通过作用于方案的适应度函数来实现，称为罚函数，即对不满足某个约束条件的方案在其适应度函数后加上1个数值，降低该方案在后续的选择操作中被选中的概率。对于多目标问题，当采用将多目标与约束条件按权重组合成1个适应度函数的权重法处理方式时，将罚函数按式(2)所示添加在适应度函数中即可。本文中直接使用多目标优化，则将罚函数依次添加在不符合某一条件的方案的每一个优化目标函数值内。

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

(2)

式中：ωi为优化目标i的权重；fi为优化目标i在当前方案下的值；P为罚函数，当前方案不满足某个约束条件时P不为0，都满足则取为0。

2 优化框架验证

2.1 旅行商问题

在对反应堆优化设计问题进行框架验证之前，先通过1个旅行商问题(TSP)验证编码方案是否合理以及DAKOTA内置选项使用是否正确。

此处设计了1个简单的TSP：若干城市以1个单位长度为间隔排列成1圈，则最短路径是依次环绕这些城市1圈，如图7所示。该模型的优点是容易确定最短路径与对应的距离，并且易于增加或减小问题的规模。

图7 简化旅行商问题Fig.7 Simplified TSP problem

本文考虑了12个城市，即进行12个整数的全排列问题搜索，解空间的大小为479 001 600个方案，最优方案数目为24个，随机生成初始方案。适应度即为当前方案下遍历12个城市的路径长度。遗传算法参数设置列于表2。方案搜索过程中，每代最短路径变化如图8所示。DAKOTA内置小生境选项，即根据当前种群个体的分散或聚集程度，通过调整个体适应度，避免个体分布过于集中，以提高种群的基因多样性从而避免陷入局部最优。测试表明，本文给出的编码方案可在无初值的条件下，可接受计算时间内，配合小生境的使用，对12个城市的TSP收敛到全局最优，从而验证了编码方案的合理性，以及DAKOTA遗传算法参数的选择的正确性。

表2 TSP的遗传算法参数设置Table 2 Parameter setting in TSP

图8 城市数为12的TSP收敛Fig.8 Convergence of 12 cities TSP

2.2 基于ABTR的优化与枚举法

本文基于先进燃烧实验堆ABTR的燃料装载优化，通过与同一问题的枚举结果对比，验证本文提出的编码方式对于十进制与组合数转换的处理。ABTR是一个用于验证嬗变技术的概念堆，堆中有两种不同富集度的燃料组件，其中高富集度燃料分布在外圈，低富集度燃料分布在内圈，周期为4个月，批数分别为12和15。此外还有一些用于进行材料辐照实验的测试组件分布在活性区中。堆芯活性区的1/3示意图如图9所示。

图9 ABTR原始活性区布置Fig.9 Original active zone layout of ABTR

对给定数量的3种燃料组件，重新布置其排列方式，目标是对keff和功率峰因子扰动最小的前提下在测试组件中获得最大的中子通量。以1/6堆为对称条件，解空间有1 260个方案，通过枚举获得每个可能方案的3个参数，作为判断优化搜索是否寻找到最优解的依据。在该问题下，基因文件具有两段信息，分别用于决定外圈高富集度燃料(未必仍然排在外圈)和内圈低富集度燃料在堆芯中的位置。当两种燃料的位置确定后，剩余位置即为测试组件的位置。对3个优化目标的函数值采用1∶1∶1的权重加和为1个适应度函数，其中keff和功率峰因子扰动量的绝对值直接取为优化目标函数值，测试组件的中子通量取负作为优化目标函数值，如式(3)所示。

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

(3)

式中：ω1、ω2、ω3为权重，ω1=ω2=ω3=1；ppf和ppf0分别为当前方案下的功率峰因子和原始方案的功率峰因子；φtest为测试组件处的中子通量。

经过对353个个体的计算，获得了一最终优化结果，堆芯活性区布置示意图如图10所示。此外对整个优化空间的1 260个方案依次进行中子学计算，并记录下keff、功率峰因子、测试组件处的中子通量密度，与最终优化结果一起绘制成三维示意图，如图11所示。用kh表示功率峰因子，T_flux表示测试组件处的中子通量密度，红色点表示搜索出的最优方案，在二维坐标上的投影如图12所示。

图10 优化框架提供的最优方案Fig.10 The best scheme provided by optimal frame

图11 最优方案在完整解空间中的分布Fig.11 Distribution of optimal schemein whole solution space

图12 最优方案在二维坐标的投影Fig.12 Mapping of optimal scheme in two-dimensional coordinate

3 优化框架在实际问题中的使用

引入中国实验快堆(CEFR)作为研究目标。CEFR是中国第一座钠冷实验快堆，运行周期为80 d，活性区布置如图13所示，是一个用于各种燃料、材料辐照实验与放射性同位素生产的优良平台[21]。快堆的用途之一是同位素焚烧和生产，可以用于焚烧随着压水堆的运行产生的大量对环境有危害的237Np。237Np可通过嬗变转换成能制作同位素电源的238Pu。基于CEFR堆芯，装载含有237Np的组件以实现238Pu在1个运行周期内产量最大化为目标，同时保证后备反应性、功率形状分布尽可能不偏离原始堆芯。237Np以NpO2的形式装载到堆芯中，包括在燃料活性区和增殖区的装载，如图14所示。通过搜索在哪个组件的活性区或增殖区中装载以及活性区中的装载量，实现表3所列的优化目标与约束条件。其中对功率波动的优化与约束是由于考虑Np的装载影响堆芯的中子价值，进而影响控制棒组的价值，而7%的功率波动限值是考虑了CEFR的反应性控制系统设计误差限给出的经验值。

图13 CEFR的活性区布置Fig.13 Layout of active zone in CEFR

表3 优化目标与约束Table 3 Optimization objective and constraint

图14 装载NpO2的组件Fig.14 Subassembly with NpO2

此处采用直接的多目标优化方法，故不构造适应度函数，而是由DAKOTA直接对每个方案的每个优化目标之间进行比较。

使用DAKOTA内部提供的rank_method作为选择算子，个体在keff和功率波动两个维度上的分布如图15所示，其中红框内是符合约束条件的个体。显而易见，符合条件的个体在这种搜索条件下并没有得到很好的保存，原因是在rank_method内，以罚函数形式出现的约束条件并不能发挥作用。随后将选择算子改为传统的权重函数加和的形式，获得的个体在keff和功率波动两个维度上的分布如图16所示。其中最优方案的堆芯布置如图17所示，对应的功率波动分布如图18所示。在该方案下，237Np的转换率为0.616%，80 d238Pu产量为0.549 kg。

图15 rank_method个体的keff和功率波动Fig.15 keff and power fluctuation of individual in rank_method

图16 权重法个体的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 总结

本文搭建了一个用于快堆设计的优化框架，可实现读入一定范围的决策变量、优化目标以及约束条件。搜索过程基于初始均匀抽样，对应提出一套对整个解空间进行编码的编码方式，可以对排列问题、组合问题进行连续编码。基于该框架，验证了简化的TSP、ABTR装料问题，从而证明编码方式的合理性和框架的可行性。最后通过一个基于CEFR实现Np-Pu生产优化的算例示范了优化框架在工程实际中的具体应用。