具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

本发明提供了一种汽车复合材料保险杠防撞梁厚度的多目标优化方法，具体过程如下：

如图1所示，设定设计变量x₁代表防撞梁后板110的厚度，x₂代表防撞梁前板120的厚度，x₃代表防撞梁肋板130的厚度，x₄代表防撞梁下板140的厚度，x₅代表防撞梁上板150的厚度，定义上下板厚度相同；设计变量采用离散取值的方式，取值的间距为0.5。

因此铺层厚度设计变量的取值范围如下：

各厚度下铺层顺序如表1所示：

表1各铺层厚度对应的铺层顺序

定义优化数学模型为：

Minimize f_m(x),f_F(x)

Subjectto g_E(x)≥400J

g_D≤60.0mm

1.0≤x₁≤2.5

1.0≤x₂≤3

1.0≤x₃≤2.5

1.0≤x₄≤2。

数学模型构建完成后，依据目标函数以及约束函数进行多目标优化设计，采用高斯过程回归的第三代非支配排序遗传序列优化，主要包括如下步骤：

步骤1：在初始设计空间中对复杂系统进行最优拉丁超方试验(OLHS)生成初始样本点，并获得初始样本点处的目标函数f(x)和约束函数g(x)，如表2所示；

其中，f(x)＝[f₁(x),f₂(x),...f_j(x)...,f_n(x)]；

g(x)＝[g₁(x),g₂(x),...g_j(x)...,g_t(x)]；

式中，n、t分别为目标函数、约束函数的个数，x为设计变量。

表2初始样本表

依据对铺层结构划分，防撞梁设计变量共有5个，因为上下层合板采用相同的铺层方式，可以将设计变量缩减至4个，变量设计范围已经给出，根据层合板铺设的单层厚度，进行50组最优拉丁超立方抽样设计。

在有限元求解后获得目标函数f(x)和约束函数g(x)：

步骤2：根据目标函数和约束函数，设置高斯随机过程核函数初始参数；模型训练方法为极大似然估计，依据贝叶斯推断超参数θ的后验密度函数为：P(Y|X)与参数θ无关，所以该分布的对于一个给定观测值的关于参数的后验分布来说相当于是一个常数，并不会影响到后验分布的均值与方差常数，此时超参数的最大后验估计最大化P(Y|X,θ)P(θ),当先验P(θ)为均匀分布时，超参数的最大后验估计就是寻找使得P(Y|X,θ)最大的θ，即极大似然估计。对边际似然取对数作为似然函数((最大化)对数似然边际LML)，对超参数求偏导：

依据初始样本进行回归模型超参数调优以及模型训练，其中一个关键的设置为高维映射核函数的选择，单一核方法对于混杂着各种数值噪声的数据往往无法得到较为准确的回归模型。故设置初始化先验均值设置为0；meanfunc＝[]；核函数为复合核，具体包括：

平方指数(SE)函数：

马特温协方差(Matern)函数：

其中，为平滑系数，l为长度尺度参数，Γ(x)为欧拉第二积分，在实数域上的定义为：

为第二类贝塞尔函数，是下列微分方程的标准解函数y(x):

有理二次(RQ)函数：(l为长度尺度参数，α为比例混合系数)；

方差核函数选择有理二次和平方指数核函数，设置为

σ初始值为2，l₁初始值为0.25，l2初始值为2，α为10^-3。

并采用共轭梯度法对超参数进行优化；

式中，gp是回归函数，hyp，Hyper为核函数超参数，包括均值函数meanfunc，协方差函数covfunc。Minimize为共轭梯度优化函数，Sample_x为原设计点，此处即为设计变量[x₁,x₂,x₃,x₄]。Sample_f为样本点的目标函数，即数据集中的[m,F],要求最小化防撞梁的质量，最小化防撞梁碰撞力。Sample_g为约束函数；为数据集中的[D,E]；约束侵入量不得大于60mm，吸能量不得小于400J。

