现今复杂产品的设计往往包含大量高维昂贵黑箱优化问题的求解，需要进行数以千计的高精度仿真分析，耗费大量计算资源^[1]。尽管目前计算机运算能力和运算速度持续增强，但数值仿真分析软件诸如有限元分析(finite element analysis, FEA)、计算流体力学(computational fluid dynamics, CFD)等同样变得越来越复杂与精细^[2]。目前最有效方法是建立代理模型来代替实际昂贵黑箱模型，从而通过对代理模型的优化得到真实模型的近似优化值。

当前广泛使用的代理模型包括多项式响应面模型(polynomial response surfaces, PRS)、径向基函数模型(radial basis functions, RBF)和克里金模型(Kriging, KRG)等。基于代理模型优化方法能够准确近似反映昂贵黑箱问题，有效降低计算成本，在工程优化领域得到深入研究与广泛应用^[3-4]。然而GOEL等^[5]指出对于不同特性优化问题，不同代理模型拥有各自优势和不足，在诸如精度、效率、鲁棒性和透明度等方面表现各异。另外，经典代理模型虽然能够较好地解决工程优化领域的低维问题，但随着优化问题维度和复杂度的增加，构建代理模型所需的计算花费也将呈指数增长，大大降低了高维问题优化效率。组合代理模型拥有选取的多个单一代理模型的优点，能够同时获取更多关于未知优化问题的信息，有效降低粗糙近似模型的消极影响，提高高维昂贵黑箱模型预测精度。Wang等^[6]提出了基于RBF和KRG组合代理模型的高维优化方法SAEAs，通过不同维度数值试验，验证了该方法在处理高维问题时的优越性。Gu等^[7]提出了基于PRS、RBF和KRG组合代理模型的高维优化方法HMBSM，在搜索过程中同时采用3种代理模型的响应值大小，优先选取对应3种近似模型目标函数值都较低的样本点，成功应用于汽车轻量化设计。

虽然使用代理模型可以一定程度上缓解高维昂贵黑箱问题的计算压力，但是维度过高仍将导致基于代理模型优化方法难以获得令人满意的优化精度和效率^[8]。设计空间缩减策略作为一种积极有效的方法被提出来克服这个难题，通过分析设计变量和目标函数之间关系，探索全局最优解可能位于的区域来提高优化质量和优化效率。Farias等^[9]提出基于信赖域流体结构耦合设计优化框架，通过自适应更新设计子空间，逐步提高近似模型在局部区域内的精度。Long等^[10]提出了智能空间探索策略，进而发展了一种改进的高效自适应响应面方法。李春娜等^[11]使用模糊聚类算法确定关键的设计子空间，然后在每个子空间内通过最大化目标函数的期望提高函数和最小化模型预测目标来增加新的样本，发展了一种适用于气动优化的高效自适应全局优化方法。

考虑到高维昂贵黑箱问题的复杂性，无论采用何种设计空间探索策略，都无法有效提高局部搜索能力。本文联合多种设计空间探索策略，提出了多层设计空间缩减策略，同时结合组合代理模型，以提高其处理高维昂贵黑箱优化问题的效率。本文提出的基于多层设计空间缩减策略的近似高维优化方法，采用多层设计空间缩减策略确定有效的设计子空间，并在缩减的设计子空间内自适应增加有效样本点。组合代理模型和缩减的设计子空间在迭代过程中自适应地更新、重建，直至收敛获得全局最优。最后通过标准优化函数和工程优化实例测试分析，验证近似高维优化方法HSRAHO的性能优势。

1 组合代理模型

组合代理模型由多个选用代理模型线性加权组合而成，其数学表达式为：

(1)

式中: 表示组合代理模型; 和w_i分别表示选用的第i个代理模型和对应的权重因子; m_s表示选用代理模型数量; x为设计变量。

在迭代过程中, 选用代理模型权重因子通过智能求解如下优化问题获得

(2)

式中: E_e和E_i分别表示构建的组合代理模型和第i个选用代理模型的预测均方根误差(prediction sum of squares of root mean square, PRESSRMS); 表示所有选用代理模型的预测均方根误差的平均值; N为训练样本数量; y(x_k)为第k个训练样本点x_k处真实响应值; 表示第i个选用代理模型(由其余N-1个训练样本点构造而成, 即留一交叉验证方法)在第k个训练样本点x_k处预测响应值。

