基于兴趣子域动态代理模型的船舶结构可靠性优化

罗文俊，王德禹

1 中远海运特种运输股份有限公司，广东 广州 510623

2 上海交通大学 海洋工程国家重点实验室，上海 200240

3 高新船舶与深海开发装备协同创新中心，上海 200240

0 引 言

基于可靠性的优化设计因在满足船舶结构经济性的同时能够保证足够的安全性而极具优势。目前，可靠性优化设计已经在航空航天、汽车和船舶等领域得到广泛应用[1-3]。

可靠性优化设计方法主要分为3 种：双循环法、单循环法和解耦法。双循环策略虽然精度高，但计算成本较大、效率较低，而后两种方法能有效解决此问题。Liang 等[4]根据库恩-塔克（KKT）最优条件代替可靠性分析过程，通过计算得到近似最小功能点，将概率约束转化为确定性约束，提出了单循环法。刘勤等[5]采用单循环法将优化迭代计算与可靠寿命的迭代求解同步进行，并通过传动箱的结构可靠性优化结果验证了方法的高效性。单循环法虽然效率较高，但只适用于线性和弱非线性问题，对于强非线性问题，无法保证精度甚至有可能出现迭代不收敛的情况。而解耦策略在保证精度的前提下，不仅能提高求解效率，同时还具有良好的适用性。Du 和Chen[6]提出的序贯优化与可靠性评估（sequential optimization and reliability assessment, SORA）法通过将可靠性分析结果转化为优化边界条件的偏移向量，使确定性优化与可靠性分析被解耦。李海燕和井元伟[7]将SORA 法应用于多学科可靠性优化，降低了协同优化计算的复杂度。解耦方法较前2 种集成策略具有很大的优势，而SORA 法是解耦方法中应用最为广泛的一种方法。

船舶结构响应属于多参数、强非线性问题，采用代理模型替代有限元模型能极大地提高可靠性优化效率，对于强非线性和高维度工程问题，Kriging 模型能实现良好的拟合。刘瞻等[8]通过结合重要抽样法和人工蜂群算法，建立了参数优化的Kriging 模型，有效提高了可靠度计算精度和效率。张干锋和王德禹[9]对常规Kriging 模型进行了分区间泛化和参数融合，其将改进后的Kriging模型用于船舶结构型线优化，在预测多维度响应时具有更高的精度。代理模型分静态代理模型和动态代理模型2 种。静态代理模型需要在设计空间中采集足够的样本点才能构造高拟合精度的代理模型，计算成本较大，而动态代理模型采用序贯抽样方法，能以较少的样本点实现对有限元模型的高度拟合。Echard 等[10]引入了一种学习函数，其以该学习函数作为筛选准则进行序贯抽样来构造动态Kriging 模型，以较少的样本点实现了对极限状态边界的高度拟合。高月华和王希诚[11]基于Kriging 代理模型提出了多点加点准则，基于样本集信息和预测函数特征添加新样本集，在寻优迭代之时自适应地提升代理模型的精度。

根据结构可靠性理论，最大可能失效点（MPP）周围区域对可靠度的贡献最大，而结构优化的最优解通常位于约束边界附近，因此可靠性优化必然存在一个可能存在最优解的区域，即兴趣子域。Zhao 等[12]通过计算影响可靠度主要区域的大小，提出兴趣子域的概念，从而构造了在MPP 附近区域实现高度近似的代理模型。龙周等[13]引入少数类合成的过采样算法(synthetic minority oversampling technique，SMOTE)对MPP 附近进行过采样，构建对极限状态边界高度拟合的BP 神经网络模型，极大地提高了计算效率和精度。

针对船舶结构因强非线性、多参数、多响应而导致计算成本较大、难以获取大量样本点的问题，本文拟提出可靠性优化策略。首先，提出并确定基于SORA 法中圆弧搜索法（arc search method，ASM）的兴趣子域范围，基于信息熵函数H制定自适应空间减缩规则，对设计空间进行不断的缩减，进而构造基于兴趣子域的自适应空间减缩序贯抽样策略，以采用尽可能少的样本点构造对兴趣子域进行高度拟合的Kriging动态代理模型，减少计算成本；其次，提出基于SORA 法的概率约束可行性检查方法，减少不必要的可靠性评估过程，提高可靠性分析效率；最后，将Kriging 模型和多岛遗传算法（multi-island genetic algorithm，MIGA）嵌入SORA 法中进行可靠性优化，以保证可靠性优化结果的准确性。

1 基于兴趣子域的自适应空间减缩及序贯抽样策略

1.1 圆弧搜索法

