基于NSGA-II的晶体学光束线多目标自动优化

张 丁 吴盈锋 何迎花 刘 科 汪启胜1, 何建华

1（中国科学院上海应用物理研究所 上海201800）2（中国科学院大学 北京100049）3（中国科学院上海高等研究院 上海201210）4（武汉大学高等研究院 武汉430072）

生物大分子晶体学是在原子分辨率水平上解析生物大分子结构最主要的技术手段之一。相比于常规实验室的X射线光源，同步辐射的高通量性使得之前绝大多数很难或者不能进行的实验都能成功开展下去，基于同步辐射的生物大分子晶体学也受益颇多。除了光源带来的诸多优越性外，随着实验硬件、软件的突破以及实验方法的创新和发展，同步辐射生物大分子晶体学也飞快发展成熟［1］。但随着生命科学领域的迅猛发展，海量的新蛋白质结构的分析测定需要依靠同步辐射生物大分子晶体学的实验方法，对同步辐射生物大分子晶体学光束线站的需求日益增加［2］，光束线必须要能够极其可靠高效地进行高通量实验来满足这种需求，重要的手段就是实现光束线的自动化运行，建立自动化程度高的集成系统已经成为几乎所有同步辐射生物大分子晶体学实验站发展的重要工作［3］。

一个完整的生物大分子晶体学实验步骤相对繁复，包括复杂的光束线优化、样品安装和准直、实验条件确定、样品衍射情况检查、实验策略的制定优化、数据采集以及后续的数据处理与分析。作为整个实验流程的第一个环节，实验前光束线上光学元件的配置优化对于用户实验是十分重要的。在此之前，已经有一些实验室提出了光束线自动优化的方案，Pugliese等［4］将模糊逻辑应用到光束线自动优化上；Olivier等［5］使用波前分析法作为优化策略，但这些方案都不能保证全局优化，直到2015年杜永华［6］将遗传算法成功应用到新加坡光源光束线全局自动优化上，随后他又采用差分进化算法，并使用LabVIEW开发出光束线自动优化程序AI-BL1.0［7］。上海光源的衍射线站借鉴新加坡光源的成功经验，将该程序应用到上海光源基于EPICS系统（Experiment Physics and Industrial Control System）的控制平台上［8］。上述案例中的优化目标都是样品位置的光通量，然而光束性能指标除了光通量之外，还包括光斑尺寸、光束位置等，因此光束线的优化问题是一个典型的多目标优化问题。非支配排序遗传算法II（Nondominated Sorting Genetic Algorithm II，NSGA-II）是一种建立在Pareto最优解理论上的多目标全局优化算法，非常适用于光束线复杂系统的全局优化问题上。同时上海光源晶体学光束线站的运动控制和状态数据获取都是基于EPICS实现的，Python和EPICS的接口已经被各光源和光束线广泛应用，因此利用Python开发优化程序，不仅能使该程序方便应用到晶体学光束线的自动优化上，同时也能使程序有效地推广和应用到其他线站的优化上。

本文提出一种基于NSGA-II的多目标优化策略，并将其应用于上海光源P2生物防护蛋白质晶体学线站（BL10U2）［9］的自动优化上，用Python实现的优化程序在BL10U2上完成了初步测试，并达到了预期目标。

1 建立优化模型

同步辐射光的性能指标包括光子能量、光通量、光斑大小、光束位置等。良好的光束质量是获得可靠实验数据的基本保证，也是提高线站运行效率的前提条件。光束线从电子储存环引出同步辐射光，再按实验要求对其进行调制，并将其传输到实验站进行实验［10］。狭缝、单色器、聚焦镜等光束线设备对同步光进行准直、单色化、聚焦和传输，将特定性能的同步光传输到实验站。光束线上每个部件的姿态都会影响实验站样品处的同步光性能，因此，光束线优化问题是一个典型的多目标全局优化问题。

1.1 多目标优化问题的数学模型

实际问题中大都具有多个目标需要同时满足，即在同一个问题模型中同时存在几个非线性目标，这些目标函数同时需要进行优化处理，并且通常是互相冲突的，此类问题被称为多目标优化问题（Multi-objective Optimization Problem，MOP）［11］。多目标优化问题的数学描述由决策向量、目标函数、约束条件组成，一般表述如下：

式中：x为决策向量；y表示目标函数，也叫目标向量；X为决策空间，由决策向量构成；Y是目标空间，由目标向量构成；e（x）≤0是约束条件，表示决策向量的可行取值范围。

