考虑高耗时约束的追峰采样智能探索方法

龙腾，毛能峰，史人赫，武宇飞，沈敦亮

1. 北京理工大学 宇航学院，北京 100081 2.北京理工大学 飞行器动力学与控制教育部重点实验室，北京 100081 3.清华大学 航天航空学院，北京 100084 4.北京宇航系统工程研究所，北京 100076

数值计算技术、学科建模技术与计算机软硬件的发展使得高精度分析模型在飞行器系统设计的工程实践中应用日益广泛，例如，采用计算流体力学(Computational Fluid Dynamics, CFD)模型、计算电磁学(Computational Electromagnetics, CEM)模型、结构有限元分析(Finite Element Analysis, FEA)模型等高精度分析模型，可以改善设计结果的可信度，但存在计算耗时长的问题。基于梯度的局部搜索算法(如增广拉格朗日乘子法、序列二次规划法和可行方向法)通常利用有限差分法计算梯度信息，需要反复调用分析模型；基于概率的全局探索算法(如遗传算法、粒子群算法和模拟退火法)同样需要成千上万次的调用分析模型[1]。为了缓解上述问题，构造高精度分析模型的代理模型，可以有效降低性能指标预测成本。神经网络是一种模拟生物神经系统结构和功能的代理模型，理论上可以逼近任意非线性函数，因此得到了广泛关注。张嘉良等[2]将神经网络与遗传算法相结合，用于求解特定飞行条件下的巡航飞行器增程优化设计问题。雷玉昌等[3]利用神经网络建立双后掠翼飞行器外形参数与性能参数之间的映射关系，用于求解气动布局多目标优化设计问题。孙伟和张呈林[4]建立直升机桨叶气动性能的神经网络预测模型，用于求解气动外形多目标优化设计问题。

基于代理模型的近似优化策略已成为飞行器设计优化领域的研究热点，以降低计算成本，缩短设计周期[1]。近似优化策略主要分为静态近似优化策略和自适应近似优化策略两类[5]。与静态近似优化策略相比，自适应近似优化策略采用动态更新代理模型的方式改善可能存在全局最优解区域的局部近似精度，具有更高的效率和更好的全局收敛性，其代理模型管理与更新策略主要包括基于空间缩减的序列采样方法和基于空间填充的序列采样方法[1]。基于空间缩减的序列采样方法的典型算法包括信赖域方法[6]、分割平面方法[7]以及重点设计空间方法[8]等。基于空间填充的序列采样方法的典型算法包括高效全局优化策略[9]、追峰采样方法[10]、基于进化操作的采样策略[11]以及基于模糊聚类的采样策略[12]等。

标准追峰采样方法(Mode-Pursuing Sampling method, MPS)适用于求解复杂黑箱函数的全局最优解，优化效率高、鲁棒性强且支持并行计算。但是Duan等[13]指出标准MPS方法求解高维优化问题效率不高，对此Cheng等[14]提出了基于信赖域的追峰采样方法(Trust Region based Mode Pursuing Sampling, TRMPS)。针对约束问题，Kazemi等[15]引入约束松弛条件，提出了一种约束追峰采样优化方法(Constraint-importance Mode Pursuing Sampling, CiMPS)。此外粟华等[16]提出了一种改进的动态追峰采样方法(Dynamic Mode Pursuing Sampling, DMPS)，缓解算法过早收敛、陷入局部最优解的问题。Wu等[17]提出了一种应用有偏坐标扰动的追峰采样方法(Mode Pursuing Sampling method using Discriminative Coordinate Perturbation, MPS-DCP)，进一步改善高维问题的优化效率，但其复杂约束处理能力有待提升。然而实际工程优化问题中，约束计算成本高的现象普遍存在[15]。CiMPS的约束处理机制忽略约束函数计算耗时，其约束函数调用次数通常比目标函数大3个数量级[18]。上述问题限制了追峰采样方法在飞行器设计优化领域的工程应用。

