Processing math: 0%

基于梯度增强Kriging方法的水下航行器结构优化设计

陈力铭, 邱浩波, 高亮

陈力铭, 邱浩波, 高亮. 基于梯度增强Kriging方法的水下航行器结构优化设计[J]. 中国舰船研究, 2021, 16(4): 79–85. DOI: 10.19693/j.issn.1673-3185.02066
引用本文: 陈力铭, 邱浩波, 高亮. 基于梯度增强Kriging方法的水下航行器结构优化设计[J]. 中国舰船研究, 2021, 16(4): 79–85. DOI: 10.19693/j.issn.1673-3185.02066
CHEN L M, QIU H B, GAO L. Structural design optimization of underwater vehicle via Gradient-enhanced Kriging[J]. Chinese Journal of Ship Research, 2021, 16(4): 79–85. DOI: 10.19693/j.issn.1673-3185.02066
Citation: CHEN L M, QIU H B, GAO L. Structural design optimization of underwater vehicle via Gradient-enhanced Kriging[J]. Chinese Journal of Ship Research, 2021, 16(4): 79–85. DOI: 10.19693/j.issn.1673-3185.02066
陈力铭, 邱浩波, 高亮. 基于梯度增强Kriging方法的水下航行器结构优化设计[J]. 中国舰船研究, 2021, 16(4): 79–85. CSTR: 32390.14.j.issn.1673-3185.02066
引用本文: 陈力铭, 邱浩波, 高亮. 基于梯度增强Kriging方法的水下航行器结构优化设计[J]. 中国舰船研究, 2021, 16(4): 79–85. CSTR: 32390.14.j.issn.1673-3185.02066
CHEN L M, QIU H B, GAO L. Structural design optimization of underwater vehicle via Gradient-enhanced Kriging[J]. Chinese Journal of Ship Research, 2021, 16(4): 79–85. CSTR: 32390.14.j.issn.1673-3185.02066
Citation: CHEN L M, QIU H B, GAO L. Structural design optimization of underwater vehicle via Gradient-enhanced Kriging[J]. Chinese Journal of Ship Research, 2021, 16(4): 79–85. CSTR: 32390.14.j.issn.1673-3185.02066

基于梯度增强Kriging方法的水下航行器结构优化设计

基金项目: 国家自然科学基金资助项目(51675198);国防基础科研计划资助项目(41423010205);华中科技大学学术前沿青年团队资助项目(2017QYTD04)
详细信息
    作者简介:

    陈力铭,男,1995年生,博士生。研究方向:基于代理模型的设计优化。E-mail:liming_chen@hust.edu.cn

    邱浩波,男,1975年生,博士,教授,博士生导师。研究方向:基于代理模型的设计优化,可靠性设计优化,工艺智能优化方法,设备健康管理和系统可靠性建模。E-mail:hobbyqiu@hust.edu.cn

    高亮,男,1974年生,博士,教授,博士生导师。研究方向:智能优化方法及其在设计制造中的应用。E-mail:gaoliang@mail.hust.edu.cn

    通讯作者:

    邱浩波

  • 中图分类号: U662.2

Structural design optimization of underwater vehicle via Gradient-enhanced Kriging

