基于自适应采样算法的芳烃异构化代理模型

谢雨珩，李智，杨明磊，杜文莉

（化学工程联合国家重点实验室，华东理工大学化工过程先进控制和优化技术教育部重点实验室，上海200237）

引 言

异构化装置是芳烃生产过程中的重要环节。随着芳烃反应动力学模型的深入研究[1]，以及计算机过程模拟技术的发展，研究人员可以利用过程模拟软件例如Aspen Hysys 对异构化环节构建机理模型，通过改变温度、压力、流量等操作参数获得关键产物收率，从而对操作参数进行实时的控制和优化。但是机理模型本身计算时间较长，实际过程中直接使用机理模型会耗费大量计算时间，特别是在进行优化时，由于需要反复调用机理模型，所以往往出现优化时间过长、效率过低等问题，很大程度上限制了机理模型在实际优化中的应用。

代理模型的出现，成为解决这一问题的关键技术之一[2]。代理模型是在保证一定精度的情况下，使用计算复杂度较小的数学模型对计算复杂度较大的原始机理模型进行模拟。一方面保证了模型的精度，另一方面大幅提高了计算效率[3]。常用代理模型有：多项式响应面模型[4-5]、径向基函 数 模 型[6-7]、Kriging 模 型[8-9]、支 持 向 量 回 归 模型[10]、多元自适应样条回归模型[11-12]、人工神经网络[13-14]等。

代理模型的建立，主要依赖于采样点的选择，常用的一次采样方法包括均匀采样（uniform design）[15]、拉 丁 超 立 方 采 样（Latin hypercube sampling，LHS）[16]、增量拉丁超立方采样（inherited Latin hypercube sampling）[17]等。但是一次采样方法只是使得样本点在采样空间内尽量均匀分布，而实际过程中，模型的某些区域较为复杂，需要较多采样点，而有些区域实际响应值较为平缓，并不需要过多的采样点。代理模型与实际模型难免有误差，必须采取合理的模型管理策略来保证模型的准确性[18]，所以如何使用较少的原始适应度函数评估次数得到机理模型的最佳采样点分布，进而建立更准确的代理模型成了研究的主要方向。国内外学者提出了很多自适应采样的方法，Wang 等[19]利用基于委员会的主动学习寻找各子代理模型相差最大的点作为新采样点，Eason 等[20]提出了基于交叉验证误差和最小距离的混合自适应采样方法，Wei 等[21]提出了基于最大曲率和最小距离的顺序采样方法（MCMPD），Garud 等[22]提出了基于样本拥挤度和偏离函数的智能采样算法。Müller 等[23]在当前最优点附近生成新的候选点，对这些候选点的距离准则和响应面准则进行加权获得新采样点。与LHS 等一次性采样方法相比，自适应采样算法的计算效率更高，采样点分布更加合理，构建出的代理模型精度也更高。实际应用过程中，由于计算维度往往较高，已有的自适应算法在建模时效率仍然较低，因此需要寻找新的适用于实际大型应用的自适应采样算法。

为了提高构建代理模型的计算效率，减少原始适应度函数评估次数，本文提出了一种新的计算稀疏度的方法，并结合最邻近期望提出了一种基于稀疏度和最邻近期望的自适应代理模型采样算法（adaptive surrogate sampling algorithm-based sparsity and nearest neighbor，ASSA-SNN）。该算法能够在模型变化较大的区域以及未探索区域自适应的分配采样点，从而使得代理模型的精度得到显著的提高。使用多个测试函数验证了该方法的准确性和高效性，并将其应用于芳烃异构化环节建模中，实现对异构化产物的准确预测。

1 Kriging代理模型

Kriging 模型是对空间分布的数据求线性最优、无偏内插估计的一种插值方法，被广泛用于非线性问题代理模型的构造。最早于1951 年由南非地质学家Krige[8]提出，早期主要用来预测矿产储量的分布，之后Sacks 等[24]将Kriging 理论进一步推广到确定性计算机实验的设计和分析领域，并给出了一种较实用的Kriging 模型。Kriging 模型近些年来在化工领域也有较多应用，例如苯乙烯流程优化[25]以及碳二加氢反应器建模[26]。

已知n维的训练样本X=[x1,x2,…,xN]T对应的响应值为Y=[y1,y2,…,yN]T。那么，Kriging 的y(x)由一个全局模型和一个局部偏差构成

式中，f(x)是确定性部分；Z(x)是随机过程，并满足统计特性：期望为0，方差为σ2，协方差为σ2R(xi,xj)。其中，R为相关函数，通常采用Gauss函数