MOP由多个目标函数组成，它会产生多个解，这些解不能像单目标优化问题那样直接比较出它们之间的优劣，被称为目标函数的非支配解（Nondominated Solutions），也叫做Pareto最优解（Pareto optimal solutions）。多个Pareto最优解组成Pareto最优解集，求得Pareto最优解集就是MOP要达到的最终目的，决策者根据实际情况和自身需要从解集中选取合适的Pareto最优解。以下为多目标优化的重要定义：

定义1：Pareto支配

对于最小化多目标优化问题，对于n个目标函数fi（x），i=1，…，n，任意给定两个决策向量xa、xb，如果有以下两个条件成立，则称xa支配xb。

1）对于∀i∈1，2，…，n，都有fi（xa）≤fi（xb）成立。

2）∃i∈1，2，…，n，有fi（xa）＜fi（xb）成立。

定义2：Pareto最优解

如果一个解x*不被其他任意一个解支配，则被称之Pareto最优解。

定义3：Pareto最优解集

多目标优化问题中所有Pareto最优解组成的解集称为Pareto最优解集（Pareto Set）。

定义4：Pareto前沿

Pareto Set中每个解对应的目标向量组成的集合称之为Pareto前沿（Pareto Front），简称为PF。

1.2 光束线优化模型

光束线的优化主要是通过改变决定光学元件姿态的运动电机位置，来调节光束性能（光通量、光斑位置等）。对照多目标优化问题数学模型来看，一组电机位置对应着一个决策向量，该组电机位置决定的光束性能则是对应的目标向量，约束条件则是各个电机的行程范围。光束线优化问题可以抽象成寻找多个目标函数的极值问题，多个目标函数无法兼顾，最终优化得到的结果是这些目标函数的Pareto最优解集，也就是若干组非支配的电机位置，决策者根据实际需要选择合适的电机位置。

2 算法概述

2.1 算法简介

1994年，Deb等［12］在遗传算法的基础上对其进行了改进，从而提出了NSGA，它是一种基于Pareto最优解的多目标优化算法，该算法对种群中的个体进行分层，处在同层的非支配个体共享一个虚拟适应度值，并采用共享函数法来保证个体多样性。尽管NSGA突破了传统的遗传算法，可以求解任意多的目标函数问题，并且在实际生产生活中针对于多目标求解问题得到了广泛的应用，但是它还是存在一些问题，主要体现在支配排序的计算复杂度高、缺少精英机制、需要指定共享参数这三个方面。2000年，Deb又在NSGA的基础上提出了一种带有精英机制的快速非支配排序遗传算法（NSGA-II）［13］，该算法采用了快速非支配排序算法和拥挤比较算子，计算复杂度比NSGA极大地降低；通过计算拥挤距离和定义拥挤偏序，使得解的分布性更佳，保持了种群的多样性；引入了精英保留策略，有效防止了精英个体的丢失。

2.2 算法理论

NSGA-II的关键点在于快速非支配排序、拥挤度和拥挤算子以及精英策略，以下是它们简单的介绍［14］。

2.2.1快速非支配排序

找出初始种群中的Pareto最优解集，将它们的Pareto等级设置为1，然后将它们从种群中删除，在余下种群中Pareto最优解集，将它们的Pareto等级设置为2，以此类推，通过支配关系对整个种群进行排序，得到所有个体的Pareto等级。

2.2.2拥挤度和拥挤算子

遗传算法有自动收敛的性质，设置拥挤度nd能够使同一Pareto等级中的个体相互间分开，从而保证了种群的多样性。先令同一非支配层里所有个体的nd为0，然后对每一个目标函数fm，根据目标函数值的大小对个体进行排序，排序后两个边界个体的拥挤度nd为无穷大，其余个体的拥挤度：

累加个体在所有目标函数的拥挤度得到该个体最终的拥挤度。

利用每个个体的Pareto等级和拥挤度可以得到其拥挤比较算子，概括来说就是：两个不同Pareto等级的个体，Pareto等级更低的个体更优；两个相同Pareto等级的个体，拥挤度更大的个体更优。

2.2.3精英策略

NSGA-II采用精英策略将优良个体遗传到下一代，具体操作为：在父代种群Pt基础上通过遗传算子来产生子代种群Ct，将父子种群混合组成新种群Rt，对Rt进行非支配排序，按Pareto等级由低到高依次将非支配集F1，F2，…，Fm放入新的父代种群Pt+1中，直到Pt+1的大小超出原有种群大小，按拥挤度从小到大从Pt+1中剔除Fm中的个体直到Pt+1的大小和原有种群大小相等。