步骤3：生成参考点信息并初始化个体信息；进行模拟二进制交叉以及多项式变异生成新个体。

首先生成初始参考点，参考点生成方式如下：其中，M为目标函数的维度，H为分段数；

1)产生个维数为M-1的初始参考点集；此处维数为2。

X＝(1/H)*nchoosek(0:H+M-2,M-1)；

其中，nchoosek为抽样函数，第一个参数表示抽样总体，第二个参数表示抽样的个体数；

2)对于进行如下换:a_ij＝x_ij-(j-1)/H,获得变换集A；

3)利用变换集A,参考点依照下式生成:

通常，因划分段数H以及维数M增大会使参考点集过大，故采用一种内外层机制：H1为边界分段数，H2为内层分段数；

H1,H2大小应满足N为初始生成个体总数，且内层数应小于等于边界分段数；边界层遵照上述3个过程生成参考点S1，内层首先生成参考点S(3个过程)，取映射系数为0.5,S2由下列公式产生,并合并S1,S2；

对于2目标而言，H1＝149，H2＝0。设置种群个体数初始化为150，设置迭代次数200，

用GPR模型获取目标函数和约束函数的预测值以及对应的方差，并构建适应度信息；

基于核方法将贝叶斯线性回归映射成高维空间的非线性回归，对于函数f(x)视为下式所示服从高斯分布的高斯随机过程：

P(f_*|X,y,X_*,σ²)～N(σ^-2φ(X_*)^TA^-1φy,φ(X_*)^TA^-1φ(X_*))

μ_*＝k(X_*,X)(K+σ²I)^-1y

σ_*＝k(X_*,X_*)-k(X_*,X)(K+σ²I)^-1k(X,X_*)

获取2中超参数，结合分布信息估计未知个体位置处的目标函数、约束函数信息以及模型的误差，并利用惩罚机制构建用于非支配比较的适应度函数；

式中，a₁₁,a₁₂,…,a_n1,…a_nt为惩罚系数，PopDec是初始化个体位置信息，sf(x),sg(x)为目标函数和约束函数的预测值的误差函数，F为适应度函数；依照上述步骤进行种群个体目标函数，约束函数预测值和方差估计：PopDec是经过遗传操作产生的个体位置信息，维度与变量相同。

[f1,sf1]＝gp(hyp1,'infGaussLik',[],covfunc2,@likGauss,Sample_x,Sample_f1,PopDec)；

[f2,sf2]＝gp(hyp2,'infGaussLik',[],covfunc2,@likGauss,Sample_x,Sample_f2,PopDec)；

[g1,sg1]＝gp(hyp3,'infGaussLik',[],covfunc2,@likGauss,Sample_x,Sample_g1,PopDec)；

[g2,sg2]＝gp(hyp4,'infGaussLik',[],covfunc2,@likGauss,Sample_x,Sample_g2,PopDec)；

f12＝[f1 f2]；

f12_interval＝[sf1 sf2]；

g12＝[-g1+400g2-60]；

g12_interval＝[sg1 sg2]；

z＝[f12 g12]；

f_penalty＝100*max(0,-g1+400)+100*max(0,g2-60)；

y＝f12+f_penalty；

y为构造的适应度函数，Samplex是原始训练数据中的设计变量，Sample_f1、Sample_f2、Sample_g1、Sample_g2对应质量、碰撞力、吸能量、侵入量响应。f₁为防撞梁质量模型预测值，g₁为防撞梁吸能模型预测值，f₂为防撞梁碰撞力模型预测值，g₂为防撞梁吸能模型预测值。Sf,Sg分别是目标函数和约束函数预测值方差。

根据遗传操作进行个体位置的产生：

ILower＝repmat(Lower,N,1)

IUpper＝repmat(Upper,N,1)

Indiv_Pos＝unifrnd(ILower,IUpper)

其中，Upper、Lower为变量上、下限，

lower＝[1 1 1 1]；upper＝[2.5 3 2.5 2]；