其中，n为X的维数，θk是待定的模型参数。任意x处的Kriging 模型预测值ŷ(x)和预测方差ŝ2表示为

其中，F=[1,1,…,1]T，rT(x)为待预测点x和训练样本点X=[x1,x2,…,xN]T之间的相关向量，定义如下

2 基于稀疏度与最近邻期望的自适应采样算法

2.1 稀疏度

代理模型的建立需要在机理模型上获得初始样本集X，然后根据已有信息找到下一个采样点。新采样点需要同时满足两个条件：能反映复杂区域的响应值和未搜索区域的响应值。所以需要综合考虑新采样点的位置信息和模型响应值信息，前者主要负责全局搜索，后者负责局部搜索。

对于新采样点的位置信息，本文提出了稀疏度的概念。根据新采样点xnew与已有采样点的欧式距离以及每一维数据的大小排序，得到在新采样点在整个搜索空间的稀疏度。已有采样点X={x1,x2, …,xN}，搜索空间的上下限分别为LU={lu1,lu2,…,luD}和LD={ld1,ld2,…,ldD}，D为样本维度。搜索空间内的任意一点xnew的稀疏度Sparsity(xnew)定义如式(7)所示，具体流程如图1所示。

当函数为2维时，样本点的稀疏度是一个面积，以图2 为例，在搜索空间x1,x2∈[-4,4]区间内先随机生成20 个采样点。矩形表示[-2.5,-2.5]、[-2.5,3]和[3,3]三个点的稀疏度。明显看出，[-2.5,3]附近采样点较多，所以稀疏度较小，而[-2.5, -2.5]和[2.5,3]附近点较少，稀疏度比较大。

当输入样本为三维时，稀疏度就可以类比成三维体积，输入样本为四维时稀疏度为四维体积。传统的使用欧式距离计算新采样点与已有采样点的位置信息时，往往只考虑了新采样点与最近已有采样点的位置信息，忽略了它与其他相邻采样点的信息。本文中的稀疏度则综合考虑了新采样点附近已有采样点的位置信息。稀疏度越大，说明这个点附近的采样点较少，应该着重搜索。稀疏度越小，说明这个点附近采样点已经很多，不再需要过多的采样点了。但是，如果只考虑稀疏度就忽略了局部复杂区域的信息，例如最优值附近往往梯度较大，应该需要更多的采样点。所以这里还需要考虑代理模型的响应值。

2.2 最邻近期望

代理模型的建立，不仅需要考虑全局信息，还需要对模型的关键区域进行重点采样，模型的关键区域往往梯度较大，本文采用最邻近期望（nearest neighbor expected,NNE）来近似梯度。

已有采样点X={x1,x2,…,xN}，对应的实际输出为Y={y1,y2,…,yN}，根据X和Y构建的代理模型为ŷ。最邻近期望NNE为

ynearest为样本集X={x1,x2,…,xN}中距离xnew最近的样本点的对应的真实响应值。NNE 越大说明新采样点附近的近似梯度较大，即函数波动较大，需要着重进行采样。传统的交叉验证的主要优点是不需要新的样本点来衡量构建的代理模型的准确性[27]，但是会消耗大量计算时间，而NNE既可以反映函数的局部信息，还可以大量缩短计算时间，尤其是在维数较高的情况下，交叉验证往往需要耗费大量的计算时间，而这时，使用最邻近期望进行局部搜索可以减少大量不必要的计算，使得搜索更有效率。

2.3 自适应采样策略

图1 稀疏度计算流程图Fig.1 Flow chart of sparsity calculation

代理模型建模过程中，最主要的问题是如何利用尽量少的原始函数评估次数，使得样本点尽可能多地分布在具有关键信息的区域[28]，从而缩短总体建模时间。即样本点个数相同的情况下，使获取的关于未知设计空间的信息最大化[29]，提升代理模型精度。

利用自适应采样策略构建代理模型的通用过程步骤如下：

（1）选用适当的取样方法生成初始样本点集，并求得其对应的原始适应度函数。

（2）基于上述初始样本点及其原始适应度函数，建立目标问题的代理模型。

（3）根据已有样本和已构建的代理模型，寻找能反映原始函数关键信息的新采样点，评估其原始适应度函数并添加到样本集中。

（4）判断是否符合终止条件，符合终止条件则输出已有的代理模型，否则返回步骤（2）。