2.3 算法流程

算法流程如图1所示。NSGA-II具体过程如下：

图1 NSGA-II流程图Fig.1 Flowchart of NSGA-II

1）初始化种群。随机产生规模大小为Np的初始父代种群：

式中：rand（0，1）为0到1的随机数；n为决策向量里的变量个数；uj、lj分别为变量的上、下限。

2）快速非支配排序。

3）计算拥挤度。

4）交叉、变异得到子代种群。

5）混合父代种群和子代种群，确定混合种群中个体的Pareto等级和拥挤度。

6）采用精英策略，选择混合种群优秀的个体生成新的父代种群。

7）重复步骤4）～6），直到达到设定的进化代数，即可获得Pareto最优解集。

3 程序设计

NSGA-II中种群个体（决策向量）的Pareto等级和拥挤度都是由其对应的目标向量决定的，因此，在算法实现中，必须保证决策向量和目标向量的是一一对应的，也就是要保证一组电机位置和由其决定的光束状态信息的对应关系。上海光源晶体学光束线站的运动控制和光束状态信息采集都是基于EPICS系统［15］实现的，该系统的核心是IOC软件，IOC包括分布式运行数据库（Distributed Runtime Database）和输入输出支持模块。EPICS利用IOC运行数据库中的记录完成对输入和输出信号量的控制，记录值的标识称为过程变量（Process Variable，PV），记录值在数据交换中充当基本单元，在客户端与服务器端之间建立通信连接。客户端与服务器端模型建立通道访问协议（基于TCP/IP协议），并提供相应的应用接口函数库，以供数据读写、连接监控和访问监控等服务。

电机的运动控制采用EPICS中的Motor模块［16］，一个Motor记录代表一个电机实体，记录中的域对应着电机属性，其中域VAL为控制电机运动的位置，域RBV为电机位置的回读值。光通量和光束位置的获取和处理则是由QBPM［17］和数字信号处理器组合完成。QBPM探测电极上获取到光电流信号，然后模拟信号被转换为数字信号后提供给IOC分析，从而得到最终光通量和光束位置。因此，电机运动、电机位置获取以及光束状态信息的获取都可以通过操作PV来实现。

NSGA-II的实现由Python完成，优化程序通过PyEpics接口与EPICS连接，PyEpics不仅有对EPICS自带接口函数的封装，还有专门的电机控制模块和PV操作模块，使用者可以自由选择最方便的模块进行操作。优化程序接入EPICS后，可以控制电机运动，获取电机位置回读值和光束信息，并在此基础上执行算法流程。自动优化的软硬件结构如图2所示。

图2 光束线自动优化软硬件结构Fig.2 Software and hardware structure of automatic beamline optimization

程序开始后，在设定的电机行程范围内随机生成初始父代种群，根据种群个体使各个电机运动到相应位置，然后获取此时的光通量和光束位置作为该个体的目标向量，确定种群个体的Pareto等级和拥挤度，父代种群通过遗传操作得到子代种群，根据子代种群信息控制电机运动和获取相应光束状态信息，混合父子种群并确定该混合种群中个体的Pareto等级和拥挤度，随后依据精英策略，在混合种群中选择个体生成新的父代种群，循环更新父代种群，直到达到设定的进化代数，迭代结束，从得到的Pareto最优解集选出所需要的最终结果。图3为优化程序流程图。

图3 自动优化程序流程图Fig.3 Flowchart of automatic optimization program

4 在线测试

本程序部署在上海光源P2生物防护蛋白质晶体学光束线（BL10U2）。该线站具有P2生物安全防护功能，可满足我国对传染性病毒晶体结构研究的需要，同时它也具有与其他的生物大分子晶体学光束线相同的功能，是上海光源里具有代表性的生物大分子晶体学光束线。光束线主要设备包括白光狭缝（Slit1）、水平偏转镜（M1）、双平晶液氮冷却单色器（DCM）、水平预聚焦镜（M2）、次级光源狭缝（Slit2）、垂直聚焦镜（M3）以及水平聚焦镜（M4）；除了以上光学设备，还包括束线调试及诊断必不可少的荧光靶、位置灵敏探测器及光强检测等设备。布局如图4所示。

