CN103336869A - 一种基于高斯过程联立mimo模型的多目标优化方法 - Google Patents

Abstract

一种基于高斯过程联立MIMO模型的多目标优化方法，该方法通过实验设计方法得到样本，利用高斯过程联立MIMO模型近似建立设计变量与待考察响应之间的映射关系，采用高斯变异混合遗传算法、蚁群算法等对高斯过程联立MIMO模型进行多目标优化，得到关于设计变量组合的Pareto前沿，进而判断Pareto前沿中解集的分布质量，最后根据具体要求选择Pareto前沿中某一设计变量组合进行高精度分析求解，当得到的结果满意，则进行物理实验。本发明将实验设计、高精度分析求解、代理模型技术以及优化算法集成应用于优化设计中，能够大大缩减优化设计的耗时和计算成本，显著提高工作效率。

Description

一种基于高斯过程联立MIMO模型的多目标优化方法

技术领域

本发明涉及一种代理模型技术结合优化算法的多目标优化方法，尤其涉及一种高斯过程联立建模结合优化算法的多目标优化设计方法。

背景技术

代理模型技术的发展为优化技术的引入提供了有效的桥梁，通过先进实验设计方法得到样本数据，利用代理模型技术近似建立设计变量与待考察响应之间的映射关系，最后运用优化算法求解该关系模型得到最优设计变量组合。这种将实验设计、高精度分析求解、代理模型技术以及优化算法集成应用于优化设计中的方法，不仅能够极大减少变量修改中的人为因素，而且能够在复杂的设计方案中得到尽可能完善的设计变量组合。

现有基于代理模型技术的多目标优化方法中，有以下几种：

（1）基于Kriging模型的多目标优化方法，如Kaymaz I在：“Application of kriging methodto structural reliability problems（Kriging模型在结构可靠性问题中的应用）”（Structural Safety,2005,27（2）:133-151）（结构安全性，2005年27期第2卷，133-151页）进行优化时，核函数参数的寻优过程会造成大量的计算成本，这在高维的各向异性问题中尤为明显，这一缺陷导致建立Kriging模型所耗费的时间相对其他建模方法而言要多很多。

（2）基于传统独立高斯过程回归模型的多目标优化方法，如Jin Yuan在：“Reliablemulti-objective optimization of high-speed WEDM process based on Gaussian process regression（基于高斯过程回归的高速线切割加工可靠性多目标优化）”（International Journal of MachineTools&Manufacture,2008（48）:47–60）（机床与制造国际期刊，2008年48期，47-60页）进行多目标优化时，实际上是采用传统高斯过程回归模型分别独立地对每一个输出端建立高斯过程模型，在对同一组未知输入向量的各个响应进行拟合或预测时，不同输出端对应的各模型所利用的估计参数也是不同的，因此，传统高斯过程模型本质上可以看作一组相对独立的多输入—单输出（multiple inputs and single output,MISO）模型，无法充分利用数据信息。

（3）基于限定性相关矩阵的代理模型的多目标优化方法，如McMillian N J在：“Analysisof Protein Activity Data by Gaussian Stochastic Process Models（基于高斯随机过程模型的蛋白质活动分析）”（Journal of Biopharmaceutical Statistics,1999（9）:145-160）（生物制药统计学期刊，1999年9期，145-160页）预测结果中，限定性相关矩阵赋予类别输入变量某些特定的假设参数化关系，虽然能够大大简化相关计算的复杂程度，但是限定性相关矩阵通常只能用于特定的过程或者系统，因此缺乏捕捉各输出端之间不同类型相关性的柔性，如乘积相关矩阵只能捕捉类别输入变量水平间的正相关性，而缺乏捕捉其负相关性的能力。

综上所述，现有代理模型多目标优化方法中，不足主要体现在以下四方面：1、存在计算耗时长的问题；2、在假设所有的输入信息仅有量化了的设计变量（数值输入变量）的前提下，通过计算未知点与已知点的相对位置关系，确定未知点在数值输入变量空间中的具体位置（数值输入变量的空间相关性），继而进行模型的拟合或预测，然而事实上，需要考察的输出端个数也是已知的，且每个输出端在性质上是分类的（类别输入变量）并能决定代理模型的一些重要性质，如果把系统中的已知条件都看作建模时的输入信息，那么传统多目标建模方法没有综合考虑这两种类型的输入变量，因此常常未能充分利用数据信息；3、许多工程问题的多目标模型中的输出端通常来说都承载了彼此之间的一些相关性，传统多目标建模方法忽略了在同一组数值输入向量下，不同输出端的响应之间存在的相关性（类别输入变量的交叉相关性），由于没有对类别输入变量的交叉相关性进行建模，因此无法利用该相关性来提高模型对未知点的拟合或预测精度；4、缺乏捕捉各输出端之间不同类型相关性的柔性。

发明内容

