0.   引　言

随着船舶结构复杂程度的不断提升，高精度的有限元仿真软件已成为船舶结构分析的主流工具，而直接调用耗时的仿真分析程序进行优化已难以满足船舶设计快速高效的需求，因此，发展基于代理模型的高效优化方法是船舶结构优化设计的重点研究领域之一^[1]。所谓代理模型，是指对一系列通过仿真分析获得的输入/输出样本进行拟合得到的统计近似模型。在代理模型提供的信息引导下的优化可以极大地减少耗时的仿真次数，达到缩短设计周期、节省设计成本的目的。目前，经典的代理模型建模与优化方法已被广泛用于船舶领域，例如对船舶板架强度和稳定性^[2]、潜艇振动声辐射^[3]、水下结构物基座阻抗特性^[4]等响应的快速预报，以及对深潜器的多球交接耐压壳^[5]、多用途船货舱段^[6]的结构优化和对海上发射船^[7]、集装箱船^[8] 等的船型优化。

考虑将梯度信息纳入到用于构建代理模型的样本中，通常可以显著提高模型的精度^[9]。而在众多的梯度增强代理模型中，梯度增强Kriging (Gradient-enhanced Kriging, GEK)模型^[10]由于其对非线性问题具有良好的插值能力和独特的概率模型属性，所以在一些领域的应用效果明显，例如飞机机翼/翼型的气动优化设计^[11]、悬臂梁与简支梁的结构优化设计^[12]、轴孔过盈配合的接触应力分析^[13]、核反应堆最大燃料温度的不确定性量化^[14]等。然而，在船舶优化领域，对实现高效梯度计算方法（例如伴随法^[15]）的工具开发尚不成熟，使用梯度增强Kriging模型分析的案例并不多见，因而，利用该模型的优势进行设计优化的方法仍有待发展。

本文拟将基于梯度增强Kriging代理模型的优化方法^[16]扩展至梯度信息只能通过有限差分法获取的场景，实现对船舶结构的高效优化设计。即提出基于缩减型梯度增强Kriging模型的加点策略，仅对有实际改进的采样位置进行梯度计算，以此减少仿真调用的次数。在加点过程中，首先采用多起点局部优化算法搜索改进期望函数的若干局部最优解作为候选加点位置，然后计算相应的近似驻点概率，并根据改进期望值和近似驻点概率值的一致程度来确定加点位置，达到提高优化效率的目的。最后，以某水下航行器的结构优化设计为对象，对所提方法进行验证。

1.   梯度增强Kriging模型

1.1   标准型梯度增强Kriging模型

Kriging模型的基本假设为：将确定性响应$y\left( {\boldsymbol{x}} \right)$视为高斯随机过程$Y\left( {\boldsymbol{x}} \right)$。不失一般性，描述为