处理约束优化问题的传统方法主要分为基于罚函数的方法和不依赖罚函数的方法。罚函数法[19]通过构造增广目标函数将原问题转变为无约束优化问题进行求解，存在计算量大、收敛速度慢与罚系数选取困难等缺陷。对于不依赖罚函数的约束处理方法，Crossley[20]提出了结合KS方程(Kreisselmeier-Steinhauser function)的遗传算法，用于求解旋翼系统的多目标约束优化问题。史人赫等[21]提出了基于过滤器的序列径向基函数优化方法，综合考虑目标函数值与最大约束违背度构建过滤器，并通过支持向量机辨识可行性区域指导有偏采样，能够有效求解复杂约束优化问题。本文将不依赖罚函数的约束处理方法引入MPS框架，提高其求解高耗时约束优化问题的效率。

针对涉及高耗时黑箱分析模型的飞行器约束优化问题，本文提出了一种基于过滤器的MPS-DCP设计空间智能探索方法(Filter-based Mode Pursuing Sampling intelligent exploring method using Discriminative Coordinate Perturbation, FMPS-DCP)。通过训练径向基函数网络(Radial Basis Function Network, RBFN)预测目标函数与约束条件响应值，并应用KS方程构造约束违背度函数，综合考虑目标函数值与约束违背度函数值，基于Pareto支配概念构造或更新过滤器，筛选优质简单样本点，再根据预测值准则与距离准则从过滤器接受样本点中选择出新增样本点，引导优化过程快速向全局可行最优解收敛。最后通过标准约束测试算例和工程优化案例，验证FMPS-DCP方法的性能。

1 无约束追峰采样方法简介

1.1 径向基函数网络简介

RBFN是一种高效的单隐层前馈神经网络，拓扑结构简单、理论上可以逼近任意非线性函数。其隐层神经元的激活函数选用径向基函数，输出层为隐层输出的线性组合。RBFN的回归预测模型如式(1)所示[22]：

(1)

式中：q为隐层神经元数量；wi为第i个隐层神经元的连接权重；φ(x,xi,c)为径向基函数；xi为对应隐层神经元的中心；c为径向基函数的形状参数。

RBFN训练时，首先确定神经元中心xi，再确定权重wi和形状参数c。由于实际复杂飞行器系统设计中训练样本规模有限，本文选取优化过程中当前所有样本点作为隐层神经元中心，并按式(2)计算连接权重wi[1]。

w=A-1F

A=

(2)

式中：w=[w1,w2,…,wn]T为连接权重系数向量；F=[y1,y2,…,yn]T为训练样本处的真实模型响应；A为径向基函数矩阵；n为训练样本规模。

1.2 标准追峰采样方法简介

标准MPS方法[10]初始样本点数计算式为

(3)

式中：nv为设计变量个数。计算目标函数值加入高精度样本点集Y并构造RBFN。

虽然标准MPS求解低维问题具有显著优势，但是求解高维问题的效率较低[13]，而且缺乏处理高耗时约束的能力。

1.3 基于有偏坐标扰动的追峰采样方法简介

MPS-DCP[17]将有偏坐标扰动的思想与标准MPS相结合，以缓解其高维问题求解效率低的缺陷。MPS-DCP初始采样过程与标准MPS方法相同。迭代过程中选取Y内目标函数值最小的10·nv个高精度样本点，构造目标函数的RBFN。

综合考虑各设计变量的灵敏度指标[17]与优化改善情况构造灵敏度影响因子，并结合当前迭代次数计算设计变量的扰动概率p。由此对当前最优点施加有偏坐标扰动，得到nc个简单样本点，通过定制的预测值准则与距离准则，选择ns个RBFN预测值较小且距离现有高精度样本点较远的优质简单样本点作为新增样本点。此外，根据优化改善情况动态调整有偏坐标扰动步长σn。若当前最优解与上一步迭代结果相同，则缩减步长，加强局部搜索能力；若当前最优解小于上一步迭代结果，则扩大步长，保证算法的全局探索性。

与TRMPS等改进算法相比，MPS-DCP求解高维问题的收敛性能与优化效率有明显优势[17]。然而MPS-DCP中构造约束函数RBFN的处理方式难以有效求解包含复杂约束的实际工程优化问题。因此，亟待将高效的约束处理机制引入MPS框架，以进一步提升MPS算法的工程实用性。

2 约束追峰采样智能探索方法

2.1 约束处理方法

本文采用不依赖罚函数的约束处理方法，具体内容叙述如下。