本文中的SORA 法采用ASM 法作为可靠性分析方法，该方法由Du 等[14]提出，用于将逆可靠性分析法中的优化过程转化为数学迭代计算，从而简化计算，提高可靠性分析效率。ASM 法的原理如下：逆可靠性分析法是在期望可靠度指标的等值面（||u||= βa）上寻找概率约束性能函数G(u)（其中，G＞0 为可行性条件）的最小值，约束条件为可靠度指标期望值 βa。根据KKT 条件，可以将逆可靠性分析的优化过程推导为如下迭代公式：

ASM 法的迭代原理如图1 所示。图中：x1，x2为设计变量；Gk为概率约束性能函数第k次迭代取值。

图1 圆弧搜索法[15]Fig. 1 Arc search method[15]

当式（1）迭代收敛至uk+1，即功能函数与期望可靠度指标面相切时，若G(uk+1)＜G(uk)，说明样本点朝着功能函数G减小的方向前进，uk+1即为可靠性分析的逆最大可能失效点；当出现G(uk+1)≥G(uk)的情况时，进行第2 阶段寻优，寻优空间为期望可靠度指标面上uk+1和uk之间的圆弧。采用以下数学优化模型进行寻优：

通过第2 阶段的寻优迭代策略，可以保证新的迭代收敛点的功能函数值比上一个最优点的更小，实现对原迭代策略的补充完善。

本文充分利用ASM 法的寻优过程原理，提出基于兴趣子域的自适应空间减缩及序贯抽样策略，实现了对兴趣子域的高度拟合；而对于非兴趣子域，则采用空间减缩技术，无需对这些区域进行精确拟合，因而减少了可靠性优化的循环次数，可在保证拟合精度的同时提高效率。

1.2 兴趣子域范围的确定

ASM 法是在半径为可靠度指标β 的圆（球面或超球面）上寻找功能函数的最小值，如图2 所示，分别为当功能函数为凹函数和凸函数时的搜索过程（G＞0 为可行域）。

图2 圆弧搜索法寻优过程[15]Fig. 2 Optimization process of arc search method[15]

1.3 自适应空间减缩及序贯抽样策略

1.3.1 自适应空间减缩规则

图3 兴趣子域示意图Fig. 3 Schematic diagram of interest subdomain

构建对兴趣子域进行高度拟合的动态代理模型，可同时提高确定性优化和可靠性分析过程的计算效率。本节将采用蒙特卡罗抽样生成整个设计空间的样本点集S，采用Lv 等[16]提出的AK-LS方法中基于信息熵的主动学习函数H实现序贯抽样，在提高Kriging 模型全局拟合精度的同时，采用自适应的空间减缩技术不断删减样本点集S中位于非兴趣子域的样本点，从而缩小设计空间的范围直至兴趣子域。学习函数H如下：

由于初始Kriging 模型的全局拟合精度较低，因此非兴趣子域中样本点的删减条件应当远离兴趣子域。随着Kriging 模型拟合精度的不断提高，非兴趣子域中样本点的删减条件自适应地不断逼近兴趣子域边界，即在进行序贯抽样更新动态代理模型的同时不断减缩设计空间。最后，只保留兴趣子域及其附近的部分样本点，通过学习函数H构造代理模型来实现对兴趣子域的局部高精度拟合。考虑到随着样本点的不断增加，整个设计空间中样本点的Hmax值逐渐减小，因此，本文在序贯抽样的同时采用Hmax作为筛选因子。

本文采用的自适应空间减缩规则如下：

其中，筛选条件中下确界的设定参考了Cox和John[17]提出的用于优化的下置信界（lcb）函数，下置信界函数如下所示：

1.3.2 自适应空间减缩及序贯抽样流程

基于兴趣子域的自适应空间减缩及序贯抽样策略是通过不断减缩全局设计空间，直至保留兴趣子域，同时通过具备高效提高全局精度特点的学习函数H建立Kriging 动态代理模型来拟合兴趣子域。

Lv 等[16]通过大量的算例结果分析指出，当H函数对应的迭代停止准则为Hmax≤0.5 时，代理模型具有足够的拟合精度。因此，为保证Kriging模型在兴趣子域内具有足够的拟合精度，本文拟采用Hmax≤0.5 作为自适应空间减缩和序贯抽样的迭代停止准则。

基于兴趣子域的自适应空间减缩及序贯抽样策略流程如图4 所示。

图4 自适应空间减缩及序贯抽样策略Fig. 4 Adaptive spatial reduction and sequence sampling strategy

1） 在设计空间中采用拉丁超立方抽样或蒙特卡罗抽样生成样本集S，S集中的所有样本点均无需调用有限元计算。