式中，$f\left( {\boldsymbol{x}} \right)$为包含若干回归项的向量；${\boldsymbol{\beta }}$为待估计的回归系数向量；平稳随机过程$Z\left( \cdot \right)$的均值为0，协方差$Cov\left[ {Z\left( {\boldsymbol{x}} \right),Z\left( {{\boldsymbol{x}}'} \right)} \right] = {\sigma ^2}Corr\left[ {Z\left( {\boldsymbol{x}} \right),Z\left( {{\boldsymbol{x}}'} \right)} \right]$，其中${\sigma ^2}$为过程方差，$Corr\left[ {Z\left( {\boldsymbol{x}} \right),Z\left( {{\boldsymbol{x}}'} \right)} \right]$为相关函数。对于梯度增强Kriging，给定采样矩阵${\boldsymbol{X}}$，其增广后的样本数据向量${\tilde{\boldsymbol{y}}}$中既包含了响应也包含了梯度，记为

式中：$n$为样本数；$k$为设计变量数。本文统一用带括号的上标表示样本个体编号，下标表示设计变量序号。根据初始样本${{S}} = \left\{ {{\boldsymbol{X}},{\tilde{\boldsymbol{y}}}} \right\}$，对${\boldsymbol{\beta }}$，${\sigma ^2}$和隐含在相关函数中的超参数${\boldsymbol{\theta }}$进行极大似然估计。似然函数为多元正态分布的联合概率密度函数，即

式中：${\tilde{\boldsymbol{F}}}$，${\tilde{\boldsymbol{R}}}$分别为增广的回归矩阵和相关函数矩阵，分别为：

作为$\left( {n + nk} \right) \times \left( {n + nk} \right)$对称矩阵的${\tilde{\boldsymbol{R}}}$，可依照式中所示分块建立。对于$n \times n$的子矩阵${\boldsymbol{R}}$，$n \times nk$的子矩阵${{\partial {\boldsymbol{R}}} / {\partial {\boldsymbol{x}}}}$以及$nk \times nk$的子矩阵${{{\partial ^2}{\boldsymbol{R}}} / {\partial {{\boldsymbol{x}}^2}}}$，其元素分别为响应之间、响应与偏导数之间和偏导数之间的相关函数。当采用高斯型相关函数时，具体表达式可详见文献[16]。

令式(3)关于${\boldsymbol{\beta }}$和${\sigma ^2}$的偏导分别为0，可得到极大似然估计量分别如下：

由于超参数${\boldsymbol{\theta }}$没有解析解，所以将式(6)和式（7）代回到式(3)，通过优化算法搜索似然函数最大值所对应的超参数作为估计值。为了简化描述，直接给出梯度增强Kriging在任意采样点${\boldsymbol{x}}$处的预测式和均方预测误差，并分别表示如下：

式中，${\tilde{\boldsymbol{ u}}} = {{\tilde{\boldsymbol{F}}}^{\rm{T}}}{{\tilde{\boldsymbol{R}}}^{ - 1}}{\tilde{\boldsymbol{r}}} - f\left( {\boldsymbol{x}} \right)$，其中${\tilde{\boldsymbol{r}}}$为$\left( {n + nk} \right) \times 1$的相关函数向量，其元素的表达式可参见文献[16]。

1.2   缩减型梯度增强Kriging模型

考虑仅将一部分梯度数据包含到样本数据向量中，可构建缩减型梯度增强Kriging模型。针对提高构建梯度增强Kriging模型效率的问题，Han等^[17]建立了一系列考虑一组梯度信息的子模型，并加权求和；Chen等^[18]筛选出了与响应显著相关的变量所对应的偏导数，将其包含到样本中。本文仅考虑在部分样本点处计算梯度，将其包含到样本数据向量中，建立的缩减型梯度增强Kriging模型可视为文献[17]所述子模型的普遍形式。

将采样矩阵${\boldsymbol{X}}$按行重新排列为

式中，${{\boldsymbol{X}}_0}$和${{\boldsymbol{X}}_1}$分别为${n_0} \times k$和${n_1} \times k$的矩阵，且${n_0} +$${n_1} = n$。仅在${{\boldsymbol{X}}_1}$中包含的采样位置处计算梯度，相应地，增广的样本数据向量${\tilde{\boldsymbol{y}}'}$写为

构建缩减型梯度增强Kriging模型的过程及表达形式与标准型梯度增强Kriging模型基本一致，本文不再赘述。但需要注意的是，此时增广的相关函数矩阵${\tilde{\boldsymbol{R}}'}$为$\left( {n + {n_1}k} \right) \times \left( {n + {n_1}k} \right)$，增广的极大似然估计量${\hat \sigma '^2}$的分母为$\left( {n + {n_1}k} \right)$，增广的相关函数向量${\tilde{\boldsymbol{ r}}'}$为$\left( {n + {n_1}k} \right) \times 1$。增广的协方差矩阵${\hat \sigma '^2}{\tilde{\boldsymbol{R}}'}$仍保持正定，因为对于任意非零向量${\boldsymbol{v}} \in {{\bf{R}}^{\left( {n + {n_1}k} \right)}}$，有

式中，${\tilde{\boldsymbol{ Y}}'}$为增广的样本数据向量${\tilde{\boldsymbol{y}}'}$所对应的随机变量向量。

2.   基于梯度增强Kriging的优化方法

2.1   近似驻点概率

工程问题对应的真实函数通常是光滑的^[19]，即函数在设计空间内处处可导。对于任意非单调的函数，若其全局最优不在设计空间的边界上，则必然位于某驻点处。由于梯度增强Kriging方法可以对梯度进行最优线性无偏预测，所以可用来估计任意点为真实函数的驻点概率。根据式（1）可知，对于梯度，有

这里，以真实函数关于第$l$维设计变量($l = 1, \cdots ,k$)的偏导为例，考虑线性预测式如下：

式中，${{\tilde{\boldsymbol{ c}}}_l}\left( {\boldsymbol{x}} \right)$为$\left( {n + nk} \right) \times 1$的权重系数向量。为得到${{\partial y}/ {\partial {x_l}}}\left( {\boldsymbol{x}} \right)$的最优线性无偏预测式，应求得一组${{\tilde{\boldsymbol{ c}}}_l}\left( {\boldsymbol{x}} \right)$，使预测的均方误差（mean-squared error, MSE）最小。

同时，满足如下无偏约束：

对于求解过程，可参见文献[16]，本文不再赘述。

由此，梯度增强Kriging方法对偏导数的预测式和预测均方误差可分别表示为：

式中，${{\boldsymbol{u}}_l} = {{\tilde{\boldsymbol{F}}}^{\rm{T}}}{{\tilde{\boldsymbol{R}}}^{ - 1}}{{\tilde{\boldsymbol{r}}}_l} - {{\partial f} /{\partial {x_l}}}\left( {\boldsymbol{x}} \right)$。

根据Kriging的概率模型属性，预测的梯度服从正态分布

因此，第$l$维偏导${{\partial y} / {\partial {x_l}}}\left( {\boldsymbol{x}} \right)$($l = 1, \cdots ,k$)属于区间$\left[ { - {\varepsilon _l}, + {\varepsilon _l}} \right]$的概率，表示为

式中，$f\left( {{{\partial y} / {\partial {x_l}}}\left( {\boldsymbol{x}} \right)} \right)$为式(19)所述正态分布的概率密度函数；d表示在积分式中${{\partial y}}/{{\partial {x_l}}}\left( {\boldsymbol{x}} \right)$为积分变量。图1中阴影部分的面积表示${P_l}$的值，${{\partial \hat y} / {\partial {x_l}}}\left( {\boldsymbol{x}} \right) = 0$的情况类似（虚线对称轴与y轴重合），故省略。

图  1  概率${P_l}$的几何表示

Figure  1.  Illustrations for the probability ${P_l}$

下载: 全尺寸图片 幻灯片

当选取接近0的正实数${\varepsilon _l}$($l = 1, \cdots ,k$)作为区间边界时，近似驻点概率(approximate probability of stationary point, APSP)可定义为

特别地，对于已知的样本个体$\left( {{{\boldsymbol{x}}^*},{y^*}} \right) \in {{S}}$，令

以保证近似驻点概率在整个设计空间内有连续定义。同样地，对于缩减型梯度增强Kriging模型，近似驻点概率的计算过程类似，仅需替换相应的预测式和均方预测误差。

2.2   加点准则

对于梯度增强Kriging代理模型的优化方法，经典的EGO （efficient global optimization）算法^[20]仍然适用。然而，此时利用梯度信息的情况仍限于建模阶段，若样本容量较少，梯度增强Kriging方法对响应的预测结果可能会出现一定程度的偏差，影响改进期望（expected improvement, EI）函数的判断。这里，以使用基于梯度增强Kriging模型的EGO算法对如下一维函数进行全局优化为例进行说明。

通过拉丁超立方采样^[21]得到4个初始样本，计算各样本点处的梯度，建立梯度增强Kriging代理模型，其预测值和均方预测误差如图2所示。

图  2  梯度增强Kriging对一维函数值的预测

Figure  2.  Prediction of the one-dimensional function by Gradient-enhanced Kriging

下载: 全尺寸图片 幻灯片

根据EGO算法的原理，应选择改进期望函数的最大值所在位置作为下一轮采样点。然而，从图3(a)中可见，当前改进期望函数的最大值位于菱形标识处，而真实函数的全局最优解位于远离于此的星形标识处。因此，当前的加点效率因改进期望函数的偏离而受限。若同时考虑近似驻点概率值的大小，如图3(b)所示，则真实函数的全局最优解所对应的近似驻点概率值将高于改进期望最大值的所在位置，即有机会选择到全局最优解附近加点，从而使优化效率得到提高。

图  3  改进期望和近似驻点概率在当前加点阶段选择更新点时的决策对比

Figure  3.  Comparison of the expected improvement and APSP in selecting the infill point at the current infill stage

下载: 全尺寸图片 幻灯片

通常，改进期望函数具有多峰特性。Zhan等^[22]提出搜索改进期望函数的多个局部最优位置，以并行的方式加点。当不具备并行条件时，需要合理地选择加点位置。本文提出了基于缩减型梯度增强Kriging的优化方法，具体步骤如下：

1)使用拉丁超立方采样确定采样矩阵${\boldsymbol{X}}$，计算相应位置的响应与有限差分梯度，得到初始样本${{S}} = \left\{ {{\boldsymbol{X}},{\tilde{\boldsymbol{y}}}} \right\}$，建立缩减型梯度增强Kriging模型。

2)使用拉丁超立方采样产生$Q$个位置作为内点法^[23]的初始点，对改进期望函数进行搜索，得到包含$Q$个局部最优解（允许存在重复）位置的集合$H$。

3)对集合$H$中元素按对应的改进期望值大小降序排列，得到排序后的集合${H_{\rm{sorted}}}$。

4)计算集合${H_{\rm{sorted}}}$中元素所对应的近似驻点概率。

5)计算集合${H_{\rm{sorted}}}$的后($Q - 1$)个元素所对应的一致程度差值${\varDelta _q}$ ($q = 2, \cdots ,Q$)，计算式为

6)若所有差值${\varDelta _q}$均不小于0，则在集合${H_{\rm{sorted}}}$的第1个元素处计算响应；否则，在集合${H_{\rm{sorted}}}$中对应的最小差值${\varDelta _q}$的元素处计算响应。若新计算的响应小于当前样本中的最优解，则通过有限差分计算当前位置的梯度；否则，不计算任何位置的梯度。

