摘要
基于代理模型的优化方法因其高效的寻优能力为高精度仿真模型在优化设计中的应用提供了一种有效的技术途径。针对优化问题中高耗时黑箱约束处理问题，提出了一种基于改进可行性规则的多约束自适应采样方法，建立了精英档案驱动的非精确搜索方法和ε-约束保持的伪可行域构造方法，通过迭代过程中动态缩放可行域以接受高质量非可行样本，增强了算法对可行域边界的探索能力，提高了基于代理模型的寻优能力。对进化计算会议约束优化标准函数的仿真计算结果表明，相较于现有方法，ε-约束保持优化方法可以有效解决多约束代理模型优化问题。对固体火箭发动机后翼柱装药设计结果表明算法具有应用在复杂工程问题中的潜力。
Abstract
The surrogate model-based optimization method provides an effective technical approach for the application of high-precision simulation models in optimal design due to its efficient search capability. To address the problem of time-consuming black-box constraint processing for optimization problems, a multi-constraint adaptive sampling method based on improved feasible rules was proposed, an elite archive-driven inexact search method and an ε-constraint-preserving pseudo-feasible domain construction method were established, and the algorithm′s ability to explore the boundaries of the feasible domain was enhanced by dynamically scaling the feasible domain to accept high-quality nonfeasible samples during the iterative process, which improved the surrogate model-based optimization search ability. Simulation results of the Congress on Evolutionary Computation constraint optimization standard function indicate that the ε-constraint maintenance optimization method is effective in solving the multi-constraint surrogate model optimization problem compared with the existing methods. The results for the solid rocket motor rear wing pillar charge design show that the algorithm has the potential to be applied to complex engineering problems.
在以飞行器设计为代表的复杂工程设计问题中，为提高设计性能、降低研制成本、缩短设计周期，通常需要对产品性能、成本等设计指标进行优化设计[1]。同时为了提高仿真计算结果的可信度，高耗时仿真模型在实际工程设计中的应用越来越广泛，而模型仿真精度的提升会导致仿真时长指数增长，致使计算成本上升的同时优化效率低下。为解决上述问题，基于代理模型的优化方法被广泛应用于飞行器优化设计过程以达到降低设计、分析和计算消耗的目标[2]。在基于代理模型的设计方法中，自适应采样策略根据现有代理模型和样本点分布确定最佳采样点来进行高耗时模型仿真，从而获得更有利于优化集成的样本点。上述采样过程需要按照采样准则搜索得到新样本点输出并扩充具有高精度特性的样本集。根据采样准则优化指标的不同，采样准则分为局部开发、全局探索和开发/探索平衡三种策略[3]。对于黑箱约束优化问题，现有算法无法直接应用，需要与约束处理方法相结合方可得到满足约束的最优可行解。同时约束优化问题中存在的可行域边界也进一步增加了算法的困难与挑战。
在实际工程问题中，除目标函数之外，大多数约束也同样具有类似的性质。在引入约束条件之后，约束类型、不同约束个数、可行域大小、有效约束个数等复杂因素都会导致优化算法效率降低、全局最优解搜索困难等问题[4]。为此，近年来研究者们提出了针对近似优化方法的约束处理方法，使算法能够探索到空间可行域。Huang等提出的可行性规则法[5]易于实现且不需要烦琐的调参过程，在差分进化（differential evolution，DE）算法和粒子群优化（particle swarm optimization，PSO）算法等进化算法中广泛应用。但其对约束的处理过程太过严格，可能导致算法无法考虑约束边界的高质量不可行解，从而降低算法的寻优能力。
针对约束处理中存在的问题，在进行优化过程中的自适应采样时需要选取合适的约束处理方法才可能选取到高质量可行解。在基于可行性规则法处理约束的基础上，首先引入ε-约束保持法[6]，通过先接受不可行解用以探索可行域边界，再不断放缩伪可行域直至可行域，确保优化结果在约束内。基于此，本文将可行性规则法进行改进，增强其探索约束边界最优的能力，并提出约束动态松弛的近似优化（constraints dynamic relaxation approximate optimization，CDRAO）方法，最后通过数学算例和工程算例对提出的方法进行验证。
1 代理模型多约束优化方法基本框架
工程问题中的黑箱模型输出结果需要通过大量计算得出，为了提高计算效率，使用代理模型代替高耗时黑箱模型进行优化设计。基于代理模型的优化方法主要包括实验设计、代理模型构建、自适应采样三个方面，其基本流程如图1所示。
图1基于代理模型的优化方法基本流程
Fig.1Basic process of surrogate model-based model optimization method
实验设计能够为代理模型提供在设计空间上均匀分布的训练样本，通常采用拉丁超立方实验设计方法获取[7]。近似建模通过训练样本输入输出关系拟合出原始模型的代理模型，使用径向基函数近似建模方法[8]构建代理模型。自适应采样是通过现有样本点的分布以及代理模型的构建去搜索选择下一个样本点，以期在拓展样本点集中选出最优解。
Mckay提出拉丁超立方实验设计方法[9]（Latin hypercube design，LHD），能够使样本点较为均匀分布在设计空间，目前已经成为主流的实验设计方法。针对n个m维样本点的设计结果是一个矩阵，通过采用拉丁超立方实验设计方法[10]对此问题设计求解，可以在降低仿真实验次数的基础上快速得出满足均匀性要求的实验设计排列。其中样本与空间中对应位置如图2所示。
构建代理模型的方法有径向基函数（radial basis function，RBF）建模方法[6]、Kriging建模方法[11]、多项式建模方法[5]、高斯过程回归（Gaussian process regression，GPR）建模方法[12]、深度学习（deep learning，DL）建模方法[13]等。RBF模型相比其他建模方法具有设计简单、计算量小、泛化能力强等优点[14]，因此选择RBF建模方法构建代理模型。
图2拉丁超立方实验设计方法
Fig.2Latin hypercube experimental design methodology
标准RBF模型基本形式为
(1)
式中，n为训练集大小，x为待预测的点，ωi为基函数系数，φi（r）为高斯基函数。其中高斯基函数可写为
(2)
式中，ci为形状参数。
自适应采样需要基于当前代理模型和样本点分布，构造准则并对下一个样本点进行搜索选择，再通过代理模型计算新样本点的输出，实现样本集的扩充。自适应采样主要分为三种：局部开发、全局探索以及开发/探索平衡。局部开发准则基本思想是在代理模型精确逼近真实模型的假设下，寻找出代理模型的最优点，代入真实模型中验证并用来更新样本数据集。全局探索准则侧重于模型的全局最优，通过在近似模型的样本点稀疏区域选取采样点，从而更易于求出模型的全局最优解。局部开发准则易陷入局部最优，全局探索准则由于需要探索设计空间全局而易造成计算资源的浪费，因此提出非精确搜索的采样准则，通过选取既能满足结果更优又能保证选取点与已选点之间的距离较远的采样点，实现更高效率的寻优采样。
但在实际工程问题中，目标函数通常有约束限制从而成为约束优化问题，且约束多为黑箱约束。使用无约束代理模型优化算法优化约束问题时，由于约束类型、可行域大小、约束个数等多方面因素的影响，难以完成目标函数寻优。如图3所示，考虑可行域内空心采样点建立的代理模型如图中虚线所示，由此得出的最优解并不是实际的最优解。因此应开展多约束条件下的采样方法研究。
图3带约束优化存在的问题
Fig.3Problems with constraint optimization
2 多约束自适应采样方法
2.1 基于约束冲突函数的约束处理方法
单目标非线性约束优化问题可按如下形式定义：
(3)
式中：X为设计变量，通常为d维向量；f（X）为目标函数（标量函数）；XL、XU分别为设计变量的上、下限；gj（X）≤0（j=1，2，···，p）为不等式约束； hk（X）=0（k=p+1，p+2，···，m）为非设计变量上、下限的不等式约束。
在工程设计问题中多数约束计算同目标函数具有相似的复杂特征，难以找寻可行解，因此需要采用约束处理方法对约束进行处理。本文引入约束冲突函数，用以表征样本点与可行域之间的距离，使算法能够平衡约束与目标函数的关系并提升算法对可行域边界的探索能力，从而能够更好地找寻全局最优可行解。将等式约束做如下处理：
(4)
式中，hk（X）表示第k个等式约束，δ为一小量。
应用式（4），将等式约束转化为不等式约束，定义个体X违反第l个约束条件的程度为：
(5)
定义个体X违反所有约束的程度为：
(6)
由于不同约束的尺度特征不同，需进行归一化处理，以避免个别约束对G（X）的值起决定性作用：
(7)
式中： Gnor（Xi）为粒子Xi的标准化约束冲突值；Gmaxl为群体针对单个约束的最大冲突值，即
(8)
N为精英档案中样本点的个数。若Gnor（Xi）≤0，则表明个体Xi是可行个体；否则Xi为不可行个体，且Gnor（Xi）值越大说明个体Xi在群体中相对约束冲突越大。
2.2 基于ε-约束保持改进的可行性规则法
可行性规则法容易实现且优化过程中不用考虑相关参数调整，但是忽略了位于可行域边界的高质量不可行解。特别是对于代理模型优化问题，在优化初期探索阶段难以定位最优解区域，而在可行域边界外附近的高质量不可行解能够提升算法对可行域边界的搜索能力，进而提升算法的探索能力。因此提出ε-约束保持的思想[6]，在优化初期对可行域进行放大，随着算法的迭代，可行域也不断缩小，直至与真实可行域重合为止。在优化前期算法对可行域进行探索时将高质量不可行解加入样本点集，提升了算法对可行域边界的探索能力。如图4所示，在原有采样点的基础上将伪可行域内采样点考虑进去，使得代理模型更加逼近于目标函数，得出更为真实的可行最优解。算法的基本流程基本与可行性规则法相同，但将可行域改成包含可行域边界外的可放缩伪可行域。改进后三条规则如下：
1）伪可行解优于不可行解；
2）同为伪可行解，目标函数性能更理想的更优；
3）同为不可行解，约束冲突值小的更优。
图4考虑伪可行域内采样点后优化
Fig.4Optimization after considering sampling points in the pseudo-feasible domain
引入松弛参数ε，将判断可行标准由Gnor（Xi）≤0改进为Gnor（Xi）≤ε，通过调整值ε放缩伪可行域。由于Gnor（X）经过归一化处理，其取值在[0，1]之间，因此针对不同约束优化问题，松弛参数ε可以采用相同调整规律。为了确保最优点处于可行域内，ε应随迭代次数增加而逐渐减小为0，以促使伪可行域边界逐渐向可行域边界靠拢并最终重合。
2.3 基于改进可行性规则的多约束自适应采样方法
引入精英档案的思想来表征采样点性能的提升，并将精英档案初始化为实验设计样本点，经过后续迭代优化保存性能最优的若干样本点且保持数量不变。在每次迭代过程中如果有新样本点优于精英档案内存储点，则将其进行替换从而保持精英档案的最优。精英档案更新流程及其整体算法流程如图5和图6所示。结合基于ε-约束保持法改进的可行性规则法，提出改进可行性规则的多约束自适应采样方法，按照“先探索可行域，再寻优”的思想，采样过程分两阶段进行：
1）若精英档案中存在不可行解，算法应在探索到可行点后停止搜索，并在当前的种群中找到可行且密度最小的点。
2）若精英档案中全部是伪可行解，算法应在探索到优于精英档案最优的点后停止搜索，并在当前采样点中找到优于精英档案最优解且密度最小的点。
图5更新精英档案示意图
Fig.5Update elite profile diagram
图6非精确搜索采样方法流程示意图
Fig.6Non-exact search sampling method flow diagram
满足Gnor（X）≤ε的X称为伪可行解，其中ε由下式得出：
(9)
式中，cj为随算法迭代次数的增加从1减小至0的系数，定义
(10)
其中，xN代表当前精英档案样本集，s表示样本集中距离不同点的对数，di表示不同样本间的欧几里得距离。通过对样本进行归一化处理确保各个变量影响权重一致。cj=1表明精英档案为均匀分布的初始样本；cj≈0表明精英档案近似收敛到同一点上，此时算法终止。定义收敛精度e，若cj≤e，则算法收敛。
优化初始阶段设Gnor（X）≤0.1×1，即所有的约束边界被放宽10%的裕度，低于这个裕度的解可作为可行解进行阶段1采样。优化后期要求所有约束满足Gnor（X）≤0，算法只能开发可行域直到定位全局最优。值得注意的是，阶段1的采样和阶段2的采样并不具有单向性，当前采样点被判断为伪可行解之后，随着约束松弛因子ε减小，可行域不断收敛，原先的伪可行解变为不可行解，此时需要重新转到阶段1进行采样以保证种群的可行性。
由于采样过程根据采样点特征分成两种策略，精英档案的更新策略同样要适应采样策略。精英档案更新策略按如下定义：
1）若精英档案中存在不可行解，精英档案中约束违反程度最大的样本点被更新为阶段1采样获得的样本点。
2）若精英档案中全为伪可行解，精英档案中目标函数最差的解被更新为阶段2采样获得的样本点。
采样过程如图7所示。
图7基于改进可行性规则的多约束自适应采样方法
Fig.7Multi-constraint adaptive sampling method based on improved feasible rules
外部代理模型优化的最优解追踪调整为：优先追踪精英档案中满足约束且目标函数值最小的点，若精英档案中不存在满足约束的点，最优解设置为当前采样点中约束冲突最小点。若采样点都不可行，最优解可以帮助其逐渐向可行域收敛。采样点都可行时，最优解要保证其在目标函数上的优越性。面向多约束问题的非精确采样代理模型优化方法流程如图8所示，其中流程中的代理模型构建采用基于样本局部密度的径向基建模方法，搜索采样过程采用改进粒子群算法。
图8多约束问题的非精确采样代理模型优化流程
Fig.8Non-exact sampling surrogate model optimization process for multi-constraint problems
3 方法验证
3.1 数值算例
为了验证约束动态松弛优化方法的有效性，从进化计算会议（IEEE Congress on Evolutionary Computation，CEC）约束优化标准函数中选取8个约束优化进行测试，其目标函数、约束条件、自变量范围以及函数全局极小值如表1所示。CEC标准测试问题是CEC组委会面向不同问题类型发布的标准测试算例，这些算例为全局最优解已知的函数，被广泛用于优化方法性能验证[17-18]。
为充分测试方法的适用性和高效性，选取的测试问题设计变量维度涵盖2维到13维，约束形式包括不等式与等式约束。在处理等式约束时，由于难以找到严格满足等式约束的解，故将其做如下处理：
(11)
式中，δ为一小量，通常取0.001。
在进行测试算例时，针对每一个算例均独立运行25次以上，确保运行结果的稳定。将优化结果与近年来发表的约束优化混沌灰狼算法[19]（chaotic grey wolf optimization，CGWO）、扩展平衡排序方法[20]（the extended balanced ranking method，E-BRM）和全局人工蜂群算法[21]（modified global best artificial bee colony，MGABC）进行对比，得出结果如表2所示。
表1约束优化标准函数
Tab.1Constraint optimization criteria functions
根据表2中结果，将不同算法对每一个标准函数求出的最优可行解加粗，CDRAO算法在八个测试函数上的优化结果均为最优的效果。而CGWO、E-BRM和MGABC都无法完全寻优出文中测试算例的最优，且CGWO、E-BRM无法在部分算例上找出可行域，导致寻优结果无效，MGABC算法不能够优化出g05函数和g07函数的全局最优。
CDRAO算法针对上述测试函数均能得到约束全局最优，且结果稳定性能够满足要求。针对每个函数迭代运算25次，计算出迭代收敛时迭代次数的平均值，每次计算设置的样本点个数保持不变，采用CDRAO算法对8个测试函数优化的迭代次数如表3所示。
表2约束采样方法算例测试结果
Tab.2Constraint sampling method study test results
表3CDRAO算法迭代次数统计
Tab.3Statistics on the number of iterations of CDRAO algorithm
相比于文献中其他算法的收敛速率，CDRAO算法收敛速率在200次以内，且低维函数和高维函数次数相差不大，适用性较强，表明提出的CDRAO算法是性能优越的约束优化算法，拥有进行复杂工程约束优化应用的潜力。
3.2 固体发动机后翼柱装药设计算例
采用CDRAO方法对固体火箭发动机装药设计进行优化，优化目标为药柱燃烧过程燃面稳定性，约束条件为发动机药柱质量，以验证算法的有效性。
3.2.1 计算模型
仿真模型采用固体火箭发动机性能匹配框架[22]实现。药柱形状选用翼柱形，主要考虑后翼柱装药，图9展示了算例的装药构型结构。优化设计变量选择后翼柱的5个几何结构参数，其中包含4个连续变量和1个离散变量，其类型及变量范围如表4所示。除此之外，其他常数参数设置见表4。装药几何结构可由表中设计变量唯一确定。
图9后翼柱装药结构示意图
Fig.9Schematic diagram of the rear wing pillar loading structure
表4后翼柱型装药设计变量
Tab.4Rear wing column type charge design variables
为了量化燃面设计的稳定性，定义算例优化问题为：
(12)
式中，X表示装药几何构型设计变量，A（X）表示根据当前构型计算出燃面燃烧厚度变化曲线，fSTD表示计算燃面曲线的标准差，m为装药质量，m0为初始装药质量最小值。本算例中m0取240 kg。
3.2.2 计算结果
初始算法种群设置为10，优化过程中目标函数的收敛情况和优化结果的燃面曲线如图10、图11所示。
图10fSTD历史收敛曲线
Fig.10fSTD historical convergence curve
图11后翼柱装药设计燃面曲线
Fig.11Rear wing pillar charge design combustion surface curve diagram
最终优化设计结果如表5所示。
表5后翼柱型装药优化结果
Tab.5Optimized results for rear wing pillar type loading
由图10、图11得出，CDRAO方法能够快速设计出燃面稳定的后翼柱型装药结构，迭代次数仅需120多次，展现出其高效性。而优化设计结果中装药质量为240.24 kg，满足模型中关于约束条件的限制。证明CDRAO方法对约束处理的可靠性，能够应用于实际工程问题中。
4 结论
针对高耗时模型优化设计问题，对代理模型自适应采样阶段开展研究。考虑目标函数存在的各种类型约束，提出了多约束自适应采样方法，通过数学算例和工程算例对文中提出方法展开验证。
1）在基于约束冲突函数处理约束问题的基础上提出基于ε-约束保持改进的可行性规则法，制定“先探索可行域，再寻优”的二阶段采样策略，增强算法探索可行域附近高质量不可行解的能力，提升算法对约束边界的探索精度。
2）提出了基于改进可行性规则的多约束自适应采样方法，通过引入精英档案的思想迭代更新可行域的高质量解，通过缩小ε的值不断使当前域逼近于可行域，避免了优化过程中采样可行域边界精度不够的问题，从而提高了优化结果的准确性。
3）采用CEC优化约束标准函数对CDRAO方法进行测试，并与其他算法进行比对，得出算法能够较好地得出约束函数的最优值且具有较高的稳定性。将CDRAO方法应用到固体火箭发动机后翼柱装药设计算例中，得出满足设计要求的优化设计结果。表明文章提出的方法可以应用于飞行器设计方面的优化问题。