值得注意的是, 上述权重因子计算方法中包含2个未知参数α和β, 他们分别控制平均和选用代理模型的重要程度。参数α值越小, 参数β负值越大, 则最佳代理模型分配到权重因子越大; 参数α值越大, 参数β负值越小, 则代理模型平均重要程度越高。参数α和β自由选择能够使权重因子计算方法具有高度灵活性, 合理的参数值将提高组合代理模型预测精度。本文将参数α和β作为设计变量进行优化, 权重因子w_i通过公式(2)计算得到。另外, 文中采用的约束条件, 有助于提高权重因子寻优质量和效率, 序列二次规划算法(sequential quadratic programming, SQP)被用来求解上述优化问题。本文选用3种经典代理模型PRS、RBF和KRG构建组合代理模型, 3种代理模型的详细介绍可参看文献[12-13]。

2 多层设计空间缩减策略

本文提出多层设计空间缩减策略以帮助研究人员识别全局最优解潜在区域, 减少不必要的计算资源, 提高优化效率。该策略采用已有的“昂贵”样本点(指采用真实分析模型计算得到实际目标函数值)确定2个减小的设计子空间“subspace1”和“subspace2”。文中按照所有样本点对应的真实目标函数值进行排序, 选择目标函数值较小的前M₁个样本点(M₁取值为“昂贵”样本点数量一半, 其中样本点数量根据样本点数据库动态变化)。然后使用模糊C均值聚类算法^[14](fuzzy C-means, FCM)计算得到设计子空间subspace1, 具体计算公式如下

(3)

式中: S_subspace1^topM₁表示选用的真实目标函数值较小的前M₁个“昂贵”样本点; n为设计变量; X^ranki和Y^ranki分别表示函数值较小的第i个样本点的设计变量值和对应的真实目标函数值。另外, c为聚类中心数, 本文c设为2;S_fcmi^Q_i表示在第i个聚类集内包含的所有Q_i个“昂贵”样本点, Q_i大小由模糊聚类算法确定, 即Q₁+Q₂=M₁。同时包含最小目标函数值和较小的平均目标函数值的聚类集将被确定为最优聚类集。这种情况下, 设计子空间subspace1将由最优聚类集内所有样本点在每个维度上的最大和最小坐标值计算得到。否则, 设计子空间subspace1将由所有聚类集计算得到的各个子空间合并确定, 以此增加其包含实际全局最优点的概率。

相比设计子空间subspace1, 设计子空间subspace2是一个更小的立方体空间。选用真实目标函数值较小的前M₂个“昂贵”样本点(M₂设为8), 分别计算subspace2的中心C和尺寸W。

(4)

式中, S_subspace2^topM₂表示选用的真实目标函数值较小的前M₂个“昂贵”样本点。L_j和U_j分别表示特定区域的上下边界值, 用来决定设计子空间subspace2在第j维上的尺寸大小。另外, 当前最优点将作为下一次迭代过程中设计子空间subspace2的中心。一旦设计子空间subspace2的中心和尺寸确定, subspace2可由如下公式计算得到。

(5)

式中, B_j^L, B_j^U, B_j^LS和B_j^US分别表示初始全局空间(original global space, OGS)和设计子空间subspace2在第j维上的上下边界值。尽管当前最优点可能与实际全局最优点不同, 尤其在早期迭代过程中, 但当前最优点是现有条件下包含有效信息最多的样本点。在寻优过程中, 设计子空间subspace2将随着中心和尺寸的变化而逐步向实际全局最优点潜在区域移动、收缩或扩大。值得注意的是, 当前最优点可能落在初始全局空间的边界附近。因此, 设计子空间subspace2的上下边界限定为初始全局空间的边界, 即不能超越初始全局空间的边界。

HSRAHO方法通过在设计子空间subspace1和subspace2进行开发和和探索, 以此加快局部收敛速度。设计子空间subspace1通常包含更多的局部最优点, 而subspace2则更专注于探索关键区域。在迭代过程中, subspace1和subspace2可能出现不包含实际全局最优解情况。为了避免HSRAHO方法陷入局部最优, 同样在初始全局空间OGS内进行全局寻优。HSRAHO方法每隔3次迭代分别在OGS、subspace1和subspace2内进行全局寻优, 并通过补充有效样本点来更新组合代理模型和下一次迭代的设计子空间subspace1和subspace2。

3 近似高维优化方法

