0 引　言

水动力外形优化设计是水下航行器（UUV、潜艇、鱼雷等）研制的重要环节之一，对动力系统、控制系统和能源系统设计有重要影响，其结果直接影响潜器运动性能的优劣^[1]。在实际的工程优化问题中，高可信度的目标和约束函数值很难直接获得，通常需要求解计算流体动力学（CFD）模型、有限元（FEM）模型等耗时长的模拟程序，整个优化过程耗费大量的计算资源且费时费力。针对这一问题，基于代理模型的优化方法（SBO）^[2]通过试验设计选取一定数量的样本点来构造目标和约束函数的近似模型，将工程优化问题转变为数值优化问题，可以显著提高优化效率，降低优化难度，被广泛应用于航空航天、船舶等复杂系统领域。

但对于复杂的多变量水动力外形优化问题，需要大量样本点才能得到满足工程优化要求的代理模型，所需的时间和计算成本仍然较高。为了进一步提高优化效率，平衡计算效率和代理模型精度之间的矛盾，学者们进行了变保真度模型（VFM）^[3-5]，序列近似优化（SAO）^[6-8]等基于代理模型的优化方法研究。其中，SAO方法基于少量初始样本点构造的初始代理模型，采取一定的加点准则更新样本点，逐步提高代理模型对全局最优解的近似精度，可在保证全局最优的前提下尽可能减少样本点的数量，受到了广泛关注。

代理模型构造和加点准则是SAO方法的两大关键技术，对算法效率和优化精度具有显著影响^[9]，国内外学者对此开展了大量工作。谢延敏等^[10]基于EI加点准则提出一种改进并行加点策略，可以有效提高Kriging模型的精度；柳斐等^[11]提出一种一次添加多个样本带点的多点EI准则（MEI），可以显著提高气动设计优化效率；Viana等^[12]提出利用多个代理模型同时添加多个EI最大值的加点准则（MSEGO），并验证了算法有效性。陈国栋等^[13]运用信赖域更新技术，提出基于序列径向基函数的优化方法，将该方法应用于车身结构轻量化问题中。Wang D等^[14]提出了一种混合序列优化方法（HSAO），并应用于S形管的耐撞性优化问题。彭科等^[15]将潜在可行域最大距离加点准则和潜在最优加点准则相结合，以Golinski减速器和“天航二号”火箭为例，验证了所提算法具有良好的全局寻优能力与搜索效率。

本文基于最大最小距离准则（MD）和最小化代理模型预测准则（MP）提出一种多点加点准则，建立二步收敛判定准则，以平衡算法的全局探索和局部开发能力，构建改进序列近似优化的算法流程。开展二维Rosenbrock函数和Golinski减速器优化算例的研究，并与传统优化方法进行对比，验证了改进序列近似优化方法的有效性和优越性。最后将该方法应用到某型水下航行器的水动力外形优化问题中，取得了预期的效果。

1 序列近似方法及其改进 1.1 序列近似优化方法基本流程

一般地，优化问题的数学模型可描述为下列形式：