2.1.1 KS方程法

KS方程是在指数空间上约束函数集合的可微最大包络。本文采用如式(4)所示的KS方程修改形式[23]，以缓解数值奇异的问题。

(4)

以G24问题(定义见3.2.1节)为例，说明KS方程缩聚约束后可行域变化情况，原可行域与KS方程包络可行域如图1所示。其中一个主动约束处(图中黑框部分)放大图如图2所示。

由图1和图2可以看出，对于G24问题，KS方程包络可行域与原可行域大致相同，但在主动约束处，KS方程包络可行域明显更小，且ρ越大越接近原可行域。因此，KS方程存在可行域偏保守的问题，但是ρ取得过大仍会导致数值奇异。为解决上述问题，本文在优化过程中采用自适应方法[23]动态更新ρ，计算公式为

图1 G24问题可行域Fig.1 Feasible region of G24

图2 可行域放大图Fig.2 Enlarged view of feasible region

lg(ρn+1)=

(5)

式中：ρ1=ρn+1×10-3；KS′为ρ的影响因子，计算式为

(6)

2.1.2 过滤器方法

考虑到实际工程问题中仿真分析结果可能存在误差，允许样本点在一定范围内违背约束，从而更充分地探索设计空间，加强算法的全局收敛性。由此引入约束容差KSmax，将满足KS(ρn,x)

定义1：支配关系

对于样本点x(i)与x(l)，当且仅当满足式(7)时，称点x(i)支配x(l)；否则，两者互不支配。

f(x(i))≤f(x(l))∩KS(ρn,x(i))≤KS(ρn,x(l))

(7)

定义2：过滤器

过滤器是互不支配样本点的集合。过滤器的示意图如图3所示。当且仅当样本点位于图中过滤器更新或扩充区域，样本点被过滤器接受。否则，样本点被过滤器拒绝。

图3 过滤器示意图Fig.3 Schematic of filter

2.2 算法流程

假设高精度样本点的目标函数与约束函数响应值可通过调用一次分析模型同时获取，而简单样本点通过RBFN预测目标函数与约束函数值。通过模型调用次数nse衡量优化问题的计算成本。此外，认为可行点优于不可行点。对于可行点，目标函数值较小的更优；对于不可行点，约束函数值较小的更优。FMPS-DCP的算法流程如图4所示，具体步骤如下:

图4 FMPS-DCP算法流程图Fig.4 Flowchart of FMPS-DCP

步骤1采用标准拉丁超方试验设计方法选取nv+1个初始高精度样本点。判断是否存在可行点，如果不存在可行样本点，执行步骤2。如果存在，则继续采样至样本点数达到n0，计算公式如下:

(8)

计算目标函数与约束函数真实响应值并加入高精度样本点集Y，执行步骤3。

(9)

式中：xk为Y内高精度样本点；Tcoincide为距离容忍度。具体流程见2.3节。

步骤3利用Y内所有样本点训练目标函数与各约束函数的RBFN。

步骤4如果为第1次迭代，确定Y内样本点之间的支配关系，构造过滤器。否则，根据上一步迭代新增样本点与过滤器内所有样本点之间的支配关系，更新过滤器。

步骤5扰动概率p计算方式与MPS-DCP方法类似，但由于存在约束函数，则分别计算目标函数与约束函数的灵敏度指标[17]sF与sKS，构造统一灵敏度指标s=sF+sKS，并根据式(10)归一化：

(10)

步骤6首先利用得到的过滤器筛选上述简单样本点，再根据预测值准则和距离准则评价备选点的得分，选取评分最好的ns个样本点。新增样本点选择策略具体内容见2.4节。计算目标函数与约束函数真实响应值并加入高精度样本点集Y。最后与MPS-DCP方法类似，比较新增样本点与当前最优点的最优性与可行性，更新有偏坐标扰动步长。

步骤8利用序列二次规划方法对PRSM模型进行局部优化，计算最优解的目标函数与约束函数真实响应值并加入高精度样本点集Y。

步骤9如果Y内样本点数达到了给定的最大模型调用次数，则优化终止，输出当前最优解。否则，更新灵敏度指标，返回步骤3。

2.3 初始采样子优化问题