本文提出的近似高维优化方法HSRAHO联合组合代理模型和多层设计空间缩减策略, 主要包含以下4个部分: ①构造单一代理模型PRS、RBF、KRG和组合代理模型; ②采用多层设计空间缩减策略确定设计子空间subspace1和subspace2;③采用混合自适应有效样本方法(hybrid and adaptive promising sampling, HAPS^[12])增加有效样本点(HAPS方法在下面做详细介绍); ④判断是否满足收敛准则。HSRAHO方法流程如图 1所示, 具体步骤如下:

1) 建立真实分析模型, 确定设计变量和设计空间, 初始化HSRAHO方法参数, 令迭代参数k=1。

2) 使用快速优化拉丁超立方试验设计方法^[15](fast optimal Latin hypercube sampling design, FOLHD)在初始全局空间内生成8个初始样本点, 同时调用真实分析模型计算目标函数值。将初始样本点及其所对应的真实目标函数值保存到样本点数据库中。

3) 提取样本点数据库中的所有样本点及其所对应的真实目标函数值, 分别构造单一代理模型PRS、RBF和KRG。

4) 使用SQP算法求解公式(2)中的优化问题, 计算获得优化权重因子, 构造组合代理模型。公式(2)中参数α和β初始值分别设为0.05和-1。

5) 判断迭代数k是否满足公式rem(k, 3)=2。如果满足, 则计算设计子空间subspace1;否则, 转入步骤7)。rem表示计算余数。

6) 使用HAPS方法在设计子空间subspace1内产生3个有效样本点, 转入步骤10)。

7) 判断迭代数k是否满足公式rem(k, 3)=0。如果满足, 则计算设计子空间subspace2;否则, 转入步骤9)。

8) 使用HAPS方法在设计子空间subspace2内产生3个有效样本点, 转入步骤10)。

9) 使用HAPS方法在初始全局空间OGS内产生3个有效样本点。

10) 将获得的有效样本点和样本点数据库中已有样本点进行比较, 满足公式(6)的有效样本点将被剔除。保存剩余新样本点及其所对应的真实目标函数值到样本点数据库中。

(6)

式中: d表示有效样本点x_p和数据库中样本点x_e在任意维度上的最大距离; U_i^UB和L_i^LB分别表示初始全局空间在第i维上的最大和最小边界值, n为设计变量; η为距离控制系数, 较大的η值会使优化算法错失重要的新样本点, 较小的η值则会产生多余的新样本点, 浪费计算资源, 本文取η值为0.000 1。

11) 重复步骤3)~10)直到满足终止准则, 终止准则表达式如下。

(7)

式中: f_opt表示获得的最优解; f_min表示实际最优解; N_fe表示复杂仿真分析模型调用次数(number of function evaluations, NFE)。当N_fe大于阈值300时, 当前最优结果将作为获得的最优解。

本文提出的HAPS方法操作步骤如下:

1) 提取样本点数据库中的所有“昂贵”样本点, 分别构造代理模型PRS、RBF和KRG, 并分别使用和描述。

2) 基于不同的单一代理模型分别构造4种组合代理模型。为便于描述, AES_PRK表示选用3种单一代理模型PRS、RBF和KRG构造组合代理模型, AES_PR表示选用2种单一代理模型PRS和RBF, AES_PK表示选用2种单一代理模型PRS和KRG, AES_RK表示选用2种单一代理模型RBF和KRG。

3) 在搜索区域内产生数量为G的“便宜”样本点(指采用代理模型计算得到近似目标函数值), 记为X_H¹, X_H², …, X_H^G。本文G设为10 000。

4) 分别使用代理模型PRS、RBF和KRG计算获得“便宜”样本点对应的预测响应值, 分别记为 。

5) 选择预测响应值较小的前g个“便宜”样本点, 分别记为X_H^P1, X_H^P2, …, X_H^Pg; X_H^R1, X_H^R2, …, X_H^Rg; X_H^K1, X_H^K2, …, X_H^Kg。本文g设为100。

6) 分别将样本点X_H^P1, X_H^P2, …, X_H^Pg; X_H^R1, X_H^R2, …, X_H^Rg; X_H^K1, X_H^K2, …, X_H^Kg保存到样本集 。样本集A, B, C中的样本点均取自X_H¹, X_H², …, X_H^G, 因此出现在某一个样本集中的点可能同时出现在另一个样本集中, 甚至包含在所有3个样本集中。