知识共享许可协议
基于梯度增强Kriging方法的水下航行器结构优化设计陈力铭,采用知识共享署名4.0国际许可协议进行许可。
  • 摘要:
      目的  船舶结构优化设计过程通常涉及对高精度数值仿真进行响应分析,其耗时特性决定了可调用的仿真次数十分有限,使得优化过程受到限制。为了探索基于梯度增强Kriging代理模型的高效设计优化方法,缩短设计周期,节省设计成本,提出基于缩减型梯度增强Kriging的加点策略,仅在有实际改进的采样位置进行梯度计算,以减少仿真调用次数。
      方法  首先,使用多起点局部优化算法搜索改进期望函数的若干局部最优解作为候选加点位置;然后,计算相应的近似驻点概率,并根据改进期望值和近似驻点概率值的一致程度来确定加点位置,从而提高优化效率;最后,针对某水下航行器结构进行优化设计,以提高其水下无约束自由振动时的第7阶固有频率为目标,对所提方法的可行性进行验证。
      结果  结果表明,优化后的固有频率值与基准型相比提升了14.6%,方法的可行性得到验证。
      结论  所提方法可以将基于梯度增强Kriging代理模型的优化方法泛化至梯度信息只能通过有限差分法获取的场景。
    Abstract:
      Objectives  The structural optimization of ships usually involves the use of high-fidelity numerical simulations which are time-consuming and thus difficult to evaluated frequently, and this intrinsic property hinders the optimization process. To promote efficient design optimization, this paper explores the use of Gradient-enhanced Kriging (GEK) surrogate mode in order to shorten the design loop and save design cost. A reduced GEK-based infill criterion is proposed to decrease the number of simulations by calculating the gradients only for sample locations where improvement occurs.
      Methods  A multi-start local optimization algorithm is employed to search the local optima of the "expected improvement" function and locate candidate infill points. The associated "approximate probability of stationary point (APSP)" values are also evaluated, and infill decisions are made according to the extent of consistency between these two quantities, thereby improving optimization efficiency. The proposed method is then applied to the structural optimization of an underwater vehicle to increase the seventh-order natural frequency under unconstrained free vibration in an underwater environment, and the validity is verifed.
      Results  The result shows that, compared with the baseline, the optimized design achieves a 14.6% improvement.
      Conclusions  The proposed GEK-based optimization method can be generalized to cases when gradients can only be evaluated by finite difference.
  • 随着船舶结构复杂程度的不断提升,高精度的有限元仿真软件已成为船舶结构分析的主流工具,而直接调用耗时的仿真分析程序进行优化已难以满足船舶设计快速高效的需求,因此,发展基于代理模型的高效优化方法是船舶结构优化设计的重点研究领域之一[1]。所谓代理模型,是指对一系列通过仿真分析获得的输入/输出样本进行拟合得到的统计近似模型。在代理模型提供的信息引导下的优化可以极大地减少耗时的仿真次数,达到缩短设计周期、节省设计成本的目的。目前,经典的代理模型建模与优化方法已被广泛用于船舶领域,例如对船舶板架强度和稳定性[2]、潜艇振动声辐射[3]、水下结构物基座阻抗特性[4]等响应的快速预报,以及对深潜器的多球交接耐压壳[5]、多用途船货舱段[6]的结构优化和对海上发射船[7]、集装箱船[8] 等的船型优化。

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

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

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

    Y\left( {\boldsymbol{x}} \right) = {f^{\rm{T}}}\left( {\boldsymbol{x}} \right){\boldsymbol{\beta }} + Z\left( {\boldsymbol{x}} \right), (1)

    式中,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}}}中既包含了响应也包含了梯度,记为

    \begin{split} & {\tilde{\boldsymbol{y}}} = \left[ y\left( {{{\boldsymbol{x}}^{\left( 1 \right)}}} \right), \cdots ,y\left( {{{\boldsymbol{x}}^{\left( n \right)}}} \right),\frac{{\partial y}}{{\partial {x_1}}}\left( {{{\boldsymbol{x}}^{\left( 1 \right)}}} \right), \cdots ,\right.\\&\qquad \left.\frac{{\partial y}}{{\partial {x_1}}}\left( {{{\boldsymbol{x}}^{\left( n \right)}}} \right), \cdots ,\frac{{\partial y}}{{\partial {x_k}}}\left( {{{\boldsymbol{x}}^{\left( n \right)}}} \right) \right]^{\rm{T}} \end{split} (2)

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

    \begin{split} & L\left( {{\boldsymbol{\beta }},{\sigma ^2},{\boldsymbol{\theta }}\left| {{\tilde{\boldsymbol{y}}}} \right.} \right) = \frac{1}{{{{\left( {2{\text{π}} \cdot {\sigma ^2}} \right)}^{\frac{{n + nk}}{2}}}{{| {{\tilde{\boldsymbol{R}}}} |}^{\frac{1}{2}}}}}\\&\exp \left[ { - \frac{{{{\left( {{\tilde{\boldsymbol{y}}} - {\tilde{\boldsymbol{ F}}}{\boldsymbol{\beta }}} \right)}^{\rm{T}}}{{{\tilde{\boldsymbol{R}}}}^{ - 1}}\left( {{\tilde{\boldsymbol{y}}} - {\tilde{\boldsymbol{ F}}}{\boldsymbol{\beta}}} \right)}}{{2{\sigma ^2}}}} \right] \end{split} (3)

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

    \tilde{\boldsymbol{F}}\!\!=\!\!{\left[\!\!{f}^{\rm{T}}({\boldsymbol{x}}^{(1)}),\!\cdots \!,{f}^{\rm{T}}({\boldsymbol{x}}^{(n)}),\partial {f}^{\rm{T}}\!/\!\partial {x}_{1}({\boldsymbol{x}}^{(1)}),\!\cdots\! ,\partial {f}^{\rm{T}}\!/\!\partial {x}_{k}({\boldsymbol{x}}^{n})\!\right]}^{\rm{T}} (4)
    {\tilde{\boldsymbol{R}}} = \left[ {\begin{array}{*{20}{c|c}} {\boldsymbol{R}}& {\dfrac{{\partial {\boldsymbol{R}}}}{{\partial {\boldsymbol{x}}}}} \\ \hline \\[-10pt] {{{\left( {\dfrac{{\partial {\boldsymbol{R}}}}{{\partial {\boldsymbol{x}}}}} \right)}^{\rm{T}}}}& {\dfrac{{{\partial ^2}{\boldsymbol{R}}}}{{\partial {{\boldsymbol{x}}^2}}}} \end{array}} \right] (5)

    作为\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,可得到极大似然估计量分别如下:

    {\hat{\boldsymbol{\beta }}} = {\left( {{{\tilde {\boldsymbol{F}}}^{\rm{T}}}{{\tilde {\boldsymbol{R}}}^{ - 1}}{\tilde {\boldsymbol{F}}}} \right)^{ - 1}}{{\tilde {\boldsymbol{F}}}^{\rm{T}}}{{\tilde {\boldsymbol{R}}}^{ - 1}}{\tilde {\boldsymbol{y}}} (6)
    {\hat \sigma ^2} = \frac{{{{\left( {{\tilde {\boldsymbol{y}}} - {{\tilde{\boldsymbol{F}}{\hat {\boldsymbol{\beta}} }}}} \right)}^{\rm{T}}}{{\tilde {\boldsymbol{R}}}^{ - 1}}\left( {{\tilde {\boldsymbol{y}}} - {{\tilde{\boldsymbol{F}}{\hat {\boldsymbol{\beta}} }}}} \right)}}{{n + nk}} (7)

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

    \hat y\left( {\boldsymbol{x}} \right){\rm{ = }}{f^{\rm{T}}}\left( {\boldsymbol{x}} \right){\hat{\boldsymbol{\beta }}} + {{\tilde{\boldsymbol{r}}}^{\rm{T}}}{{\tilde{\boldsymbol{R}}}^{ - 1}}\left( {{\tilde{\boldsymbol{y}}} - {\tilde{\boldsymbol{F}}{\hat{\boldsymbol{\beta }}}}} \right) (8)
    {s^2}\left( {\boldsymbol{x}} \right){\rm{ = }}{\hat \sigma ^2}\left[ {1 - {{{\tilde{\boldsymbol{r}}}}^{\rm{T}}}{{{\tilde{\boldsymbol{R}}}}^{ - 1}}{\tilde{\boldsymbol{r}}} + {{{\tilde{\boldsymbol{ u}}}}^{\rm{T}}}{{\left( {{{{\tilde{\boldsymbol{F}}}}^{\rm{T}}}{{{\tilde{\boldsymbol{R}}}}^{ - 1}}{\tilde{\boldsymbol{F}}}} \right)}^{ - 1}}{\tilde{\boldsymbol{ u}}}} \right] (9)

    式中,{\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]。

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

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

    {\boldsymbol{X}}' = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{X}}_0}} \\ {{{\boldsymbol{X}}_1}} \end{array}} \right] (10)

    式中,{{\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}}'}写为

    \begin{split} & {\tilde{\boldsymbol{y}}'} = \left[ y\left( {{{\boldsymbol{x}}^{\left( 1 \right)}}} \right), \cdots ,y\left( {{{\boldsymbol{x}}^{\left( n \right)}}} \right),\frac{{\partial y}}{{\partial {x_1}}}\left( {{{\boldsymbol{x}}^{\left( {{n_0} + 1} \right)}}} \right), \cdots ,\right.\\&\qquad \left.\frac{{\partial y}}{{\partial {x_1}}}\left( {{{\boldsymbol{x}}^{\left( n \right)}}} \right), \cdots ,\frac{{\partial y}}{{\partial {x_k}}}\left( {{{\boldsymbol{x}}^{\left( n \right)}}} \right) \right]^{\rm{T}} \end{split} (11)

    构建缩减型梯度增强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)}},有

    {{\boldsymbol{v}}^{\rm{T}}}\left( {{{\hat {\sigma} }'^{2}}{\tilde{\boldsymbol{R}}'}} \right){\boldsymbol{v}} = Var\left[ {{{\boldsymbol{v}}^{\rm{T}}}{\tilde{\boldsymbol{ Y}}'}} \right] > 0 (12)

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

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

    \frac{{\partial Y}}{{\partial {\boldsymbol{x}}}}\left( {\boldsymbol{x}} \right) = \frac{{\partial {f^{\rm{T}}}}}{{\partial {\boldsymbol{x}}}}\left( {\boldsymbol{x}} \right){\boldsymbol{\beta }} + \frac{{\partial Z}}{{\partial {\boldsymbol{x}}}}\left( {\boldsymbol{x}} \right) (13)

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

    \frac{{\partial \hat y}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right) = {\tilde{\boldsymbol{ c}}}_l^{\rm{T}}\left( {\boldsymbol{x}} \right){\tilde{\boldsymbol{y}}} (14)

    式中,{{\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)最小。

    MSE\left[ {\frac{{\partial \hat y}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right)} \right] = E\left[ {{{\left( {{\tilde{\boldsymbol{ c}}}_l^{\rm{T}}\left( {\boldsymbol{x}} \right){\tilde{\boldsymbol{Y}}} - \frac{{\partial Y}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right)} \right)}^2}} \right] (15)

    同时,满足如下无偏约束:

    E\left[ {{\tilde{\boldsymbol{ c}}}_l^{\rm{T}}\left( {\boldsymbol{x}} \right){\tilde{\boldsymbol{Y}}}} \right] = E\left[ {\frac{{\partial Y}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right)} \right] (16)

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

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

    \frac{{\partial \hat y}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right){\rm{ = }}\frac{{\partial {f^{\rm{T}}}}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right){\hat{\boldsymbol{\beta }}} + {\tilde{\boldsymbol{r}}}_l^{\rm{T}}{{\tilde{\boldsymbol{R}}}^{ - 1}}\left( {{\tilde{\boldsymbol{y}}} - {\tilde{\boldsymbol{F}}}{\hat{\boldsymbol{\beta }}}} \right) (17)
    \begin{split} & \frac{{\partial {s^2}}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right){\rm{ = }}{\hat \sigma ^2}\left[ Corr\left[ {\frac{{\partial Z}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right),\frac{{\partial Z}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right)} \right] -\right.\\&\qquad\left. {\tilde{\boldsymbol{r}}}_l^{\rm{T}}{{{\tilde{\boldsymbol{R}}}}^{ - 1}}{{{\tilde{\boldsymbol{r}}}}_l} + {\boldsymbol{u}}_l^{\rm{T}}{{\left( {{{{\tilde{\boldsymbol{F}}}}^{\rm{T}}}{{{\tilde{\boldsymbol{R}}}}^{ - 1}}{\tilde{\boldsymbol{F}}}} \right)}^{ - 1}}{{\boldsymbol{u}}_l} \right] \end{split} (18)

    式中, {{\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的概率模型属性,预测的梯度服从正态分布

    \frac{{\partial Y}}{{\partial {\boldsymbol{x}}}}\left( {\boldsymbol{x}} \right) \sim N\left( {\frac{{\partial \hat y}}{{\partial {\boldsymbol{x}}}}\left( {\boldsymbol{x}} \right),\frac{{\partial {s^2}}}{{\partial {\boldsymbol{x}}}}\left( {\boldsymbol{x}} \right)} \right) (19)

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

    {P_l}\left( { - {\varepsilon _l} \leqslant \frac{{\partial Y}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right) \leqslant + {\varepsilon _l}} \right) = \int_{ - {\varepsilon _l}}^{ + {\varepsilon _l}} {f\left( {\frac{{\partial y}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right)} \right)} {\rm{d}}\frac{{\partial y}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right) (20)

    式中,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)可定义为

    APSP\left( {\boldsymbol{x}} \right) = \prod\limits_{l = 1}^k {{P_l}\left( { - {\varepsilon _l} \leqslant \frac{{\partial Y}}{{\partial {x_l}}}\left( {\boldsymbol{x}} \right) \leqslant + {\varepsilon _l}} \right)} (21)

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

    APSP\left({\boldsymbol{x}}^{\rm{*}}\right)=\left\{ \begin{array}{l}1,\;\;{\text{若}}\partial \hat{y}/\partial {\boldsymbol{x}}\left({\boldsymbol{x}}^{*}\right)=\vec 0\\ 0,\;\;{\text{其他}}\end{array} \right. (22)

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

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

    f\left( x \right) = {\left( {6x - 2} \right)^2}\sin\left( {12x - 4} \right),\;x \in \left[ {0,1} \right] (23)

    通过拉丁超立方采样[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),计算式为

    \begin{split} & {\varDelta _q} = \frac{{EI\left( {H_{\rm{sorted}}^{\left( 1 \right)}} \right) - EI\left( {H_{\rm{sorted}}^{\left( q \right)}} \right)}}{{EI\left( {H_{\rm{sorted}}^{\left( 1 \right)}} \right)}} -\\& \frac{{APSP\left( {H_{\rm{sorted}}^{\left( q \right)}} \right) - APSP\left( {H_{\rm{sorted}}^{\left( 1 \right)}} \right)}}{{APSP\left( {H_{\rm{sorted}}^{\left( 1 \right)}} \right)}} \end{split} (24)

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

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

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

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

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

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

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

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

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

    Figure  1.   Illustrations for the probability {P_l}

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

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

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

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

    图  4   设计变量示意图

    Figure  4.   Diagram of the design variables

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

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

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

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

    表  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
  • [1] 赵留平, 詹大为, 程远胜, 等. 船舶结构优化设计技术研究进展[J]. 中国舰船研究, 2014, 9(4): 1–10. doi: 10.3969/j.issn.1673-3185.2014.04.001

    ZHAO L P, ZHAN D W, CHENG Y S, et al. Review on optimum design methods of ship structures[J]. Chinese Journal of Ship Research, 2014, 9(4): 1–10 (in Chinese). doi: 10.3969/j.issn.1673-3185.2014.04.001

    [2] 郑少平, 陈静, 程远胜, 等. 代理模型技术及其在船舶板架强度和稳定性计算中的应用[J]. 中国造船, 2013, 54(1): 40–51. doi: 10.3969/j.issn.1000-4882.2013.01.007

    ZHENG S P, CHEN J, CHENG Y S, et al. Surrogate models and their application in calculation of strength and stability of ship grillage[J]. Shipbuilding of China, 2013, 54(1): 40–51 (in Chinese). doi: 10.3969/j.issn.1000-4882.2013.01.007

    [3] 郭明慧, 黎胜. 运用代理模型方法预测潜艇结构模型的振动声辐射[J]. 中国舰船研究, 2013, 8(6): 69–74. doi: 10.3969/j.issn.1673-3185.2013.06.012

    GUO M H, LI S. A surrogate model for structural vibration and acoustic radiation of underwater submarine structures[J]. Chinese Journal of Ship Research, 2013, 8(6): 69–74 (in Chinese). doi: 10.3969/j.issn.1673-3185.2013.06.012

    [4] 夏志, 刘均, 程远胜. 基于代理模型的水下结构物基座阻抗特性快速预报[J]. 中国舰船研究, 2020, 15(3): 81–87. doi: 10.19693/j.issn.1673-3185.01536

    XIA Z, LIU J, CHENG Y S. Fast prediction of mechanical impedance of an underwater foundation based on surrogate models[J]. Chinese Journal of Ship Research, 2020, 15(3): 81–87 (in Chinese). doi: 10.19693/j.issn.1673-3185.01536

    [5] 苟鹏, 崔维成. 基于Kriging模型的深潜器多球交接耐压壳结构优化[J]. 船舶力学, 2009, 13(1): 100–106. doi: 10.3969/j.issn.1007-7294.2009.01.013

    GOU P, CUI W C. Structural optimization of multiple intersecting spherical pressure hulls based on Kriging model[J]. Journal of Ship Mechanics, 2009, 13(1): 100–106 (in Chinese). doi: 10.3969/j.issn.1007-7294.2009.01.013

    [6] 袁野, 王德禹, 李喆. 基于支持向量机的船舶结构优化方法[J]. 舰船科学技术, 2013, 35(7): 12–17, 31. doi: 10.3404/j.issn.1672-7649.2013.07.003

    YUAN Y, WANG D Y, LI Z. Optimization of ship structure based on support vector machine[J]. Ship Science and Technology, 2013, 35(7): 12–17, 31 (in Chinese). doi: 10.3404/j.issn.1672-7649.2013.07.003

    [7] 郝浩浩, 韩端锋, 高良田, 等. 基于近似模型的海上发射船船型优化[J]. 船舶工程, 2017, 39(7): 1–5, 54.

    HAO H H, HAN D F, GAO L T, et al. Hull form optimization of the sea-launch catamarans based on approximation model[J]. Ship Engineering, 2017, 39(7): 1–5, 54 (in Chinese).

    [8] 常海超. 近似理论在船型优化中的应用研究[D]. 武汉: 武汉理工大学, 2014.
    [9]

    LAURENT L, LE RICHE R, SOULIER B, et al. An overview of gradient-enhanced metamodels with application[J]. Archives of Computational Methods in Engineering, 2019, 26(1): 61–106. doi: 10.1007/s11831-017-9226-3

    [10]

    MORRIS M D, MITCHELL T J, YLVISAKER D. Bayesian design and analysis of computer experiments: use of derivatives in surface prediction[J]. Technometrics, 1993, 35(3): 243–255. doi: 10.1080/00401706.1993.10485320

    [11] 韩忠华, 许晨舟, 乔建领, 等. 基于代理模型的高效全局气动优化设计方法研究进展[J]. 航空学报, 2020, 41(5): 623344.

    HAN Z H, XU C Z, QIAO J L, et al. Recent progress of efficient global aerodynamic shape optimization using surrogate-based approach[J]. Acta Aeronautica et Astronautica Sinica, 2020, 41(5): 623344 (in Chinese).

    [12]

    LEARY S J, BHASKAR A, KEANE A J. A derivative based surrogate model for approximating and optimizing the output of an expensive computer simulation[J]. Journal of Global Optimization, 2004, 30(1): 39–58. doi: 10.1023/B:JOGO.0000049094.73665.7e

    [13]

    LAURENT L, BOUCARD P A, SOULIER B. Generation of a cokriging metamodel using a multiparametric strategy[J]. Computational Mechanics, 2013, 51(2): 151–169. doi: 10.1007/s00466-012-0711-0

    [14]

    LOCKWOOD B A, ANITESCU M. Gradient-enhanced universal Kriging for uncertainty propagation[J]. Nuclear Science and Engineering, 2012, 170(2): 168–195. doi: 10.13182/NSE10-86

    [15]

    GILES M B, PIERCE N A. An introduction to the adjoint approach to design[J]. Flow, Turbulence and Combustion, 2000, 65(3): 393–415.

    [16]

    CHEN L M, QIU H B, GAO L, et al. Optimization of expensive black-box problems via Gradient-enhanced Kriging[J]. Computer Methods in Applied Mechanics and Engineering, 2020, 362: 112861. doi: 10.1016/j.cma.2020.112861

    [17]

    HAN Z H, ZHANG Y, SONG C X, et al. Weighted gradient-enhanced Kriging for high-dimensional surrogate modeling and design optimization[J]. AIAA Journal, 2017, 55(12): 4330–4346. doi: 10.2514/1.J055842

    [18]

    CHEN L M, QIU H B, GAO L, et al. A screening-based gradient-enhanced Kriging modeling method for high-dimensional problems[J]. Applied Mathematical Modelling, 2019, 69: 15–31. doi: 10.1016/j.apm.2018.11.048

    [19]

    FORRESTER A I J, KEANE A J. Recent advances in surrogate-based optimization[J]. Progress in Aerospace Sciences, 2009, 45(1–3): 50–79. doi: 10.1016/j.paerosci.2008.11.001

    [20]

    JONES D R, SCHONLAU M, WELCH W J. Efficient global optimization of expensive black-box functions[J]. Journal of Global Optimization, 1998, 13(4): 455–492. doi: 10.1023/A:1008306431147

    [21]

    MCKAY M D, BECKMAN R J, CONOVER W J. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J]. Technometrics, 1979, 21(2): 239–245.

    [22]

    ZHAN D W, QIAN J C, CHENG Y S. Balancing global and local search in parallel efficient global optimization algorithms[J]. Journal of Global Optimization, 2017, 67(4): 873–892. doi: 10.1007/s10898-016-0449-x

    [23]

    BYRD R H, HRIBAR M E, NOCEDAL J. An interior point algorithm for large-scale nonlinear programming[J]. SIAM Journal on Optimization, 1999, 9(4): 877–900. doi: 10.1137/S1052623497325107

  • 期刊类型引用(2)

    1. 孔德冬,夏茂龙,李永正,张海华,刘双,张宇阳. 声呐透声窗结构的减振降噪优化设计. 中国造船. 2024(02): 63-73 . 百度学术
    2. 龚俊斌,王鹏九,汪皓,杨萌. 基于智能模糊推理的UUV艇型参数生成方法研究. 中国舰船研究. 2024(06): 56-63 . 本站查看

    其他类型引用(8)

  • 其他相关附件

图(6)  /  表(1)
计量
  • 文章访问数:  754
  • HTML全文浏览量:  260
  • PDF下载量:  122
  • 被引次数: 10
出版历程
  • 收稿日期:  2020-08-10
  • 修回日期:  2020-10-27
  • 网络出版日期:  2021-03-19
  • 刊出日期:  2021-08-09

目录

/

返回文章
返回