为改善优化结果的可行性，需要保证初始采样存在可行样本点，子优化问题求解的具体步骤如下:

步骤3计算样本点与Y内所有样本点之间的最小欧氏距离，构造初始采样子优化问题的约束函数。

步骤4利用序列二次规划法求解初始采样子优化问题，计算最优解的目标函数与约束函数真实响应值并加入高精度样本点集Y，更新模型调用次数。判断Y内是否存在可行点或达到最大模型调用次数，若满足，执行步骤5。否则，返回步骤1。

步骤5如果Y内样本点数小于最大模型调用次数，执行步骤6。否则，返回集合Y和KS方程控制参数ρ，子优化终止。

步骤6利用拉丁超方试验设计方法生成n0-nsub个样本点，其中nsub为Y内样本点数。计算目标函数与约束函数真实响应值并加入高精度样本点集Y，子优化终止。

初始采样子优化问题求解过程伪代码如算法1所示。

算法1 初始采样子优化问题求解方法

2.4 基于过滤器的新增样本点选择策略

为高效处理高耗时约束，首先基于过滤器筛选优质简单样本点，再根据定制的预测值准则与距离准则选择新增样本点，引导采样快速趋向全局可行最优解，具体步骤如下:

步骤2如果过滤器接受点数小于拟新增样本点数ns，则偏向可行性，取约束违背度函数预测值为预测值准则评价指标。否则，偏向最优性，取目标函数预测值为评价指标。

步骤3计算样本点与Y内所有样本点之间欧氏距离的最小值，作为距离准则的评价指标。

步骤4根据当前迭代次数niter、优化持续未改善次数Cstall与ns确定权重系数集合ws。

步骤5如果已新增样本点数nadd小于ns，则根据nadd确定权重系数w，执行步骤6。否则，执行步骤7。

步骤6按式(11)计算总评分。

(11)

如果样本点与Y内所有样本点之间距离的最小值小于阈值Tcoincide，则去除该样本点。选择总评分最小的样本点为新增样本点，并加入Y与新增样本点集Ya，更新nadd。如果nadd小于ns，返回步骤5。否则，执行步骤7。

步骤7样本点选择结束，输出新增样本点集Ya与更新后的高精度样本点集Y。

基于过滤器的样本点选择策略伪代码如算法2所示。

算法2 基于过滤器的样本点选择策略Algorithm 2 Filter based sample point selection strategy

2.5 算法收敛性分析

对FMPS-DCP方法的收敛性开展理论分析并加以证明。根据文献[17]，定义如下引理。

条件a，算法采用精英策略。

(12)