2） 生成训练集T。采用最优拉丁超立方抽样生成少量样本点，作为训练集T。

3） 调用功能函数或有限元软件计算，得到训练集响应值G。

4） 根据训练集T，同时对约束条件和目标函数使用Matlab 自带的DACE 工具箱建立多个Kriging 模型。

5） 分别用Kriging 模型预测S集中样本点对应的响应值。

6） 采用学习函数H筛选最优样本点ui，具体准则见1.3.1 节，同时根据空间缩减规则选择样本点集D。

7） 进行序贯抽样迭代停止准则判断。若满足所有功能函数的迭代停止条件，则输出代理模型，否则，将筛选的新样本点ui加入训练集，转到步骤3）。若某个功能函数对应的代理模型已达到精度要求，则无需继续对该代理模型进行最优样本点筛选，同时从样本集S中剔除样本集D中的样本点，从而更新S集，转到步骤5），循环直到所有动态代理模型均满足序贯抽样迭代停止准则。

2 基于兴趣子域的可靠性优化方法

2.1 基于SORA 法的概率约束可行性检查方法

SORA 法作为最理想的可靠性优化解耦方法，在工程实际中得到了广泛应用。不过该方法的提出者Du 和Chen[6]指出，该方法值得改进之处在可靠性评估寻找逆最大可能失效点的过程中，并非所有的概率约束都起作用，因此研究识别概率约束可行性检查方法，以此避免不必要的可靠性分析过程可以极大地减少计算量，提高计算效率。如图5所示，图中：G1，G2，G3为3 个功能函数，x1，x2为2 个设计变量，由于最优点H2的位置远离第2 个概率约束的极限状态边界，在对该最优点进行可靠性分析时，第2 个概率约束是不起作用的，所以无需对第2 个概率约束进行可靠性分析。

图5 概率约束条件偏移示意图Fig. 5 Schematic diagram of probability constraint deviation

本节采用ASM 法对SORA 法中的可靠性进行评估，在确定性优化计算得出最优解样本点之后，以最优解样本点为中心进行蒙特卡罗抽样，生成若干个随机样本点，然后分别计算各个概率约束条件的失效概率以及整体失效概率。若某个概率约束条件的失效概率高于设计值，则采用ASM 法计算逆最大可能失效点，从而计算偏移向量，在下一循环的确定性优化中按该偏移向量移动该概率约束条件；若某个概率约束条件的失效概率低于设计值，则直接设置偏移向量为0 的向量，即无需移动该概率约束条件。

本节提出的概率约束可行性检查方法操作简单、易于实现，能极大地减少不必要的可靠性评估过程，提高了计算效率。

2.2 基于兴趣子域的可靠性优化流程

MIGA 算法具备遗传算法的优点，即适用于全局的、离散的、非光滑的优化问题，同时，还克服了传统遗传算法易早熟、局部优化能力弱且计算耗时的缺点。本节将具有全局寻优性能的MIGA算法和基于兴趣子域的Kriging 模型嵌入SORA解耦方法中，通过SORA 法，将可靠性分析解耦成了确定性优化和可靠性分析的优化过程。本节提出的基于兴趣子域的可靠性优化设计流程如图6所示，归纳如下：

图6 基于兴趣子域的可靠性优化流程Fig. 6 Reliability optimization process based on interest subdomain

1） 在确定性空间中，采用自适应空间减缩及序贯抽样策略构造基于兴趣子域的Kriging 模型，具体过程详见1.3.2 节。

2） 设置所有概率约束条件的偏移向量初值为0 的向量。

3） 在确定性优化部分，根据偏移向量V调整约束条件位置，采用MIGA 算法作为优化算法，通过每次调用Kriging 模型计算各个功能函数的响应值，并判断各约束响应是否满足约束条件，然后输出符合约束条件的最优样本点ui。

4） 基于最优样本点ui，根据2.1 节中的概率约束可行性检查方法，筛选出需要进行可靠性分析的概率约束条件，然后，采用ASM 法调用Kriging模型寻找逆最大可能失效点，计算偏移向量，以作为下一次循环过程确定性优化中对应约束条件的偏移向量。

5） 根据2.1 节中计算的整体失效概率Pf判断是否停止迭代，若得到满足约束条件要求的最优可行解，可靠性优化结束；否则，转到步骤3）进行下一次循环。

2.3 数学算例效果校验

本文采用文献[18]中具有非线性概率约束的数学算例，来验证本文所提基于兴趣子域的自适应空间减缩及序贯抽样策略的有效性。该算例的数学模型如下：