本发明的目的在于提出一种基于高斯过程联立MIMO（multiple inputs and multipleoutputs）模型的多目标优化方法，凭借高斯过程联立MIMO模型的模型相关参数较少、易实现、样本规模较小的优势，解决优化计算耗时长、系统信息利用不充分、多目标代理模型无法利用类别输入变量的交叉相关性来提高模型对未知点拟合或预测的精度以及无法捕捉输出端之间不同类型相关性的问题。

本发明实现上述目的所采取的技术方案是：一种基于高斯过程联立MIMO模型的多目标优化方法，包括实验设计方法抽样，标准高斯过程建模及多目标优化算法寻优与验证，具体按照如下步骤进行：

步骤1，根据用户给定的初始条件，通过增强型平移传播拉丁超立方ETPLHD抽样方法进行实验设计，得到设计变量组合X=[x₁,x₂,...,x_N]^T，并对设计变量组合分别进行高精度分析求解，得到相应的响应值y=[y₁,y₂…,y_N]^T，与设计变量组合一并组成模型的训练数据集（X,y）；

步骤2，以步骤1得到的训练数据集，建立基于非限定性相关矩阵的高斯过程联立MIMO模型，确定基于非限定性相关矩阵的高斯过程联立MIMO模型中的协方差函数，该协方差函数的计算公式如下公式（1）：