条件b，对于非全局最优点的任意样本点x′，对应水平集L(x′)={x|f(x)

(13)

3 算法性能测试

为验证FMPS-DCP方法的性能，采用7个中低维、4个高维标准约束测试算例和1个实际工程约束优化案例进行测试，将测试算例中的目标函数和约束函数视作高耗时黑箱模型。对于中低维问题，本文认为达到最优解目标值即得到了最优解，算法终止，通过统计不同算法在最大的模型调用次数限制下达到给定最优解目标值所需的模型调用次数，并与CiMPS、MPS-DCP和改进随机径向基函数方法(Extended ConstrLMSRBF)进行对比分析，评价算法的性能。由于CiMPS难以求解涉及高耗时约束的优化问题，Extended ConstrLMSRBF方法未提供源代码，对于高维问题，本文与约束差分进化算法(Kriging assisted Constrained Differential Evolution, KRG-CDE)[24]进行对比；对于实际工程约束优化案例，本文与改进自适应响应面优化策略(Adaptive Response Surface Method using Intelligent Space Exploration Strategy, ARSM-ISES)进行对比。为降低随机因素的影响，对各数值算例，不同算法分别连续运行30次。各算法在模型调用次数大于最大模型调用次数时，算法终止。

3.1 参数设置

对于中低维问题，MPS-DCP、FMPS-DCP和Extend ConstrLMSRBF目标函数与约束函数调用次数相等，最大调用次数设为500。CiMPS约束函数调用次数通常多于目标函数，最大目标函数与约束函数调用次数分别设为500与10 000。对于30维C0130和C0230问题，KRG-CDE和FMPS-DCP最大的目标函数与约束函数调用次数均设为600；对于50维C0150和C0250问题，各算法最大的目标函数与约束函数调用次数均设为1 500。FMPS-DCP约束容差设为0.05。

3.2 问题描述

3.2.1 标准约束测试算例

选用的中低维标准约束测试算例来自文献[25-29]，高维标准约束测试算例来自文献[30]，基本信息如表1和表2所示。

表1 中低维标准约束测试算例

表2 高维标准约束测试算例

3.2.2 工程案例

本文选用文献[31]中的全电推进卫星平台多学科设计优化问题(如图5所示)，验证FMPS-DCP方法的工程实用性。全电推进卫星平台多学科设计优化案例在满足轨道转移时间、轨道位置保持精度、整星功率需求以及基频等约束条件下，以整星质量最小为优化目标，涉及轨道转移、位置保持、供配电、热控、姿控和结构6个学科。为了实现学科解耦，本文采用定点迭代法实现多学科分析(Multidisciplinary Analysis, MDA)[31]，获取满足系统一致性设计约束的样本点。在Intel i5 6300U 2.40 GHz CPU和8 GB内存的电脑上，利用定点迭代法进行1次多学科分析平均所需时间约10 min。

图5 全电推进卫星示意图[31]Fig.5 Schematic of all-electric satellite[31]

全电推进卫星平台多学科设计优化问题的优化数学模型如式(14)所示。

findX=[α,β,φ,dT,dN,Asa,Cs,Ar,Hw,

SH,CH,TBH,SP,CSP,TBP]T

(14)

式中：α、β、φ分别为第一阶段偏航、俯仰角和第二阶段俯仰角；dT和dN分别为推力器T和N方向安装位置；Asa为太阳能帆板面积；Cs为电池容量；Ar为散热器面积；Hw为动量轮角动量；SH和CH分别为服务舱和通信舱芯子厚度；TBH为中心承力筒芯子厚度；SP和CSP分别为服务舱和通信舱铺层厚度；TBP为中心承力筒铺层厚度。Msatellite为全电推进卫星系统总质量；mi为第i个学科/分系统的质量预算；mfuel为总推进剂质量预算；mothers=746.2 kg为系统固定质量；tf为总轨道转移时间；λmax和imax分别为东西和南北位置保持精度；PBOL和PEOL分别为寿命初、末期功率；DOD为放电深度；T0为稳态温度；cAC为动量轮容量余量；fX和fY分别为X、Y方向一阶弯曲频率。

3.3 优化结果

3.3.1 标准约束测试算例优化结果

1) 优化结果对比分析

对于中低维标准约束测试算例，3种算法优化结果如表3所示，Extend ConstrLMSRBF优化结果引自文献[32]；由于MPS-DCP无法在给定的最大模型调用次数限制下获取PVD4、SR7、G7和G9等问题可行解，本研究仅针对G1和G24问题，对比了MPS-DCP和FMPS-DCP优化过程迭代曲线，如图6所示，图中与X轴平行虚线代表最优解目标值。对于高维标准约束测试算例，优化结果如表4所示。上述表中加粗部分代表各算例的最优结果。

表3 中低维标准约束测试算例优化结果Table 3 Optimization results of middle and low dimension constrained numerical benchmark problems

图6 中低维算例优化过程迭代曲线Fig.6 Convergence history curves for middle and low dimension constrained numerical examples

表4 高维标准约束测试算例优化结果Table 4 Optimization results of high dimension constrained numerical benchmark problems

对于中低维测试算例，在给定的模型调用次数限制下，FMPS-DCP都存在达到最优解目标值的优化结果；而CiMPS只有算例G24可以达到最优解目标值，Extended ConstrLMSRBF除SR7与G1外，其余算例都存在达到最优解目标值的优化结果。FMPS-DCP与CiMPS和Extended ConstrLMSRBF相比，收敛性更好，且远优于CiMPS。CiMPS通常在目标函数调用次数未达到最大次数前，约束函数调用次数已达到最大值，无法收敛到最优解，表明CiMPS难以求解高耗时约束优化问题。FMPS-DCP的模型调用次数均远少于CiMPS算法，且小于Extended ConstrLMSRBF，特别是SR7、G6与G24问题，FMPS-DCP算法优化效率具有显著优势；对于PVD4问题，FMPS-DCP未达到最优解目标值的次数略多于Extended ConstrLMSRBF，但是上述差异在工程优化领域可以忽略。在鲁棒性方面，因为对于大多数问题，CiMPS均无法达到最优解目标值，本文只对比FMPS-DCP与Extended ConstrLMSRBF。对于所有的算例，FMPS-DCP算法的标准差均小于Extended ConstrLMSRBF，说明FMPS-DCP具有更好的鲁棒性。上述优化结果表明，本文提出的FMPS-DCP算法达到给定最优解目标值所需的模型调用次数更少且鲁棒性更好，从而证明了FMPS-DCP设计空间智能探索方法具有更好的优化效率与鲁棒性。