7) 计算获得同时出现在任意2个样本集或3个样本集的样本点, 3个样本集A, B, C将被划分为7个样本子集E~K(如图 2所示), 数学表达式如下

(8)

8) 在所有样本子集E~K中选取数量为3个有效样本点, 各子集分配到的样本点数量S_i由(7)式决定

(9)

式中: i表示各样本子集E~K序号; round表示对结果进行四舍五入; λ_i为反映不同样本子集重要性的权系数, 由各样本子集本身重要性指数l_i和包含的样本点数量r_i共同决定。本文将l_i设为样本子集涉及的代理模型数量, 即l₁=3, l_2-4=2, l_5-7=1。

9) 采用相对应的组合代理模型AES_PRK、AES_PR、AES_PK、AES_RK和单一代理模型 和 分别计算7个样本子集中样本点的预测响应值。从7个样本子集中依次选择预测响应值较小的前S_i个样本点, 作为新增的有效样本点。

4 测试结果与分析 4.1 标准优化函数

选用6个经典的标准优化函数来测试HSRAHO方法的优化性能, 同时选取近似优化方法ESGOHSR^[3]和HAM^[4]进行比较分析。ESGOHSR和HAM方法同样使用了多个代理模型和设计空间缩减策略, 能够有效处理不同类型工程优化问题。Shan和Wang^[1]将设计变量n≥10定义为高维问题。

1) Paviani函数(Paviani), n=10

(10)

2) Trid函数(Trid), n=10

(11)

3) Ellipsoid函数(Ellipsoid), n=12

(12)

4) Sum Squares函数(SumS), n=15

(13)

5) 16维函数(F16), n=16

(14)

式中, α_ij取值参考文献[4]。

(6) Sphere函数(Sphere), n=20

(15)

在实际工程应用中, 优化时间基本上取决于调用真实分析模型次数N_fe, 本文采用N_fe来表征优化效率, 同时采用收敛全局最优解f_opt反映优化精度和鲁棒性。为了减小测试过程中的随机误差, 分别使用3种优化方法HSRAHO、ESGOHSR和HAM连续运行50次, 并对优化结果做统计分析。优化精度与优化效率测试结果对比如表 1至2所示, 分别包括f_opt和N_fe最小值、最大值和平均值, 表中f_opt和N_fe最小平均值采用加粗表示。另外, 表 2中符合“>”表示在N_fe超过阈值max m时, 优化算法仍无法获得满足收敛条件的优化解, 其中括号内数字表示50次重复数值试验过程中无法获得满意优化结果次数。

分析表 1优化精度测试结果可知, HSRAHO方法能够成功捕捉到所有测试函数的收敛优化解, 除Trid和Ellipsoid函数外优化精度均好于比较方法ESGOHSR和HAM。ESGOHSR方法在Trid和Ellipsoid函数中获得的f_opt平均值最接近实际全局最优解, 但在函数Paviani、SumS、F16和Sphere中, 分别有数十次无法获得满足收敛条件的优化解。HAM方法在函数Paviani、Trid、Ellipsoid和F16中, 获得的优化结果平均值最大, 优化精度最低。为了进一步展示HSRAHO方法鲁棒性, 图 3给出了各函数50次优化结果的分布图。从图中可以清楚看出, HSRAHO方法获得的优化解分布均匀, 鲁棒性最强, ESGOHSR方法其次, HAM鲁棒性最弱。优化结果表明, ESGOHSR和HAM方法在搜索过程中容易陷入局部最优, 而HSRAHO方法采样策略很好地平衡了局部开发和全局探索, 有效避免“落入早熟陷阱”。

分析表 2优化效率测试结果可知, HSRAHO方法在除Trid函数外均能以最小计算资源获得收敛优化解。以Paviani函数为例, HSRAHO方法平均只需要调用真实分析模型129.7次, 比ESGOHSR方法平均少超过109.1次, 比HAM方法平均少超过163.7次。另外, ESGOHSR方法在函数Paviani、SumS、F16和Sphere中, 分别有35, 21, 29和34次无法获得收敛优化解。HAM方法则在所有测试函数中仅有数次获得令人满意的优化解。结果表明, 本文提出的多层设计空间缩减策略能够准确识别包含实际全局最优解的有效设计子空间, 而ESGOHSR和HAM方法采用的设计空间缩减策略则无法准确定位。