$\mathrm{Cov}(y\left(h_r\right),y\left(h_s\right))={\mathrm{\τ}}_{z_r,z_s}\left\{v_1\exp \right[-\frac{1}{2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d{(x_p^d-x_q^d)}^2]+{\mathrm{\δ}}_{\mathrm{pq}}{v}_0\}---\left(1\right)$

式中h=（x^t,z_c）^t，其中x^t是数值输入变量，z为类别输入变量，并以c=1,...,m作为标记，c表示模型输出端的个数，

是输出端z_r和z_s之间的交叉关联度大小，定义非限定性相关矩阵T={τ_r,s}，设θ=[w₁,...,w_D,v₀,v₁]^T是指数平方平稳协方差函数的超参数，v₀是白噪音估计方差，v₁是纵坐标上的估计方差，D是设计变量的维数，δ_pq是Kronecker算子；

并对所述的协方差函数中的模型参数

θ、v₀和v₁进行估计和优化，估计和优化的方法是通过最大化边界似然函数法，转化为求解最小化边界log-似然函数；

步骤3，采用高斯变异混合遗传算法GGA-h、多目标粒子群优化算法和蚁群算法对步骤2所建立的基于非限定性相关矩阵的高斯过程联立MIMO模型进行多目标优化；

步骤4，判断步骤3中多目标优化得到的Pareto前沿中解的分布是否均匀，解的个数是否足够以及边界处的解是否接近单目标极值而具有良好拓展性，当前沿质量达标时，保存该Pareto前沿，否则返回重新设定多目标优化的算法参数；

步骤5，在步骤4得到符合要求的Pareto前沿的前提下，根据用户具体要求选择Pareto前沿中某一设计变量组合进行高精度分析求解，当得到的结果满意时，则为优化结果，否则返回步骤1中重新进行ETPLHD抽样实验设计。

步骤2所述的确定基于非限定性相关矩阵的高斯过程联立MIMO模型中的协方差函数，包括有效定义的数值输入变量空间相关性和类别输入变量交叉相关性，首先根据用户要求的具体设计变量维数D赋予相应超参数的初值个数为D+2，即增加v₀和v₁，然后根据具体输出端个数c确定非限定性相关矩阵T的规模，即c×c，并赋予初值，具体方法为该矩阵主对角线元素的自相关系数为1，其余矩阵元素为0.01，高斯过程联立MIMO模型相关参数个数为（D+2）+c（c-1）/2。

所述的非限定性相关矩阵T，其超球面分解参数化步骤如下：

步骤A，对一个非限定性相关矩阵T进行乔里斯基分解（Cholesky-type decomposition），得到公式（2）：

T=LL^t                    （2）

式中L={l_r,s}是一个主对角线元素严格正定的下三角矩阵；

步骤B，每个行向量（l_r,1,...,l_r,r）根据一个r-维单位超球面体的曲面坐标点来定义，根据球坐标系，当r=1时，l_1,1=1，r=2,...,m，得到公式（3）：

$\left\{\begin{array}{c}{l}_{r,1}=\cos \left({\mathrm{\φ}}_{r,1}\right),\\ l_{r,s}=\sin \left({\mathrm{\φ}}_{r,1}\right)\mathrm{\Λ}\sin \left({\mathrm{\φ}}_{r,s-1}\right)\cos \left({\mathrm{\φ}}_{r,s}\right),s=2,K,r-1,\\ l_{r,r}=\sin \left({\mathrm{\φ}}_{r,1}\right)\mathrm{\Λ}\sin \left({\mathrm{\φ}}_{r,r-2}\right)\sin \left({\mathrm{\φ}}_{r,r-1}\right),\end{array}\right.---\left(3\right)$

式中φ_r,s∈(0,π)，所有的φ_r,s记作Φ，且φ_r,s的取值约束在（0,π）之间，即将优化问题中的矩阵正定性不等式约束参数化转变到易于计算的边界约束。

步骤2所述的将模型参数

θ、v₀和v₁进行估计和优化的具体步骤如下：

步骤a，将高斯过程联立MIMO模型中的协方差矩阵K_m表示成Kronecker形式

为标准高斯过程协方差矩阵，根据克罗内克积的特性，K_m转变为一个Nm×Nm的分块矩阵；

步骤b，按照公式（4）、（5）和（6）对协方差矩阵K_m进行一系列矩阵计算和变换；

${K_m}^{-1}=T^{-1}\⊗K^{-1}---\left(4\right)$

${{\mathrm{EK}}_m}^{-1}=E(T^{-1}\⊗K^{-1})=(E\⊗1)(T^{-1}\⊗K^{-1})=\left({\mathrm{EK}}^{-1}\right)\⊗T^{-1}---\left(5\right)$

$\mathrm{tr}\left({{\mathrm{EK}}_m}^{-1}\right)=\mathrm{tr}({\mathrm{ET}}^{-1}\⊗K^{-1})=\mathrm{tr}\left({\mathrm{EK}}^{-1}\right)\mathrm{tr}\left(T^{-1}\right)---\left(6\right)$

从公式（6）可知，非限定性相关矩阵T和标准高斯过程协方差矩阵K是相互独立的，而T经过超球面分解与Φ一一对应，通过公式（7）和（8）可分别得到参数θ和Φ。

$\frac{\∂L}{\∂\mathrm{\θ}}=-\frac{1}{2}\mathrm{tr}\left({k_m}^{-1}\frac{{\∂K}_m}{\∂\mathrm{\θ}}\right)+\frac{1}{2}{y}^T{K_m}^{-1}\frac{\∂K_m}{\∂\mathrm{\θ}}{K_m}^{-1}y---\left(7\right)$

$\frac{\∂L}{\∂\mathrm{\Φ}}=-\frac{1}{2}\mathrm{tr}\left({k_m}^{-1}\frac{{\∂K}_m}{\∂\mathrm{\Φ}}\right)+\frac{1}{2}{y}^T{K_m}^{-1}\frac{\∂K_m}{\∂\mathrm{\Φ}}{K_m}^{-1}y---\left(8\right)$

与现有技术相比，本发明具备的有益效果是：提出一种简单灵活的基于高斯过程联立MIMO模型的多目标优化方法，以标准高斯过程回归模型的重要概念和模型的公式实现为基础，分别以距离基函数定义数值输入变量的空间相关性，以相关矩阵定义类别输入变量的交叉相关性，充分利用各类过程或系统信息，构建基于非限定性相关矩阵的高斯过程联立MIMO模型中协方差函数的一般方法；进而采用超球面分解参数化了非限定性相关矩阵，保持应用柔性的同时，将优化问题中复杂的正定性约束转化为较为简单的边界约束，大大降低计算复杂度，显著减少计算耗时。

附图说明

图1为本发明具体实施例1中IG建模方法得到的Pareto前沿。

图2为本发明具体实施例1中SG建模方法得到的Pareto前沿。

图3为本发明具体实施例1中解析式传回法得到的Pareto前沿。

图4为本发明具体实施例中二维超球面分解示意图。

图5为本发明具体实施例中三维超球面分解示意图。

图6为本发明具体实施例二中汽车桥壳后盖剖视图。

图7为本发明具体实施例二中汽车桥壳后盖正二测视图。

图8为本发明具体实施例二中优化流程图。

图9为本发明具体实施例二中本基于高斯过程联立MIMO模型的多目标优化方法得到的Pareto前沿。

图10为本发明具体实施例二中优化后的成形极限图。

图11为本发明具体实施例二中材料流入情况图。

图12为本发明具体实施例二中料厚的样件值与模拟值。

具体实施方式

下面结合实施例及其附图详细说明基于高斯过程联立MIMO模型的多目标优化方法的技术方案，但本发明的具体实施方式不局限于下述的实施例。

在实施例1中将通过一个数值算例来证明本发明所提出的基于高斯过程联立MIMO模型的多目标优化方法的有效性和可行性。为了达到比较的目的，将考虑下述两种建模方法：

A）独立高斯过程模型（individual Gaussian process），记作IG。该方法实际上就是通过标准高斯过程模型分别对每一个需要考察的目标建模。

B）联立高斯过程模型（simultaneous Gaussian process），记作SG。即通过本发明所述的高斯过程联立MIMO模型对目标进行多目标联立建模。

实施例1

解决多目标问题的数值算例

一种基于高斯过程联立MIMO模型的多目标优化方法，包括实验设计方法抽样，标准高斯过程联立MIMO模型及多目标优化算法寻优，包括如下步骤：

步骤1，根据用户给定的数值多目标问题，其数学描述如公式（9），通过拉丁超立方实验设计方法进行实验设计，在变量x∈[0,1]中抽样得到20个样本点，并代入公式（9）的解析式f₁，g（x₁,x₂）中，得到相应的响应值，与变量x的样本集一并组成模型的训练数据集；

f₁=x₁

f₂=g(x₁,x₂)h(x₁,x₂)

g(x₁,x₂)=11+x₂ ²-10cos(2px₂)                    （9）

$h(x_1,x_2)=\left\{\begin{array}{c}1-\sqrt{f_1/g,}{f}_1<g\\ 0,\mathrm{else}\end{array}\right.$

(x∈(0,1))

步骤2，以步骤1得到的训练数据集，分别通过IG和SG建模方法对本实施例中的解析式f₁，g（x₁,x₂）构建代理模型，其中包括分别确定模型中的协方差函数，IG和SG建模方法的协方差函数分别如公式（10）和公式（11）所示，并对两个协方差函数中的模型相关参数（

θ,v₀,v₁）估计和优化，估计和优化的方法是通过最大化边界似然函数法，转化为求解最小化边界log-似然函数；

$\mathrm{Cov}(y\left(x_p\right),y\left(x_q\right))=\left\{v_1\exp \right[-\frac{1}{2}\underset{d=1}{\overset{1}{\mathrm{\&Sigma;}}}{w}_d{(x_p^d-x_q^d)}^2]+{\mathrm{\&delta;}}_{\mathrm{pq}}{v}_0\}---\left(10\right)$

$\mathrm{Cov}(y\left(h_r\right),y\left(h_s\right))={\mathrm{\&tau;}}_{z_1,z_2}\left\{v_1\exp \right[-\frac{1}{2}\underset{d=1}{\overset{1}{\mathrm{\&Sigma;}}}{w}_d{(x_p^d-x_q^d)}^2]+{\mathrm{\&delta;}}_{\mathrm{pq}}{v}_0\}---\left(11\right)$

式中h=（x^t,z₂）^t，其中x^t是数值输入变量，z为类别输入变量，并以c=2作为标记，c表示模型输出端的个数，

是输出端z₁和z₂之间的交叉关联度大小，定义T={τ_1,2}，设θ=[w₁,v₀,v₁]^T是指数平方平稳协方差函数的超参数，v₀是白噪音估计方差，v₁是纵坐标上的估计方差，变量的维数为1，δ_pq是Kronecker算子；

步骤3，采用第二代遗传算法NSGA II分别对步骤2中IG和SG建模方法所建立的代理模型进行多目标优化，算法参数设定如下：种群大小为50，交叉概率为0.8，变异概率为0.3，迭代次数为200，如图1为IG建模方法得到的Pareto前沿，图2为SG建模方法得到的Pareto前沿；

步骤4，在使用第二代遗传算法解决本实施例中的多目标问题时，算法寻优过程中，是将每一次迭代产生的新种群传回到公式（9）的解析式f₁，g（x₁,x₂）中，根据响应值进行相应的适应度初值计算并进行非劣排序（本算例中称为解析式传回法），使用解析式传回法得到的标准Pareto前沿，如图3。对比图1、图2、图3，可知步骤3中应用IG建模方法时得到的图1中Pareto前沿在边界处出现了很大的误差，而且解的分布性也不均匀；而图2中应用SG建模方法时得到的Pareto前沿边界附近的解比较准确，虽然边界附近的分布性稍稍差于图3，但是前沿整体的吻合程度要比IG高得多，且解的分布更加均匀，在Pareto前沿中最常使用的前沿中段部分，其解的分布性甚至优于图3；

步骤5，由步骤4得到的Pareto前沿结果分析中，可以证明本发明提出的基于高斯过程联立MIMO模型的多目标优化方法比基于标准高斯过程建模方法具有更高的预测精度；且能够充分利用数据中的信息；并在解决多目标问题中有着更加良好的表现。

在步骤2中，SG建模方法的高斯过程联立MIMO模型中的协方差函数，此时的具体变量维数为1赋予相应超参数的初值个数为1+2，即增加v₀和v₁，且具体输出端个数2，可确定非限定性相关矩阵T₂的规模，即2×2，并赋予初值，如公式（12），具体方法为该矩阵主对角线元素的自相关系数为1，其余矩阵元素为0.01，高斯过程联立MIMO模型相关参数个数为4。

$T_2=\left[\begin{array}{cc}1& 0.01\\ 0.01& 1\end{array}\right]---\left(12\right)$

非限定性相关矩阵T₂，其超球面分解参数化步骤如下：

步骤A，对一个非限定性相关矩阵T进行乔里斯基分解（Cholesky-type decomposition），得到公式（13）

式中L={l_r,s}是一个主对角线元素严格正定的下三角矩阵；

步骤B，每个行向量（l_2,1,l_2,2）根据一个2-维单位超球面体的曲面坐标点来定义，此时输出端为2维，则其分解结果是一个半球面，如图4所示；若输出端为3维，则其分解结果是一个半球体，如图5所示；若输出端超过3维，则其分解结果为超球面，结果无法以图形形式表示。本实施例中根据如下球坐标系，当r=1，l_1,1=1，r=2，得到公式（14）

$\left\{\begin{array}{c}{l}_{21}=\cos \left({\mathrm{\&phi;}}_{21}\right)\\ l_{22}=\sin \left({\mathrm{\&phi;}}_{21}\right)\end{array}\right.---\left(14\right)$

式中φ_r,s∈(0,π)，所有的φ_r,s记作Φ，且φ_r,s的取值约束在（0,π）之间，即将优化问题中的矩阵正定性不等式约束参数化转变到易于计算的边界约束。

在步骤2所述的模型相关参数（

θ,v₀,v₁）估计和优化步骤如下：

步骤a，将高斯过程联立MIMO模型中的协方差矩阵K₂表示成Kronecker形式

K为标准高斯过程协方差矩阵，根据克罗内克积的特性，K₂转变为一个40×40的分块矩阵；

步骤b，对K₂进行一系列矩阵计算和变换，如公式（15）、（16）和（17）所示：

${K_2}^{-1}=T^{-1}\&CircleTimes;K^{-1}---\left(15\right)$

${{\mathrm{EK}}_2}^{-1}=E(T^{-1}\&CircleTimes;K^{-1})=(E\&CircleTimes;1)(T^{-1}\&CircleTimes;K^{-1})=\left({\mathrm{EK}}^{-1}\right)\&CircleTimes;T^{-1}---\left(16\right)$

$\mathrm{tr}\left({{\mathrm{EK}}_2}^{-1}\right)=\mathrm{tr}({\mathrm{ET}}^{-1}\&CircleTimes;K^{-1})=\mathrm{tr}\left({\mathrm{EK}}^{-1}\right)\mathrm{tr}\left(T^{-1}\right)---\left(17\right)$

从公式（17）可知，非限定性相关矩阵T和标准高斯过程协方差矩阵K是相互独立的，而T经过超球面分解与Φ一一对应，则有公式（18）和（19）：

$\frac{\&PartialD;L}{\&PartialD;\mathrm{\&theta;}}=-\frac{1}{2}\mathrm{tr}\left({K_2}^{-1}\frac{{\&PartialD;K}_2}{\&PartialD;\mathrm{\&theta;}}\right)+\frac{1}{2}{y}^T{K_2}^{-1}\frac{\&PartialD;K_2}{\&PartialD;\mathrm{\&theta;}}{K_2}^{-1}y---\left(18\right)$

$\frac{\&PartialD;L}{\&PartialD;\mathrm{\&Phi;}}=-\frac{1}{2}\mathrm{tr}\left({K_2}^{-1}\frac{{\&PartialD;K}_2}{\&PartialD;\mathrm{\&Phi;}}\right)+\frac{1}{2}{y}^T{K_2}^{-1}\frac{\&PartialD;K_2}{\&PartialD;\mathrm{\&Phi;}}{K_2}^{-1}y---\left(19\right)$

因此模型相关参数参数θ和Φ可以通过（18）和（19）两个公式在估计和优化过程中分别得到。

实施例2

基于高斯过程联立MIMO模型的多目标优化方法在板料拉深中的应用

采用基于高斯过程联立MIMO模型的多目标优化方法优化某汽车桥壳后盖，如图6和图7所示，拉深工艺参数多目标优化，其优化过程流程图如图8所示，包括如下步骤：

步骤1，根据客户制订的工艺文件，选定待优化的成形工艺参数，并确定其取值范围如下：80≤压边力BHF（kN）≤500，275≤坯料直径D（mm）≤325，0.10≤板料与凹模的静摩擦系数μ₁≤0.15，0.10≤板料与压边圈的静摩擦系数μ₂≤0.15，通过增强型平移传播拉丁超立方（ETPLHD）抽样方法进行实验设计，得到30组工艺参数组合，并利用板料成形CAE软件根据每一组样本所对应的工艺参数组合分别建立有限元模型进行数值模拟，根据成形极限图FLD，定义得到拉裂与起皱目标函数，如公式（20）（21）所示，再根据每一组工艺参数组合下的数值模拟结果数据分别计算目标函数值D_C和D_W，各组响应值与工艺参数组合一并组成高斯过程联立MIMO模型的训练数据集，见表1；

$D_C=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left\{\begin{array}{cc}{({\mathrm{\&epsiv;}}_1^i-\mathrm{\&phi;}\left({\mathrm{\&epsiv;}}_2^i\right))}^2& ,{\mathrm{\&epsiv;}}_1^i>\mathrm{\&phi;}\left({\mathrm{\&epsiv;}}_2^i\right)\\ 0& ,{\mathrm{\&epsiv;}}_1^i<\mathrm{\&phi;}\left({\mathrm{\&epsiv;}}_2^i\right)\end{array}\right.---\left(20\right)$

$D_W=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left\{\begin{array}{cc}{(\mathrm{\&eta;}\left({\mathrm{\&epsiv;}}_2^i\right)-{\mathrm{\&epsiv;}}_1^i)}^2& ,{\mathrm{\&epsiv;}}_1^i<\mathrm{\&eta;}\left({\mathrm{\&epsiv;}}_2^i\right)\\ 0& ,{\mathrm{\&epsiv;}}_1^i>\mathrm{\&eta;}\left({\mathrm{\&epsiv;}}_2^i\right)\end{array}\right.---\left(21\right)$

表1ETPLHD抽样样本及其响应

序号	BHF（kN）	D（mm）	μ1	μ2	DC	DW
							1	110.6	314	0.115	0.129	0.0493	6.1528
2	431.3	317	0.109	0.143	1.2988	3.0485
							3	176.2	278	0.105	0.112	0.0193	12.7052
4	81.2	308	0.102	0.122	0.0359	7.6124
							5	463.1	300	0.125	0.108	0.687	3.5504
6	302.6	290	0.14	0.101	0.3651	5.8154
							7	356.1	301	0.11	0.14	0.3099	4.2335
8	282.3	318	0.136	0.11	0.9198	3.9821
							9	248.6	309	0.123	0.114	0.2748	4.879
10	472.5	297	0.124	0.106	0.7327	3.5862
							11	455.8	286	0.12	0.131	0.7507	4.1728
12	486.4	311	0.127	0.135	1.5908	2.54
							13	192.6	289	0.128	0.144	0.0983	7.8043
14	269.2	294	0.139	0.102	0.1724	5.9838
							15	390.8	298	0.1	0.118	0.1743	4.5515
16	190	276	0.146	0.145	0.0706	10.3085
							17	424.1	324	0.149	0.132	2.7094	2.4166
18	221.2	284	0.147	0.147	0.2274	7.0869
							19	404.9	282	0.112	0.12	0.1629	5.684
20	332.6	306	0.134	0.115	0.7532	3.9816
							21	322.8	323	0.107	0.127	0.5417	3.518
22	387.5	319	0.103	0.104	0.2413	3.7359
							23	124.1	321	0.119	0.149	0.1258	5.1912
24	206.9	293	0.114	0.133	0.0885	7.1054

25	140.5	313	0.131	0.138	0.1432	5.7626
							26	363.4	288	0.137	0.139	0.3734	4.604
27	153.6	279	0.13	0.107	0.0327	12.3305
							28	312.5	303	0.143	0.123	0.8352	4.3467
29	238.8	281	0.117	0.119	0.0568	7.2013
							30	105.6	305	0.142	0.125	0.1867	7.4352

步骤2，以步骤1得到的训练数据集，建立基于非限定性相关矩阵的高斯过程联立MIMO模型，其中包括确定基于非限定性相关矩阵的高斯过程联立MIMO模型中的协方差函数，如公式（22）所示，并对该协方差函数中的模型相关参数（

θ,v₀,v₁）估计和优化，估计和优化的方法是通过最大化边界似然函数法，转化为求解最小化边界log-似然函数；

$\mathrm{Cov}(y\left(h_1\right),y\left(h_2\right))={\mathrm{\&tau;}}_{z_1,z_2}\left\{v_1\exp \right[-\frac{1}{2}\underset{d=1}{\overset{4}{\mathrm{\&Sigma;}}}{w}_d{(x_p^d-x_q^d)}^2]+{\mathrm{\&delta;}}_{\mathrm{pq}}{v}_0\}---\left(22\right)$

式中h=（x^t,z₂）^t，其中x^t是数值输入变量，z为类别输入变量，并以c=2作为标记，c表示模型输出端的个数，

是输出端z₁和z₂之间的交叉关联度大小，定义T={τ_1,2}，设θ=[w₁,...,w₄,v₀,v₁]^T是指数平方平稳协方差函数的超参数，v₀是白噪音估计方差，v₁是纵坐标上的估计方差，设计变量的维数为4，δ_pq是Kronecker算子；

步骤3，采用高斯变异混合遗传算法（GGA-h）算法参数设定如下：种群大小为50，交叉概率为0.9，变异概率为0.1，迭代次数为300，如图9所示，为基于高斯过程联立MIMO模型的多目标优化方法得到的Pareto前沿；

步骤4，由步骤3中的Pareto前沿，解的分布均匀情况良好，边界处的解比较接近单目标极值，说明解集有良好的拓展性。与此同时，能够知道本例中拉裂和起皱存在非线性关系，而交叉关联度

则说明两者存在负相关，但不是完全负相关，即并非拉裂目标函数值减小越快，起皱目标函数值就增大越快，当受压达到一定程度时，板料成形失稳，由面内应变转变为弯曲变形时，材料流动受阻，也会导致拉裂缺陷增加，这也与实际情况相符。为了验证本发明中的优化方法得到的结果的准确性，表2为图9所示Pareto前沿中的某3组工艺参数组合的模型算法预测值和计算机辅助工程（CAE）验证值，由对比数据可知模型算法预测值与CAE验证值十分接近，虽有少许偏差，但其变化趋势一致，因此优化结果可以指导该桥壳后盖的生产，保存该Pareto前沿；

表2对比数据验证

步骤5，在步骤4得到了符合要求的Pareto前沿的前提下，表3是从该Pareto前沿中选择的与初始样本点中选出的拉裂目标函数值相近的3组数据作对比，Pareto前沿中非劣解的起皱目标函数值均比初始样本中的该值有不同程度的降低。当从该Pareto前沿和初始样本中选择起皱目标函数值相近的数据时，也可得到相同的结论。这表明本发明的优化方法能够同时优化拉裂和起皱两个质量指标。根据用户制订的工艺文件，选取如图9所示的P点作为满意解，其D_C和D_W值分别为0.01307和4.518，设计变量BHF、D、μ₁和μ₂所对应的值分别为114.7kN，286mm，0.127和0.119，利用该工艺参数组合进行数值模拟，根据得到的成形极限图图10所示，能够看出零件的拉裂和起皱趋势都很小。材料的流入情况如图11所示，可以看出拉深成形过程中，料流状况良好，材料流入均匀，最大最小流入量差值为3.22mm。根据优化的如图9所示的P点工艺参数组合，做相应的微调，使其满足实际生产设置要求，样件成形质量良好。对样件在对称轴处进行线切割，量取样件对称轴上圆角、直边等特征起始位置以及中段均分点处24个点的料厚值，并在数值模拟结果中提取相应位置的料厚值，样件值与模拟值如图12，两者的分布和变化趋势基本是一致的，同时说明了样件的料厚分布情况良好。根据用户制订的检测计划，对试制的5个样件进行尺寸测量，结果见表4。

表3拉裂目标函数值D_C相当时，起皱目标函数D_W优化值和样本值对比数据

表4拉深成形样件检验报告

在步骤2中，本发明的优化方法中的高斯过程联立MIMO模型中的协方差函数，此时的具体设计变量维数为4，赋予相应超参数的初值个数为4+2，即增加v₀和v₁，且具体输出端个数2，可确定非限定性相关矩阵T₂的规模，即2×2，并赋予初值，如公式（12），具体方法为该矩阵主对角线元素的自相关系数为1，其余矩阵元素为0.01，高斯过程联立MIMO模型相关参数个数为7。

非限定性相关矩阵T₂，其超球面分解参数化步骤如下：

步骤A，对一个非限定性相关矩阵T进行乔里斯基分解（Cholesky-type decomposition），得到公式（13），式中L={l_r,s}是一个主对角线元素严格正定的下三角矩阵；

步骤B，每个行向量（l_2,1,l_2,2）根据一个2-维单位超球面体的曲面坐标点来定义，此时输出端为2维，则其分解结果是一个半球面，如图4所示；若输出端为3维，则其分解结果是一个半球体，如图5所示；若输出端超过3维，则其分解结果为超球面，结果无法以图形形式表示。本实施例中根据如下球坐标系，当r=1，l_1,1=1，r=2，得到公式（14），式中φ_r,s∈(0,π)，所有的φ_r,s记作Φ，且φ_r,s的取值约束在（0,π）之间，即将优化问题中的矩阵正定性不等式约束参数化转变到易于计算的边界约束。

在步骤2所述的模型相关参数（θ,v₀,v₁）估计和优化步骤如下：

步骤a，将高斯过程联立MIMO模型中的协方差矩阵K₂表示成Kronecker形式为标准高斯过程协方差矩阵，根据克罗内克积的特性，K₂转变为一个40×40的分块矩阵；

步骤b，对K₂进行一系列矩阵计算和变换，如公式（15）（16）（17）所示，从公式（17）可知，非限定性相关矩阵T和标准高斯过程协方差矩阵K是相互独立的，而T经过超球面分解与Φ一一对应，则有公式（18）（19），因此模型相关参数参数θ和Φ可以通过（18）（19）两个公式在估计和优化过程中分别得到。

由此可见，高斯过程模型有严格的统计学理论基础，对小样本、非线性以及高维问题有着良好的适应性，且实现过程相对简单。因此本发明在标准高斯过程回归模型的基础上，提出了一种简单灵活的基于高斯过程联立MIMO模型的多目标优化方法。分别以距离基函数定义了数值输入变量的空间相关性，以相关矩阵定义了类别输入变量的交叉相关性，构建了基于非限定性相关矩阵的高斯过程联立MIMO模型中协方差函数的一般方法；进而采用超球面参数化分解方法构建非限定性相关矩阵，保持了应用柔性的同时，将优化问题中复杂的正定性约束转化为较为简单的边界约束，大大降低计算复杂度。通过两个实施例，证明基于高斯过程联立MIMO模型的多目标优化方法在模型精度、输出端相关性捕捉以及在解决多目标优化问题中具有非常明显的优势。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于高斯过程联立MIMO模型的多目标优化方法，包括实验设计方法抽样，标准高斯过程建模及多目标优化算法寻优与验证，其特征在于，具体按照如下步骤进行： 

步骤1，根据用户给定的初始条件，通过增强型平移传播拉丁超立方ETPLHD抽样方法进行实验设计，得到设计变量组合X=[x₁,x₂,...,x_N]^T，并对设计变量组合分别进行高精度分析求解，得到相应的响应值y=[y₁,y₂...,y_N]^T，与设计变量组合一并组成模型的训练数据集（X,y）； 

步骤2，以步骤1得到的训练数据集，建立基于非限定性相关矩阵的高斯过程联立MIMO模型，确定基于非限定性相关矩阵的高斯过程联立MIMO模型中的协方差函数，该协方差函数的计算公式如下公式（1）： 

式中h=（x^t,z_c）^t，其中x^t是数值输入变量，z为类别输入变量，并以c=1,...,m作为标记，c表示模型输出端的个数，是输出端z_r和z_s之间的交叉关联度大小，定义非限定性相关矩阵T={τ_r,s}，设θ=[w₁,...,w_D,v₀,v₁]^T是指数平方平稳协方差函数的超参数，v₀是白噪音估计方差，v₁是纵坐标上的估计方差，D是设计变量的维数，δ_pq是Kronecker算子； 

并对所述的协方差函数中的模型参数θ、v₀和v₁进行估计和优化，估计和优化的方法是通过最大化边界似然函数法，转化为求解最小化边界log-似然函数； 

步骤3，采用高斯变异混合遗传算法GGA-h、多目标粒子群优化算法和蚁群算法对步骤2所建立的基于非限定性相关矩阵的高斯过程联立MIMO模型进行多目标优化； 

步骤4，判断步骤3中多目标优化得到的Pareto前沿中解的分布是否均匀，解的个数是否足够以及边界处的解是否接近单目标极值而具有良好拓展性，当前沿质量达标时，保存该Pareto前沿，否则返回重新设定多目标优化的算法参数； 

步骤5，在步骤4得到符合要求的Pareto前沿的前提下，根据用户具体要求选择Pareto前沿中某一设计变量组合进行高精度分析求解，当得到的结果满意时，则为优化结果，否则返回步骤1中重新进行ETPLHD抽样实验设计。 

2.根据权利要求1所述的基于高斯过程联立MIMO模型的多目标优化方法，其特征在于，步骤2所述的确定基于非限定性相关矩阵的高斯过程联立MIMO模型中的协方差函数，包括 有效定义的数值输入变量空间相关性和类别输入变量交叉相关性，首先根据用户要求的具体设计变量维数D赋予相应超参数的初值个数为D+2，即增加v₀和v₁，然后根据具体输出端个数c确定非限定性相关矩阵T的规模，即c×c，并赋予初值，具体方法为该矩阵主对角线元素的自相关系数为1，其余矩阵元素为0.01，高斯过程联立MIMO模型相关参数个数为（D+2）+c（c-1）/2。 

3.根据权利要求2所述的基于高斯过程联立MIMO模型的多目标优化方法，其特征在于，所述的非限定性相关矩阵T，其超球面分解参数化步骤如下： 

步骤A，对一个非限定性相关矩阵T进行乔里斯基分解（Cholesky-type decomposition），得到公式（2）： 

T=LL^t                    （2） 

式中L={l_r,s}是一个主对角线元素严格正定的下三角矩阵； 

步骤B，每个行向量（l_r,1,...,l_r,r）根据一个r-维单位超球面体的曲面坐标点来定义，根据球坐标系，当r=1时，l_1,1=1，r=2,...,m，得到公式（3）： 

式中φ_r,s∈(0,π)，所有的φ_r,s记作Φ，且φ_r,s的取值约束在（0,π）之间，即将优化问题中的矩阵正定性不等式约束参数化转变到易于计算的边界约束。 

4.根据权利要求1或2所述的基于高斯过程联立MIMO模型的多目标优化方法，其特征在于，步骤2所述的将模型参数

θ、v₀和v₁进行估计和优化的具体步骤如下： 

步骤a，将高斯过程联立MIMO模型中的协方差矩阵K_m表示成Kronecker形式 

K为标准高斯过程协方差矩阵，根据克罗内克积的特性，K_m转变为一个Nm×Nm的分块矩阵； 

步骤b，按照公式（4）、（5）和（6）对协方差矩阵K_m进行一系列矩阵计算和变换； 

从公式（6）可知，非限定性相关矩阵T和标准高斯过程协方差矩阵K是相互独立的， 而T经过超球面分解与Φ一一对应，通过公式（7）和（8）可分别得到参数θ和Φ。 

。

Images