整个流程的关键就是步骤（3）中如何找到反映原始函数关键信息的采样点。在搜索新采样点过程中代理模型应与原始适应度函数一起使用[30]。新采样点一方面要考虑全局信息，即新采样点应该在采样点稀疏的空间内，以保证代理模型在整个搜索空间的泛化性，同时要考虑响应值的局部信息，即波动较大的区域应该有较多的采样点以保证代理模型在此区域的准确性，本文根据稀疏度和最邻近期望的乘积最大来寻找新采样点。

图2 二维稀疏度示意图Fig.2 Diagram of 2D sparsity

稀疏度Sparsity负责控制全局搜索，最邻近期望NNE负责局部关键信息搜索。R为采样空间定义域。初始样本点通过拉丁超立方采样获得，终止条件为新增采样点个数达到设定值。使用模式搜索法[31]求解式(9)。考虑到搜索过程中需要计算欧式距离，但是实际应用中，很多情况下不同维的数据的量纲不同时，欧式距离可能会有较大误差，所以在构建代理模型前，将所有数据根据采样空间的上下限进行归一化，这样保证欧式距离的计算不会受到量纲的影响。具体流程如图3所示。

基于稀疏度与最近邻期望的自适应采样算法（ASSA-SNN）在寻找新采样点时可以有效地权衡全局搜索和局部搜索，从而在迭代过程中寻找到对代理模型精度影响较大的点，不断充实样本集，提高代理模型的精度。

3 测试函数验证

本文首先以二维PK 函数为例，PK 函数表达式为

图3 算法流程图Fig.3 Flowchart of algorithm

初始采样点20，新增采样点80，通过ASSASNN 得到的样本点分布如图4 PK 函数采样点分布情况所示，可以看出在新增的采样点主要集中在上下两个极值点附近，而其余较为平坦的区域采样点较少。使用均方根误差（root mean square error，RMSE）来对算法的精度进行评估，RMSE表达式为

图4 PK函数采样点分布情况Fig.4 Sampling point distribution of PK function

图5 PK函数RMSE变化情况Fig.5 Variation curve of RMSE of PK function

表1 PK函数自适应采样与拉丁超立方采样结果对比Table 1 Comparison of PK function adaptive sampling and Latin hypercube sampling results

本文共选取了8 个测试函数，并与混合自适应采样算法[20]（MASA）、基于最近邻和马氏距离的自适应采样算法[32]（ASA-NNMD）、智能采样算法[22]（SSA）进行对比。为了保证对比的有效性，统一使用Kriging 模型，其中，f(x)为零阶多项式回归函数，相关函数为高斯函数[式(2)]。初始采样点均使用拉丁超立方（LHS）采样获得。由于模型的维数越高意味着模型的复杂程度增加，所以针对不同维数的模型，所选取的初始采样点个数与新增采样点个数不同，具体采样方式如表2所示。

为消除不同初始采样点以及算法寻优时的随机性，每种算法在每个测试函数上运行10 次取平均值以消除随机性带来的影响。最终结果如表3 所示。从结果中可以看出，8 个测试函数中，ASSA-SNN 有6 个测试函数的RMSE 占优，Hartman3 函数的结果也与最小值相差较小，只有在AlpineN.2 函数的结果不如MASA 和SSA。这是由于该函数的极值点过多，局部波动较大，导致搜索到的点过于集中在各个极值附近，影响了整体的效率。

表2 采样方式Table 2 Sampling configurations

表3 各测试函数的平均RMSETable 3 Average RMSE obtained by four sampling approaches for test cases

为了更好地说明算法的计算效率，表4 列出了每种测试函数在不同方法下运行的平均时间。从表格中可以看出由于SSA 使用了逐一交叉验证，MASA 使用了K 折交叉验证，导致计算时间随着维数增加而迅速上升。ASA-NNMD 由于需要通过计算大量矩阵计算马氏距离，以及本身的优化过程较为耗时，导致计算时间也较长。而ASSA-SNN 只有基本的欧式距离计算，通过最邻近期望代替交叉验证，使用模式搜索代替耗时的智能优化算法，达到了精度与效率的双赢，不仅精度有提高，计算时间也大幅缩短。

综上所述，本文提出的ASSA-SNN 在大部分测试函数中的表现优于其他两种算法，并且维数大于3 后ASSA-SNN 计算时间要明显优于其他的算法，函数测试结果说明了ASSA-SNN 的建模精度和建模效率要优于其他算法。

图6 异构化过程流程图Fig.6 Flow chart of isomerization process