对于G1和G24问题，MPS-DCP收敛速度随优化进行明显降低，甚至无法收敛到最优解目标值；FMPS-DCP在优化初期性能优于MPS-DCP且可快速收敛到最优解目标值。上述结果表明在MPS框架中仅通过简单罚函数法难以有效处理高耗时约束，本文FMPS-DCP方法采用的约束处理策略与新增样本点选择策略可以引导优化过程快速向全局可行最优解收敛，从而验证了FMPS-DCP中高耗时约束处理策略和代理模型更新管理机制的有效性与各部分策略之间的相互适应性。

对于高维测试算例，本文提出的FMPS-DCP算法可以快速收敛到全局最优解，最优性明显优于KRG-CDE，最优解标准差更小，表明FMPS-DCP可用于求解高维问题且鲁棒性较好。

2) 调节参数影响性分析

影响算法性能的调节参数，包括序列新增样本点数ns、阈值参数Tcoincide和Tfail。对于PVD4、G1、G7和SR7问题，采用不同调节参数优化的迭代曲线如图7所示，图中与X轴平行虚线代表最优解目标值。

图7 调节参数影响性分析Fig.7 Tuning parameters influence analysis

3.3.2 工程案例优化结果

FMPS-DCP与ARSM-ISES优化结果如表5所示，其中ARSM-ISES优化结果引自文献[31]。FMPS-DCP与ARSM-ISES优化后的设计变量与约束条件如表6与表7所示。

由表6与表7可以看出，与初始设计方案相比，两种算法优化后的设计变量均有明显调整，约束条件均在约束容差范围内满足要求，表明了优化结果的可行性，其中，放电深度与X/Y方向一阶弯曲频率均成为主动约束，表明了优化结果的最优性。由表5可以看出，FMPS-DCP方法优化后的整星质量为2 357.7 kg，与初始设计方案相比减少了186.2 kg，降低了约7.32%；该结果与ARSM-ISES方法优化结果2 358.6 kg相当，但是模型调用次数从1 200次降为200次，优化时间从34.3 h降为5.3 h，计算成本降低84.5%，有助于缩短全电推进卫星平台的设计周期。上述结果验证了FMPS-DCP设计空间智能探索方法求解包含高耗时约束的全电推进卫星平台多学科设计优化问题的有效性与工程实用性。

表5 全电推进卫星平台多学科设计优化问题优化结果

表6 全电推进卫星平台多学科设计优化问题设计变量优化结果Table 6 Optimized design variables of all-electric satellite MDO problem

表7 全电推进卫星平台多学科设计优化问题优化后的约束条件Table 7 Optimized constraint values of all-electric satellite MDO problem

4 结 论

1) 针对高耗时约束优化问题求解复杂，计算成本高的难题，本文提出了一种基于过滤器的MPS-DCP设计空间智能探索方法。在通过训练RBFN预测目标函数与约束条件响应值并应用KS方程法构造约束违背度函数的基础上，定制了一套新增样本点选择策略引导优化过程快速向全局可行最优解收敛，从而提高求解高耗时约束优化问题的效率。

2) 标准约束测试算例与优化结果表明，与CiMPS、Extended ConstrLMSRBF和KRG-CDE方法相比，FMPS-DCP方法具有更好的优化效率与鲁棒性。全电推进卫星平台多学科设计优化问题求解结果表明，FMPS-DCP方法具有较好的工程实用性。

3) 未来将进一步探索模型特征驱动的优化策略研究，进一步提高飞行器系统设计的优化效率。