$ \left\{ \begin{array}{*{20}{l}} {\text{min}}&f(x,y(x)),\\ {\text{s.t.}}&{g_i}(x,y(x)) \leqslant 0,\;\;i = 1, \cdots ,p,\\ {\text{ }}&{h_j}(x,y(x)) = 0,\;\;j = 1, \cdots ,q,\\ {\text{ }}& {x_L} \leqslant x \leqslant {x_U}\text{。} \end{array} \right. $

(1)

其中，x为设计变量向量；x_U和x_L分别为其上、下边界；f为目标函数向量；g为不等式约束函数向量；h为等式约束函数向量；p和q分别为不等式约束和等式约束的数目；y(x)为设计变量对应的状态响应。

对于一个实际的工程优化问题，序列近似优化方法的基本流程如图1所示，具体的优化步骤为：

步骤1　建立优化问题的数学模型，确定设计变量、设计目标、设计约束和设计空间；

步骤2　采用给定的试验设计方法选取初始样本点，并计算初始样本点的目标和约束函数值，构造初始代理模型；

步骤3　采用给定的优化算法求解当前代理模型的全局最优解，并用高可信度的模拟程序对最优解进行检验；

步骤4　如果当前最优解满足收敛判定规则，则计算终止，输出最优解，否则执行下一步；

步骤5　根据优化加点准则，生成新样本点，并计算新样本点的目标和约束函数值；

步骤6　将新样本加入到初始样本中，构造新的代理模型，然后执行步骤3，进行下一轮迭代，直至满足收敛判定准则。

1.2 加点准则

为了平衡算法的全局探索和局部开发能力，将最大最小距离准则（MD）^[16]和最小化代理模型预测准则（MP）^[17]结合起来，各取所长。

MD准则是一种均匀性准则，其基本思想是依据式(2)求解设计空间内一点x，使x与包含所有样本的点列 $ X = \{ {x_1},{x_2}, \cdots ,{x_n}\} $ 之间的最小欧式距离最大化，将式(2)的最优解作为新的样本点，最大限度保证算法探索能力。MD准则仅用到样本的空间位置信息，与其状态响应无关，因此在一次样本更新过程中，可反复迭代产生多个新样本点。

$ d(x) = \mathop {\max }\limits_{x \in {R^s}} \mathop {\min }\limits_{1 \leqslant i \leqslant n} d({x_i},x)。$

(2)

MP准则的原理是将式(3)的全局最优解 $ {\hat f_{out}} $ 作为新增样本点，以提高代理模型对最优解的近似精度。若单独采用MP准则，对于复杂的优化模型往往会陷入局部最优解。

$ \left\{ {\begin{array}{*{20}{l}} {{\hat f}_{out}} = \min \hat f(x,y(x)),\\ {\text{ s.t.}}\;\;{g_i}(x,y(x)) \leqslant 0,\;\;i = 1, \cdots ,p ,\\ \qquad{h_j}(x,y(x)) = 0,\;\;j = 1, \cdots ,q, \\ \qquad{x_L} \leqslant x \leqslant {x_U}\text{。} \end{array}} \right. $

(3)

式中， $ \hat f(x,y(x)){\text{ }} $ 为本轮迭代的代理模型。

图2给出了一维问题的MD和MP加点准则的原理图。可以看出，2种加点准则仅与样本分布和代理模型全局优化解有关，因此不受代理模型约束，具有通用性。

1.3 收敛判定准则

为了在保证全局优化精度的前提下提高算法的收敛速率，定义二步收敛判定准则如式(4)和式(5)所示。首先，同时使用MD和MP准则，每轮迭代增加n+1个样本点，以提高代理模型的整体精度和确定全局最优解所在区域为目标；当满足式(4)后，仅采用MP准则更新样本点，以代理模型最优解 $ {\hat f_{out}} $ 与最优解处的真实响应 $ {f_{out}} $ 之间的相对误差不大于容差T_e为优化终止条件。

$ \left\{\begin{array}{*{20}{l}} \mathrm{max}\dfrac{\left|\right|{\widehat{x}}_{opt}^{i}-{\widehat{x}}_{opt}^{i-1}\left|\right|}{\left|\right|{\widehat{x}}_{opt}^{i}\text+\epsilon \left|\right|}\leqslant {T}_{x}\text{，}\varepsilon \text=\left\{\begin{array}{cc}0& \left|\right|{x}_{opt}^{i}\left|\right| > {x}_{0}\text{，}\\ 0.01& \left|\right|{x}_{opt}^{i}\left|\right|\leqslant {x}_{0}\text{，}\end{array}\right.\\ {R}^{2}\geqslant {R}_{0}^{2},\\ {N}_{\mathrm{max}}\geqslant N\geqslant {N}_{\mathrm{min}}\text{。}\end{array}\right. $

(4)

式中： $ i = N - {N_{\min }},N - {N_{\min }} + 1, \cdots ,N $ ；N为当前迭代步数；N_max为最大迭代步数；N_min为最小迭代步数， $ {\hat x_{opt}} $ 为 $ {\hat f_{out}} $ 对应的样本点空间位置；R²为拟合优度，其值越接近1表明近似模型的拟合精度越高。

$ \left\{ {\begin{array}{*{20}{l}} {\dfrac{{||{f_{opt}} - {{\hat f}_{opt}}||}}{{||{f_{opt}}{\text{ + }}\delta ||}} \leqslant {T_e},{\text{ }}\delta {\text{ = }}\left\{ {\begin{array}{*{20}{c}} 0&{||{f_{opt}}|| > {f_0}} ,\\ {0.00{\text{0}}1}&{||{f_{opt}}|| \leqslant {f_0}} , \end{array}} \right.} \\ {{N_{\max }} \geqslant N{\text{ }}}\text{。} \end{array}} \right. $

(5)

1.4 改进序列近似优化方法基本流程

首先将设计空间映射s维超立方体内，对样本点进行归一化，采用优化拉丁超立方试验设计方法（OLHD）选取初始样本点，代理模型选择径向基函数（RBF）模型，得到改进序列近似优化方法的基本流程如图3所示。

2 算例测试

为比较改进序列近似优化方法与传统优化方法的优化效率和精度，选择2个经典算例进行测试。

2.1 二维Rosenbrock函数

二维Rosenbrock函数^[18]为单峰、非凸的病态函数，优化难度较大，是无约束优化的经典算例，优化的数学模型为：

$ \left\{ {\begin{array}{*{20}{c}} {\min {\text{ }}f({x_1},{x_2}) = {{(1.0 - {x_1})}^2} + 100.0{{({x_2} - {x_1}^2)}^2}} ,\\ {{{\rm{s}}} .{\rm{t}}.{\text{ }}{x_1},{x_2} \in [ - 2.0,2.0]{\text{。}}} \end{array}} \right. $

(6)

其函数图像如图4所示，最优解为(x₁^*,x₂^*)=(1.0,1.0)，f(x₁^*,x₂^*)=0.0，位于一个狭长的峡谷中。

首先将设计变量取值均转化到[0,1]区间，根据OLHD方法选取10个初始样本点，构造初始RBF代理模型。设置改进SAO方法最小迭代步数N_min为3，最大迭代步数N_max为100，最小拟合优度 $ R_0^2 $ 为0.96，容差 $ {T_x} $ 和 $ {T_e} $ 分别为0.05和0.02，MD准则每个迭代步生成2个新样本点。

图5为序列优化过程中目标函数真实值 $ {f_{{\rm{out}}}} $ 和拟合优度 $ {R^2} $ 的迭代曲线。可以看出，随着样本点增加，代理模型的拟合精度不断提高，目标函数真实值逐渐逼近全局最优解，与期望一致。但目标函数取得最优解时，优化并未终止，这是因为Rosenbrock函数的全局最优解为0，此时可能绝对误差很小，相对误差很大，因此，若可以预见优化问题的最优解在0附近，可将优化终止条件定义为绝对误差。最终改进SAO方法迭代32次后收敛，优化过程中调用实际函数88次，最优解为3.13E-5，收敛后样本在设计空间的分布如图6所示。

将各优化方法随机重复10次，平均解为去除最大值和最小值后的平均值，得到优化结果如表1所示。表中， $ {N_{Sam}} $ 为样本总数， $\Delta {{\boldsymbol{x}}_{opt\_best}}$ 为最优解对应的样本点在单位正方形内空间位置与全局最优解之间的欧式距离。可见相比SBO方法，本文方法可以较少的样本点获得较高的优化精度，提高优化效率。

2.2 Golinski减速器优化问题

Golinski减速器优化问题^[19]是NASA提出的多学科优化设计测试算例，该算例含7个变量，11个约束，理论上最优解为2994.0，优化的数学模型为：

$ \begin{split}\left\{ {\begin{array}{*{20}{l}} \min \,f({\boldsymbol{x}}) = 0.7854{x_1}{x_2}^2(3.3333{x_3}^2 + 14.9334{x_3} -\\ \qquad\qquad\;\; 43.0934) - 1.5079{x_1}({x_6}^2 + {x_7}^2)+ \\ \qquad\qquad\;\; 7.477({x_6}^3 + {x_7}^3) + 0.7854({x_4}{x_6}^2 + {x_5}{x_7}^2)，\\ {\rm{s.t.}}\;{g_1}{\text{(}}{\boldsymbol{x}}) = 27.0/({x_1}{x_2}^2{x_3}) - 1 \leqslant 0，\\ \quad\;\;{g_2}{\text{(}}{\boldsymbol{x}}) = 397.5/({x_1}{x_2}^2{x_3}^2) - 1 \leqslant 0，\\ \quad\;\;{g_3}{\text{(}}{\boldsymbol{x}}) = 1.93{x_4}^3/({x_2}{x_3}{x_6}^4) - 1 \leqslant 0，\\ \quad\;\;{g_4}{\text{(}}{\boldsymbol{x}}) = 1.93{x_5}^3/({x_2}{x_3}{x_7}^4) - 1 \leqslant 0，\\ \quad\;\;{g_5}{\text{(}}{\boldsymbol{x}}) = \sqrt {{{(745.0{x_4}/{x_2}{x_3})}^2} + 16.9 \times {{10}^6}} / \\ \;\;\qquad\qquad 0.1{x_6}^3 - 1100 \leqslant 0 ，\\ \quad\;\;{g_6}{\text{(}}{\boldsymbol{x}}) = \sqrt {{{(745.0{x_5}/{x_2}{x_3})}^2} + 157.5 \times {{10}^6}}/ \\ \;\;\qquad\qquad 0.1{x_7}^3 - 850 \leqslant 0，\\ \quad\;\;{g_7}{\text{(}}{\boldsymbol{x}}) = {x_2}{x_3} - 40 \leqslant 0 ，\\ \quad\;\;{g_8}{\text{(}}{\boldsymbol{x}}) = {x_1}/{x_2} - 5 \geqslant 0 ，\\ \quad\;\;{g_9}{\text{(}}{\boldsymbol{x}}) = {x_1}/{x_2} - 12 \leqslant 0 ，\\ \quad\;\;{g_{10}}{\text{(}}{\boldsymbol{x}}) = (1.5{x_6} + 1.9)/{x_4} - 1 \leqslant 0 ，\\ \quad\;\;{g_{11}}{\text{(}}{\boldsymbol{x}}) = (1.1{x_7} + 1.9)/{x_5} - 1 \leqslant 0 ，\\ \quad\;\;{x_1} \in {\text{[2}}{\text{.6,3}}{\text{.6], }}{x_2} \in {\text{[0}}{\text{.7,0}}{\text{.8], }}{x_3} \in {\text{[17,28]}}(Z) ，\\ \quad\;\;{x_4},{x_5} \in {\text{[7}}{\text{.3,8}}{\text{.3], }}{x_6} \in {\text{[2}}{\text{.9,3}}{\text{.9], }}{x_7} \in {\text{[5}}{\text{.0,5}}{\text{.5]}} 。\end{array}} \right. \\[-150pt]\end{split}$

(7)

本算例中含有离散变量，因此将设计变量取值均转化到[0,11]区间。选取初始样本个数为15，MD准则每迭代步新增7个样本点，优化过程中目标函数真实值的迭代曲线如图7所示。本文方法迭代16次之后收敛，搜索到全局最优解为 [3.5, 0.7, 17, 7.3, 7.715320, 3.350215, 5.286655]。

表2列出了改进SAO方法与其他文献对Golinski减速优化问题的优化结果，表明改进SAO方法也可以收敛到含约束问题的全局最优解，而且具有较高的优化效率。

3 水下航行器水动力外形优化设计 3.1 优化问题的数学模型

本文研究的某型水下航行器水动力外形由艇体、首升降舵和尾部鳍舵组成，其中首升降舵为梯形，尾鳍为切尖三角形，尾升降舵和方向舵为矩形，均采用NACA0015剖面翼型，如图8所示。

为便于优化过程中模型更新，建立水下航行器外形的参数化模型如图9所示。

其中，b_(i)，C_r(i)，C_t(i)分别为尾鳍的展长，根弦长和梢弦长，(i)取b，s，r分别表示首升降舵、尾水平鳍和尾垂直鳍。

采用Nystrom线型^[22]作为艇体进流段和去流段母线线型，得到艇体回转半径r(x) 表达式为：

$ r(x) = \left\{ {\begin{array}{*{20}{l}} R{{\left[ {1 - {{\left( {\dfrac{{{L_1} - x}}{{{L_1}}}} \right)}^{{n_n}}}} \right]}^{1/{n_n}}},&0 \leqslant x < {L_{1}} ，\\ R,&{L_1} \leqslant x < {L_1} + {L_2}，\\ R\left[ {1 - {{\left( {\dfrac{{x - {L_1} - {L_2}}}{{{L_3}}}} \right)}^{{n_t}}}} \right],&{L_1} + {L_2} \leqslant x \leqslant {L_1} + {L_2} + {L_3}。\end{array}} \right. $

(8)

在水下航行器设计过程中，需要考虑快速性、操纵性、结构、能源、推进、操纵与控制、总布置等诸多因素，对航行器的安全性能和可靠性有重要意义^[23]。一般地，巡航阻力系数越小，快速性越好。本文在前人研究工作^[24]的基础上，根据总体性能要求和任务特点，考虑艇体快速性和内部布置空间大小，以巡航阻力系数和空间损失率最小化为目标，确定水下航行器水动力外形优化的参数变量及其取值空间如表3所示。

采用线性加权方法，将多目标优化问题转化为单目标优化问题。根据表3，该优化问题包含7个变量，6个不等式约束，11个等式约束，将各变量的取值均转化到[0,1]区间上，得到其数学模型为：

$ \left\{ \begin{array}{*{20}{l}} \min &f({\boldsymbol{x}}) = {w_1}{C_d}({\boldsymbol{x}}) + {w_2}\xi ({\boldsymbol{x}})，\\ {\text{s.t.}}&{{\boldsymbol{g}}_i}({\boldsymbol{x}},y({\boldsymbol{x}})) \leqslant 0,\;\;i = 1, \cdots ,6 ，\\ {\text{ }}&{{\boldsymbol{h}}_j}({\boldsymbol{x}},y({\boldsymbol{x}})) = 0,\;\;j = 1, \cdots ,11，\\ {\text{ }}&{\boldsymbol{x}} \in [0,1],\;\;{w_1} = 0.8,{\text{ }}{w_1} 0.2。\end{array} \right. $

(9)

其中，C_d和ξ的计算表达式为：

$ {C_d} = {{2{F_d}} \mathord{\left/ {\vphantom {{2{F_d}} {\rho S{v^2}}}} \right. } {\rho S{v^2}}}, $

(10)

$ \xi {\text{ = }}1 - Vol/(L \times B \times H)\text{。} $

(11)

式中：F_d为巡航阻力；ρ为流体介质的密度；S为艇体的截面面积。

3.2 数值计算方法

水下航行器巡航阻力通过STAR-CCM+数值仿真计算获得，流体计算域及边界条件如图10所示。其中，计算域为92 m×40 m×40 m的长方体，左侧和周围侧面为速度入口边界，入口速度设定为[1.0288，0，0] m·s⁻¹，右侧为压力出口边界，出口压力设定为 0 Pa，航行器表面边界为壁面条件，为了合理布置网格，航行器周围为局部加密区域。

流域网格划分采用切割体网格生成器和棱柱层网格生成器以有效模拟航行器直航运动时的流场流动。为了验证网格无关性，确定合适的网格数量，使用多种网格基本尺寸计算来流为2 kn时的巡航阻力系数，计算结果如表4所示。可以看出，网格基本尺寸小于0.72 m时，就可以获得较高的计算精度，根据计算资源，最终选用358.0万网格进行航行器巡航阻力系数仿真计算。

3.3 优化结果与分析

应用OLHD方法从设计空间随机选取15个初始样本点，计算样本的目标和约束函数值，建立初始RBF代理模型，并将优化结果作为MP准则的初始输入。收敛参数设置同2.1节，为提高优化效率，MD准则每迭代步生成7个样本点。利用Isight集成CATIA，STAR-CCM+和Matlab等商业软件搭建改进序列近似优化仿真流程，如图11所示。

优化过程目标函数真实值的迭代曲线如图12所示，迭代至第10步，仅调用原始计算模型80次后优化收敛，大大提高了设计效率。

经过对优化结果的综合对比与分析，获得了某型水下航行器的最佳水动力外形。图13为优化前后航行器水动力外形轮廓的对比，图14为艇体中间剖面的压力分布对比。表5为优化方案与原始方案设计变量和性能参数的对比，其中，初始方案艇体设计变量取值为文献[24]中优化结果。优化后所得设计方案与初始设计方案相比，巡航阻力系数降低了9.81%，空间损失率降低了0.28%，航行器在不减小艇体内部布置空间的前提下，水动力性能得到大幅提升。

可以看出，相比初始方案，优化方案的艇体首部更钝，逆压梯度更大，巡航阻力系数增加，空间损失率减小。艇体尾部更尖锐，逆压梯度减小，巡航阻力系数减小，空间损失率增加。尾部鳍舵的展弦比和梢根比分别降低了42.4%和78.32%，其外露面积减小，前缘后掠角增大，巡航阻力系数降低。因此，优化结果表明减小巡航阻力系数，提升快速性和减小空间损失率，增大艇体内部布置空间二者是相互矛盾的，必须综合平衡。

4 结　语

本文对序列近似优化算法进行了改进，通过2个经典算例测试了算法性能，最后将该方法应用于某型水下航行器水动力外形优化问题中，得到以下结论：

1）基于最大最小距离准则（MD）和最小化代理模型预测准则（MP）提出了一种多点加点准则，兼顾了算法的全局探索能力和局部开发能力，适用于所有代理模型。

2）通过二维Rosenbrock函数和Golinski减速器优化算例，并与SBO方法及其他文献结果进行对比，验证了本文方法的有效性和优越性，可以较少的样本点获得较高的优化精度，提高优化效率。

相比初始方案，优化方案的巡航阻力系数和空间损失率分别降低了9.81%和0.28%，航行器在不减小艇体内部布置空间的前提下，水动力性能得到大幅提升。