图4 上海光源BL10U2布局示意图Fig.4 Layout of the optical equipment of BL10U2

由于P2线站还处于调试阶段，测试优化目标选为次级光源点的光通量和光束位置，水平预聚焦镜M2为被优化设备，水平预聚焦镜采用椭圆柱面反射镜，对光束进行水平方向进行聚焦，形成水平方向的一级聚焦光斑。该镜子涉及到的运动方向和对应控制电机如表1所示。次级光源点的光束性能状态则是由水平预聚焦镜后的QBPM获取。

表1 水平预聚焦镜M2的运动方向及相应控制电机Table 1 The movement directions and corresponding control motors of M2

先以不同种群大小（15、20、25）进行了3次测试，进化代数设置为50，3个种群找到的Pareto最优解的个数分别为9、13、14。种群数越大，能找到的Pareto解也越多，但后两组的差距很小，可知当种群规模到达20后，继续增大种群规模对Pareto最优解集的影响较小。在此基础上，以种群大小为20、进化代数为80增加了1次测试，得到Pareto解的个数为20。4次测试的交叉概率为0.9，变异概率为0.2，测试结果列于表2中，4组测试的Pareto前沿如图5所示。

表2 测试条件及结果Table 2 Test conditions and results of the experiment

测试结果表明：1）优化的目标是同时寻找光通量的极大值和光束实际位置到参考位置偏差的极小值。在图5中，横坐标以QBPM总探测电流来代表光通量，纵坐标是光束位置偏差，从图5中可以看出，两个优化目标之间存在着制约关系。本次测试采取折中的方法，程序预先设定的优化目标以光束位置为标准，将目标位置定为距参考位置10μm的位置，程序运行结束后会自动从Pareto最优解集中选取出离目标位置最近的位置和相应的光通量作为优化结果，4次测试的优化结果如表3所示，4次测试的光通量都在1.8×10-6A左右的强度；2）在此之前上海光源晶体学光束线的调光工作主要是由线站工程师手动完成的，每次只能对单个运动电机进行优化，即对每个电机在其行程范围内依次调整位置，直到光束状态满足实验要求，这种调光模式效率非常低，整个优化过程少则几小时，多则数日。本文开发的优化程序能同时对多个电机进行全局优化，所有测试的最长优化时间在30 min左右，相较于手动调节，优化效率大大提高；3）图5中，当进化代数为50时，种群大小为15得到的Pareto最优解个数较少且散乱分布；增大种群规模到20后，数据点增多，但其分布有所起伏，不太均匀；继续增大种群规模到25，数据点虽然没有增加太多，但其分布更加均匀；保持种群规模20不变，增加进化代数到80，得到数据点最多，且其分布也十分均匀和紧密。增大种群规模和进化代数可以得到数量更多和质量更优的Pareto最优解集，也能更好地满足多样的实验需求。

表3 测试优化结果Table 3 Optimization results of tests

图5 4组测试的Pareto前沿(a)种群大小15，进化代数50，(b)种群大小20，进化代数50，(c)种群大小25，进化代数50，(d)种群大小20，进化代数80Fig.5 Pareto front of tests (a)Population size 15,generations 50,(b)Population size 20,generations 50,(c)Population size 25,generations 50,(d)Population size 20,generations 80

5 结语

本文设计并开发了基于NSGA-II的光束线多目标自动优化程序，进一步完善了上海光源光束线自动优化方案。基于Python实现的优化程序不仅能适用于上海光源生物大分子晶体学线站，还具备向上海光源其他线站推广的潜能。自动优化程序结合上海光源P2生物防护蛋白质晶体学线站软硬件条件，以次级光源点的光通量和光束位置为优化目标完成了部署与测试。测试结果表明：该优化程序能准确找到预定的优化目标。测试耗时均在30 min以内，相较于之前的单目标优化方案，基于NSGA-II的自动优化程序不仅保证了自动优化效率，还进一步提高了光束线智能优化的能力。

测试结果还表明：为了获得更好的Pareto最优解集，则需要花费更多的优化时间。后续为了实现真正的全局优化，待优化的光学元件会增多，相应电机也会增多，电机运动时间将会大大增加，为了获得更准确的光束信息，采样时间也会增加，所以为了确保优化效果，相关硬件性能的提高也是非常有必要的。同时测试时P2线站还处于调试阶段，优化程序得到了初步成功的结果，但仍需在P2线站稳定运行状态下做进一步的验证。