7)判断是否满足终止条件，若满足，停止加点，输出优化结果；否则，将新采样的样本响应和梯度补充到样本集合中，更新梯度增强Kriging模型，重复步骤2~步骤6。

3.   算　例

水下航行器由于受到自身发动机工作的影响，外壳会发生振动并激励外场的海水介质形成辐射声场，所以分析水下航行器自身振动特性是研究其辐射声场强度分布的基础。航行器在水下的振动湿模态分析比其在空气中的振动干模态分析更复杂：一方面，水下航行器航行时在周围海水介质的作用下会产生流固耦合效应，无法对其湿模态进行解析求解； 另一方面，水下试验难度大、成本高，对于全尺寸水下航行器进行振动试验也很难完成。因此，利用有限元分析软件对水下航行器进行湿模态分析，可以高精度地模拟其在水下振动时的流固耦合特性，从而为后续的结构优化设计提供分析手段。

本文以某水下航行器的简化结构优化设计为例，验证所提优化方法的可行性。为了模拟航行器在水中的悬浮状态，不对其施加任何位移约束和力载荷，其在水下自由振动时的前6阶模态均为刚体模态（计算数值均接近于0），而将提高其第7阶固有频率作为优化目标，各设计变量如图4所示。待优化的水下航行器基准形的平行舯体直径$\phi = 2.54$ m，平行舯体长$L = 10.16$ m，艏径$R1 =$$4.52$ m，艉径$R2 = 9.42$ m，相应的取值范围如表1所示。此外，水下航行器的材料参数包括：弹性模量2×10¹¹ Pa、泊松比0.3、密度7 800 kg/m³、海水密度1 025 kg/m³、水中声速1 500 m/s。