式中：ux表示由2 个设计变量x1，x2组成的样本点； Φ为标准正态分布的分布函数；N表示正态分布；u0为初始样本点； μxi为均值。

根据3 个概率约束函数Gi(X)，在平面空间中绘制如图7 所示曲线。

图7 概率约束函数曲线Fig. 7 Curves of probability constraint function

根据1.3.2 节的自适应空间减缩及序贯抽样流程，如图8(a)所示，横、纵坐标分别表示2 个设计变量，采用拉丁超立方抽样方法，在设计空间中生成10 000 个样本点作为样本集S，生成10 个样本点作为训练集T。图中，小圆点代表样本集S，*点代表训练集T。

图8(a)～图8(f)所示为该数学算例的自适应空间缩减过程图。由图可知，随着序贯抽样的进行，代理模型精度不断提高，同时，设计空间不断缩小至兴趣子域。其中，至图8(f)所示在第15 次空间缩减后，共剩余3 265 个样本点，分布区域主要集中在兴趣子域中，训练集38 个样本点，主要分布在约束边界上。

图8 自适应空间缩减过程图Fig. 8 Flowchart of adaptive space reduction

由表1 可知，采用本文提出的基于兴趣子域的可靠性优化方法得到了该数学算例的全局最优解，该最优解与理论解的相对误差为0.066 8%，证明本文方法具有足够的计算精度。且本文方法相比其他方法极大地减少了调用功能函数的次数，其中和文献[18] 中最优的HCC 方法相比减少了40.6%的调用次数，证明本文可靠性优化方法高效。

表1 数学算例可靠度计算结果Table 1 Reliability calculation results of the math example

3 船舶舱段的可靠性优化

3.1 船舶舱段模型介绍

本文将采用基于兴趣子域的可靠性优化设计方法，对一艘618TEU 型多用途船的货舱舱段进行可靠性优化。该船主要装载矿石等散装货物和集装箱，其主尺度等参数如表2 所示。

表2 618TEU 型多用途船主尺度及航速Table 2 Main dimensions and speed of the 618TEU multipurpose vessel

本文优化的区域为该船的中间货舱。由于模型结构、载荷和边界条件都是沿中纵剖面对称分布，为节省计算成本，船宽方向仅选取舱段的1/2 作为计算模型。同时，考虑到减少边界条件的影响，船长方向选取1/2+1+1/2 舱段作为计算模型。舱段模型的垂向范围则选取为船体型深。该舱段有限元模型如图9 所示。

图9 舱段有限元模型Fig. 9 Finite element model of the cabin

舱段模型的坐标规定为：船长方向为X轴，正方向由船艉指向船艏；船宽方向为Y轴，正方向由右舷指向左舷；型深方向为Z轴，正方向由基线指向甲板。

该舱段首、尾端面载荷对称分布，可在两端面建立多点约束（multi-point constraints，MPC），将MPC 与两端面上节点耦合，从而施加约束和载荷。舱段模型边界约束如表3 所示，有限元模型上边界约束情况如图10 所示。表3 中：Ux表示约束X方向的平移；Uy表示约束Y方向的平移；Uz表示约束Z方向的平移；URy表示约束绕Y轴的转动；URz表示约束绕Z轴的转动。

图10 舱段模型约束情况Fig. 10 Constraints of the cabin model

表3 边界条件施加表Table 3 Introduction of boundary conditions

货舱舱段的载荷布置情况按中国船级社的《散货船结构强度直接计算分析指南》[19]确定。有限元模型中施加的载荷主要分为3 个部分：端面弯矩、外部水压力和舱室内货物压力，如图11 所示。

图11 舱段模型的载荷分布情况Fig. 11 Load distribution of the cabin model

3.2 船舶舱段可靠性优化数学模型

图12 舱段横剖面图及设计变量分布情况[13]Fig. 12 Cross section of the cabin and distribution of design variables[13]

次确定性优化迭代后得出的最优设计方案所对应的板厚值；而标准差 σti则设为制造加工的许用偏差，本文取为对应板材厚度均值的2%，即ti～N(μti,(0.02μti)2)。

各设计变量的编号以及对应的名称和初始值如表4 所示。

表4 设计变量参数表Table 4 Parameter list of design variables

式中：Gj(ti)为 各个概率约束条件的功能函数；gj(ti)为按照规范分类的板和梁的应力响应。gj(ti)对应的计算区域及许用应力要求如表5 所示。通过对表5 中6 个结构分类区域分别建立分组，输出每个分组中单元绝对值最大的von Mises应力值作为表5 中6 个变量的应力响应值。

表5 g i(x)对应的计算区域及许用应力要求Table 5 Analysis domain and allowable stress for gi(x)