repmat为重复构造函数，N＝150；unifrnd为随机数函数，产生一个NxN的介于Lower与Upper的随机数矩阵，保证样本的随机性、多样性；之后，按照遗传操作流程，首先确定交配池大小，并根据交配池进行筛选，其具体实现方法为：

随机产生行数为2，列数为N的个体序号(可重复，序号介于1～N),根据粒子违反约束的函数值，该值总是小于0的，且值越小对于约束的满足度越高，先随机产生两组个体，依据违约值对两个体取最优：

Tool＝randi(N,M,N)；

对违约值进行升序排序，得到对应原来位置的索引：

[～,index1]＝Sort(Popcon，[],1)；

本例中，Popcon＝0；再对索引进行升序排序，得到index2，此时index2便包含原位置个体对应的大小排序序号:

[～,index2]＝sortrows(index1)

根据index2信息，可知交配池中个体的大小排序，因为序列号越小，满足约束的可能性越大，所以对比两个个体，将对应较小个体排序序号选出；对于子代个体通过模拟二进制交叉以及多项式变异进行；模拟二进制交叉，模拟遗传算法二进制编码中的单点交叉：

主要遵循交叉前后父代与子代的均值相等：

定义传播因子(Spread Factorβ)，可得c₁与c₂:

其中，传播因子β可由下式确定：

其中，o为分布系数；同样，对于多项式变异，遵循以下公式：

indiv_k'＝indiv_k+δ(Upper_k-Lower_k)

δ₁＝(indiv_k-Lower_k)/(Upper_k-Lower_k)

δ₁＝(indiv_k-Upper_k)/(Upper_k-Lower_k)

本例中，r为一个随机算子，o,η＝20，若随机算子小于等于0.5，则按照上面式子计算传播因子，变异系数。否则按照第二个式子计算传播因子，变异系数。

步骤4：对步骤3生成的个体(总数为2*N)，依照适应度评估，划分非支配等级，定义储备集S_t(N(S_t)≥N),定义F_k为等级为k的非支配解集；对储备解集S_t进行归一化，在非支配等级划分时采用一种高效非支配排序方法ENS(EffcientNon_dominatedSort)减少计算中时间和空间复杂度，定义MaxFNo为当前最高支配等级(数字越小支配等级越高),初始值设定为0；FrontNo为前沿集合，初始设置为INF，涉及的主要步骤为：

1)依照对所有个体的适应度第一维度进行降序排序，获得排序信息；

2)定义大循环结束条件为FrontNo中不为INF的总个体数超过N，每次循环开始对MaxFNo进行加一操作，记录满足截止条件时最大非支配等级；对于当前个体，将其与之前处于前沿解集的个体进行非支配比较，比较从第二维度开始进行，若个体被支配，即进行下个个体的非支配比较；若个体未被支配，个体位置对应的FrontNo为当前MaxFNo；

3)上述步骤2结束，根据排序信息以及MaxFNo，获取非支配等级小于MaxFNo的所有个体，定义为P_t+1,非支配等级等于MaxFNo的个体定义为F_m,S_t为P_t+1和F_m的并集；

对获取到的S_t进行如下所示的归一化操作：

4)理想点位置获取，对于最小值问题，Z_min＝min(f₁,f₂,…f_n)；

5)坐标系转换，以理想点作为原点，变化后函数目标f’＝f-Zmin对于多目标问题，需要获取对应各个目标轴的极值点以及在各个目标轴的截距，其获取机制为:

XA＝1

dir为方向向量，大小为nxn,主对角线为1，其余位置均为1e-6；X为n个极值点的集合，大小为nxn，A为目标截距集合，大小为nx1，1是全为1且大小为nx1的列向量。

如果矩阵X的秩小于n，那么这n个极限点就不能构成一个n维的超平面。甚至即使超平面能够建立，也可能在某些方向上得不到截距或某些截距a_i不满足在所有上述情形下，对每个i∈{1,2,...,n}，设置为S_t中的非支配解在目标f_i上的最大值。对于S_t中各个目标函数进行归一化处理：