图  4  设计变量示意图

Figure  4.  Diagram of the design variables

下载: 全尺寸图片 幻灯片

表  1  设计变量的取值范围

Table  1.  Ranges of the design variables

设计变量	取值范围
平行舯体直径$\phi$/m	[2.286, 2.794]
平行舯体长$L$/m	[9.144, 11.176]
艏径$R1$/m	[4.068, 4.972]
艉径$R2$/m	[8.478, 10.362]

下载: 导出CSV 

| 显示表格

本文使用了ANSYS的Modal模块来实现水下航行器的振动湿模态分析。首先，导入SOLIDWORKS中创建的参数化几何模型，然后，在水下航行器四周建立长方体形的流体域作为外界流场，其尺寸为水下航行器各端缘向外扩展3 m， 并双选水下航行器和流体域建立一个整体件，以便后续划分网格时二者公共面处的节点位置重合。水下航行器的有限元网格采取自动划分方式，得到的网格尺寸为0.25 m，网格数37 388。流体域的有限元网格划分方式为空间六面体占优，网格尺寸为0.5 m，网格数30 035。本文算例采用了单向流固耦合分析方法。分析时，流体域与水下航行器的公共面为流固耦合面，定义流体域外表面的边界压力为0，插入APDL命令流，调用非对称模态算法求解结构动力学有限元方程。由于每次仿真分析耗时约6 min，直接调用仿真程序进行优化的效率低下，所以本文算例符合基于缩减型梯度增强Kriging模型的优化方法应用场景。