表4 各测试函数的平均计算时间/sTable 4 Average running time obtained by four sampling approaches for test cases/s

4 芳烃异构化代理模型构建实例

芳烃异构化反应的目的是将异构体邻二甲苯(OX)、间二甲苯(MX)和乙苯(EB)转化为价值更高的对二甲苯(PX)[33]。异构化过程主要由反应器和分离系统组成，其工艺流程如图6 所示，首先是反应器，其主要功能就是在催化剂作用下进行异构化反应，将含有贫PX 的C8芳烃混合物转化成接近热力学平衡的C8芳烃混合物。其次是分离系统，反应产物在进行分离前需要实现气液分离，将循环氢和反应液送入不同的装置。分离的液体经换热后送入脱轻组分精馏塔，实现轻组分和重组分的分离，分离出来的C8及C8+送到二甲苯精馏单元的分离塔实现C8芳烃与重芳烃分离。

采用化工过程动态模拟软件Aspen Hysys 建立其机理模型。由于机理模型的计算采用序贯模块法逐个装置求解，优化时每次迭代都要将优化后的输入数据重新送入模型中并调用模型求解，采用种群数为20 的遗传算法，迭代次数80 次，总共耗时13059.68 s，耗时较长，不利于优化进行，所以寻求利用代理模型节约计算时间。根据异构化环节的反应机理和经验知识，本文选取8个输入变量：异构化进料、循环氢流量、补充氢流量、异构化反应温度、异构化反应压力、乙苯含量、MX 含量和OX 含量。篇幅所限，本文只选取两个重要的质量指标进行说明：轻烃收率和C8+收率。利用拉丁超立方采样获得80个初始采样点和1000个测试样本点，算法最大迭代次数为320，分别用四种方法建立代理模型并计算它们的RMSE。结果如表5所示，所用建模时间如表6所示。

表5 芳烃异构化重要指标的RMSETable 5 RMSE of important indicators for isomerization of aromatics

图7 给出了异构化代理模型中MX 含量、OX 含量对出口C8+收率的影响关系，可以明显地反映出模型变量之间的复杂映射关系。图8 给出了MX含量、OX 含量样本分布情况，可以看出在各极值附近均有较多样本点分布。从表5 和表6 看出，四种关键性能指标的建模中，ASSA-SNN 的准确度高于其他三种算法，轻烃收率的RMSE 分别降低了16.49%、5.56%、8.20%，轻烃收率的计算时间分别 缩 短 了33.92%、52.41%、89.27%，C8+收 率 的RMSE 分别降低了18.91%、12.20%、7.07%，C8+收率的计算时间分别缩短了35.94%、53.60%、89.61%。以轻烃和C8+收率RMSE 变化曲线如图9 和图10 所示，可以看出ASSA-SNN 的整体RMSE 下降速率要高于其他三种方法，最终的RMSE 也小于其他方法。

表6 芳烃异构化重要指标的建模时间/sTable 6 Modeling time for important indicators of aromatic isomerization/s

图7 异构化代理模型Fig.7 Surrogate model of isomerization

图8 MX含量、OX含量样本分布Fig.8 Sample distribution of MX content and OX content

图9 异构化的轻烃收率RMSE变化曲线Fig.9 Variation curves of RMSE of isomerized light hydrocarbon

图10 异构化C8+收率RMSE变化曲线Fig.10 Variation curves of RMSE of C8+yield with number of iterations

从测试结果中看出，ASSA-SNN 相比于其他三种方法，精度有所提升，能够更加准确地对各关键性能指标进行建模，并且计算时间大幅缩短，这使得整体的建模效率大为提升，为后续的操作变量和操作条件的优化提供便利。

5 结 论

本文提出了一种新的基于Kriging 代理模型的自适应采样算法，该算法通过最大化稀疏度与最邻近期望的乘积来确定新采样点的位置，从而使得模型能够根据实际函数自适应的分配采样点位置。相比其他自适应采样算法，使用稀疏度可以更好地判断新采样点与已有采样点的位置关系，从而更好地进行全局搜索，用最邻近期望代替交叉验证，提高了局部搜索效率，使用模式搜索解决最优问题，极大地节约了计算时间，提升了寻优效率。将该算法与其他的自适应采样算法在8个测试函数上进行了验证，表明在相同采样点个数的情况下，ASSASNN 得到的采样点分布能够有效提高建立的代理模型精度并且缩短建模时间。芳烃异构化实例也证明了算法的有效性，为后续的异构化环节操作条件的优化提供了便利。