3.3 船舶舱段可靠性优化计算结果

本文通过Isight 优化平台，集成Matlab，Patran和Nastran 软件，采用2.2 节所述的基于兴趣子域的可靠性优化流程，对此舱段结构进行可靠性优化设计。即采用Isight 自带的MIGA 优化算法，采用Matlab 软件实现SORA 算法，然后采用Patran软件进行参数化建模及后处理，并调用Nastran 软件进行模型分析计算。具体参数设置和流程如下：

1） 采用最优拉丁超立方抽样技术生成180 个初始样本点，对6 个概率约束条件和目标函数分别构造初代Kriging 模型，进而采用基于兴趣子域的自适应空间减缩及序贯抽样策略构造7 个Kriging动态代理模型。

2） 采用SORA 算法结合MIGA 算法进行可靠性评估和优化，具体流程详见2.2 节。并将不考虑概率约束的优化结果与可靠性优化结果进行对比，结果如表6 和表7 所示。

表6 可靠性优化设计与无概率约束优化结果对比Table 6 Comparison between reliability-based design optimization and probability-free constrained optimization

由表7 可知，对船舶舱段进行可靠性优化后，可靠度达99.68%，满足概率约束条件；舱段总质量相比初始方案下降了6.32%，优化效果显著；舱段总质量相比确定性优化增加了0.728%，符合基于可靠性的优化设计原理，即考虑到设计变量的不确定性，通过偏移向量将边界条件朝着可靠域偏移，从而牺牲部分经济性来达到可靠度指标。通过对比文献[13]的可靠性优化结果可知，本文与文献[13]相比少调用了94 次有限元计算次数，且舱段总质量和文献[13]的结果相比减少了0.883%，证明了本文所提可靠性优化策略的高效性和适用性。

表7 可靠性优化设计与无概率约束优化结果对比Table 7 Comparison between reliability-based design optimization and probability-free constrained optimization

采用蒙特卡罗法对上述可靠性优化结果进行验证。在设计空间均匀生成106个随机样本点，采用蒙特卡罗仿真（Monte Carlo simulation，MCS）得到可靠度Pr= 99.67%，由此可知，本文的可靠性优化方案满足可靠度要求。为检验本文中基于兴趣子域的自适应空间减缩及序贯抽样策略所构建的动态代理模型的预测精度，将最终方案进行了有限元计算验证，表8 所示为概率约束条件和优化目标对比结果。

舱段质量作为优化目标，可靠性优化设计和有限元分析两种方式的计算结果分别为456.6 t 和452.536 t，相 对 误 差 为0.898 1，因 此 综 合 表8 可知，本文构建的替代舱段模型的Kriging 模型对各个约束功能函数和优化目标的预测相对误差均在1%以内，验证了本文所提可靠性优化方法的精度，以及其在船舶结构领域的适用性。

表8 可靠性优化结果与有限元模型结果对比Table 8 Comparison between result of reliability-based design optimization and FEM

4 结 论

本文针对船舶结构可靠性优化，采用序贯抽样策略构建了高精度、高效率的动态代理模型，用以达到使用尽可能少的样本点构建具有足够拟合精度的代理模型的目的，从而减少调用耗时的有限元模型计算的次数，节约计算成本。得到的主要结论如下：

1） 通过研究本文所采用SORA 算法中ASM法的寻优原理，以及确定性优化的最优解可能存在的区域，确定了基于可靠性优化的兴趣子域范围，并提出了基于兴趣子域的自适应空间减缩序贯抽样策略，通过制定自适应空间减缩规则，在采用主动学习函数进行序贯抽样的同时，逐步对设计空间进行了缩减，从而能采用尽可能少的样本点构造对兴趣子域进行局部高度拟合的动态代理模型。

2） 针对SORA 法在可靠性评估寻找逆最大可能失效点的过程中并不是所有的概率约束都起作用的问题，提出了概率约束可行性检查方法，该方法操作简单、易于实现，能极大地减少不必要的可靠性评估过程，节省了计算成本。

3） 将基于兴趣子域构建的动态Kriging 代理模型与MIGA 算法相结合嵌入SORA 算法中，构建了基于兴趣子域的可靠性优化方法，通过采用数学算例，验证了该方法的计算效率和精度。

4） 将基于兴趣子域的可靠性优化方法应用于船舶舱段结构可靠性优化设计中，结果显示舱段总质量相比初始方案下降了6.32%，相比确定性优化增加了0.883%，可靠度达到了99.68%。该方案在保证结构安全性的同时提升了经济性，证明了本文所提可靠性优化策略的高效性和适用性。