以下为优化设计流程：首先，设定初始采样成本为20次仿真调用。采用拉丁超立方采样确定$4 \times 4$的初始采样矩阵${\boldsymbol{X}}$，调用仿真得到4组响应值，并通过有限差分法计算采样位置的梯度。需注意，此时初始采样样本数据向量中共包含20组响应和4组梯度。 然后，运用基于缩减型梯度增强Kriging模型的优化方法进行序列优化，终止条件为加点过程中调用的仿真次数超过40次。对于本文算例，取$Q{\rm{ = 10}}$，${\varepsilon _l}$($l = 1, \cdots ,4$)等于样本所对应维度偏导绝对值中最大值的千分之一。最后，作为参照，在相同样本量的前提下，运用基于标准型梯度增强Kriging模型的EGO算法和经典的EGO算法进行序列优化，各方法优化过程的迭代曲线如图5所示。

图  5  水下航行器第7阶固有频率优化过程的迭代曲线

Figure  5.  Iterative curves of the optimization process of seventh-order natural frequency for the underwater vehicle

下载: 全尺寸图片 幻灯片

由图5可知，与基准形参考值相比，3种优化方法在迭代终止时得到的设计方案一定程度上均提高了水下航行器的第7阶固有频率，而基于缩减型梯度增强Kriging模型的优化方法在效率方面最高。图6所示为水下航行器优化前后的仿真分析结果。优化后的水下航行器第7阶固有频率提高了14.6%，从而验证了所提出的优化方法在船舶优化领域的可行性。

图  6  水下航行器结构优化前后对比

Figure  6.  Comparison of the underwater vehicle structural design before and after optimization

下载: 全尺寸图片 幻灯片

4.   结　语

本文针对船舶优化领域中梯度信息只能通过有限差分法获取的问题，提出了基于缩减型梯度增强Kriging代理模型的优化方法，扩展了标准型梯度增强Kriging方法的应用范围，并以某水下航行器结构优化为例对所提方法进行了验证。结果表明，该方法在船舶结构优化领域具有良好的应用前景。在基于梯度增强Kriging方法的船舶结构优化领域，未来值得进一步研究的相关问题包括：

1） 探索合理的梯度/偏导数筛选方法，以缓解船舶结构优化中因设计变量的增多而导致计算量剧增的问题；同时，致力于开发船舶结构有限元分析所对应的高效梯度计算（伴随方法）程序。

2） 探索基于梯度增强Kriging方法的船舶结构优化过程中约束（例如强度、寿命等）昂贵、优化目标（例如重量、模态等）多的解决方法。

3） 探索基于梯度增强Kriging方法的船舶结构优化过程中涉及离散（例如型材数量、规格等）/混合变量时的解决方法。