总之, HSRAHO方法采用多个代理模型和自适应空间探索方法, 在增强优化精度和鲁棒性的同时, 能够有效减轻计算压力, 优化效率明显高于比较方法ESGOHSR和HAM。

4.2 翼型设计优化问题

飞翼布局水下滑翔机作为一种效率高、续航能力强新型无人水下潜航器, 在海洋资源开发、环境监测、深海搜救, 乃至军事应用方面都具有极大研究意义。高升阻比是飞翼布局水下滑翔机的一个重要指标, 设计优化翼型剖面形状成为提高其水动力性能的一个重要手段。本文通过优化设计飞翼布局水下滑翔机的剖面翼型(11维NACA0012翼型)来提高其水动力性能, 同时进一步验证HSRAHO方法的工程实用性。采用形状类别函数变换方法(class function/shape function transformation, CST^[16])描述翼型的几何外形。CST参数化方法表达式如下

(16)

式中: x表示沿翼型弦线方向的横坐标值; y_u0, y_l0, y_u和y_l分别表示基准翼型和目标翼型的上下弧面纵坐标值。参数n表示伯恩斯坦多项式阶数, 本文取n=5。权重系数A_ui和A_li表示翼型优化设计变量。

(17)

考虑到翼型上下弧面曲线的前缘半径应保持一致, 式中A_l0=-A_u0。

以NACA0012为初始翼型, 优化翼型在马赫数Ma=0.5, 攻角α=3°, 雷诺数Re=5×10⁶工况下的升阻比(lift to drag ratio, LDR), 约束条件为翼型的最大相对厚度和面积不小于初始翼型。优化模型如下

(18)

式中: C_l和C_d分别表示翼型升力系数和阻力系数, 由翼型流场仿真软件XFOIL计算得到; s, s₀, t, t₀分别表示优化翼型和基准翼型的最大相对厚度和面积大小。各设计变量的上下边界值由NACA0008和NACA0016翼型采用最小二乘法确定。

使用HSRAHO方法对翼型设计优化问题进行求解, 并与ESGOHSR和HAM方法进行比较。本文采用罚函数^[17]将约束问题转换为非约束问题, 罚因子设为10⁶。优化过程中, 初始样本点数量设为50, 真实分析模型调用最大次数设为300。3种方法对翼型设计优化问题分别连续优化11次, 优化结果中位数见表 3。

由表中结果可知，HSRAHO、ESGOHSR和HAM 3种方法获得的优化翼型最大相对厚度和面积都满足约束，且升阻比相比初始翼型分别提升63.02%，59.61%，50.39%，其中HSRAHO优化翼型升阻比提升最为显著，达到103.65。另外，图 4描绘了3种方法的翼型优化迭代结果。从图中可以清楚看出，相比ESGOHSR和HAM方法，HSRAHO方法能够更快速地收敛到最优翼型，而且占用的计算资源更少。在N_fe到达168前，ESGOHSR方法具有一个较快的收敛速度，但在此之后升阻比并没得到明显提高。HAM方法收敛速度最慢，并且获得的升阻比最小。

图 5给出了初始翼型和HSRAHO优化翼型的外形比较。从图中可以发现，相比初始翼型外形，优化翼型的上弧面厚度增大，下弧面厚度减小，这种翼型外形的变化提升了翼型升阻比。总之，HSRAHO方法在优化效率和全局收敛性方面具有较明显的优势，表明HSRAHO方法具有良好的工程实用性。

5 结论

本文针对高维昂贵黑箱问题求解过程中，传统优化方法需要耗费大量计算时间和计算资源的问题，提出一种基于多层设计空间缩减策略的近似高维优化方法。获得结论如下：

1) 多层设计空间缩减策略打破了传统优化方法依赖经验给定设计空间的局限，通过构建减小的设计子空间，逐步提高组合代理模型在最优解附近的近似精度，最后引导优化快速收敛到全局最优解。

2) 采用6个高维标准优化函数和翼型设计优化实例对HSRAHO方法的性能进行测试，并与同样采用设计空间缩减策略的优化方法ESGOHSR和HAM进行比较分析，结果表明HSRAHO方法在全局收敛性、优化效率和鲁棒性方面具有明显优越性。

3) HSRAHO方法在处理实际高维工程应用还有待进一步研究。在未来工作中，考虑引入维度降低策略。