步骤5：如图2所示，以三维为例，计算解空间中参考线(参考点与理想点之间的连线)与归一化目标之间的距离，依照下式生成距离矩阵D：

其中，q代表S_t中个体数量，r对应参考点的数量；S_t每个个体考虑距离每个参考点距离，选择具有最小距离的参考点进行依附；

步骤6：考虑前MaxFNo-1前沿等级的情况下，根据参考点被选择的个数，选择参考点被选择较少的作为最后一前沿种群个体选择的依据，进而实现选择。个体选择依照小生境保持算子，步骤5取得的S_t所有个体对应参考点信息，定义ρ_j为P_t+1依附在参考点j的个体总数目，进行如下判断：

1)参考点选择，J_min＝{j:argminρ_j}，若J_min存在多个，随机选择

2)进行F_m(最后前沿解集)的选择，存在两种情况：

首先，若F_m中并未发现存在个体与相对应，那么将在当前迭代步不在考虑，从剩余参考点位置按照步骤1筛选新的参考点；若F_m存在个体j与关联，同样存在两种情况：第一种情况为，对应最少数目个体的参考点值为0，表示对于之前非支配前沿中并未存在个体与之对应，为了保持个体的多样性，从F_m中(如果存在多个，依照步骤5计算各个个体距离该参考点对应参考线的距离，取距离最小的个体)选取个体加入P_t+1，ρ_j的数值随着加入个体而增大；第二种情况，对应最少数目个体的参考点值不为0，表示已经存在个体与参考点与之对应，从F_m中依附于该参考点的众多个体(如果存在多个)进行随机挑选加入；ρ_j的数值随着加入个体而增大；

上述过程重复执行，直至P_t+1中的个体数达到需求的个体数N；依据设置迭代次数，个体不断演变进化，得到近似Pareto近似前沿解；

步骤7：将所述近似Pareto前沿解集作为真实函数的真实Pareto前沿解集，并获取每个个体基于欧拉距离的多目标期望提高准则：

当前最优解是由多个多维点共同组成的近似Pareto前沿：

式中，S为近似Pareto解的情况，下标表示目标个数，上标表示近似解的个数；EIM_e(x)为基于欧拉距离的期望提高准则，EI表示期望提高矩阵，EI每一行表示对于处于第j个前沿位置所有目标函数的提高值，EI每一列表示对于所有前沿位置第j个目标函数的提高值；f_i ^j(x)为近似Pareto前沿解集中第j号个体的第i个目标函数值，s-f_i(x)为第i个目标函数的高斯过程回归模型对应的预测值以及误差函数值；φ(.),Φ(.)为标准正态分布的概率密度函数和累积分布函数；

选择EIM_e值最大的个体作为新的样本点加入原样本点中生成新的样本点集；将所述新的样本点集作为初始样本点，重复步骤1-7，直至满足迭代次数，输出最终的样本点集。

按照上述步骤1～7，设置优化结果如图3所示，图3为优化Pareto前沿面，优化后依据基于熵值法的Toposis准则(C_i越接近1，个体越好)选取前沿Pareto中一个目标解：

熵值法(确定权值)：若前沿pareto解集中含有m个个体，每个个体为n维度目标函数。如下式：

对上式进行按列归一化处理:

计算n个目标函数各自的熵权值：

Toposis方法：对PY每个元素进行判断，若属于极大型指标(越大越好)，不做处理，若属于极小型指标，进行指标正向化：

得到同向化处理矩阵:

最优方案为同向化处理后每一列的最大值最劣方案为每一列的最小值

计算每个个体距离最优方案以及最劣方案的接近程度：

[x₁,x₂,x₃,x₄,x₅]＝[1,1.5262,1.4610,1.1295,1.1295]

[x₁,x₂,x₃,x₄,x₅]＝[1.0,1.5,1.5,1.0,1.0]

[E,F,D,M]＝[435.5802,10.6882,59.9988,0.8668]

优化结果满足目标以及约束条件。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。