融合径向基函数与克里金模型的高维代理模型构建方法

1.本发明涉及机械结构输入输出近似关系的构建，尤其适用于输入变量数目在30～200维的高维问题。本发明可以为高维度可靠性分析或优化设计提供有力的计算工具。


背景技术：

2.目前，复杂机械结构的可靠性分析、优化设计需要多次调用结构的响应函数。而结构的响应函数往往需要调用商用有限元软件、多体动力学软件、计算流体力学软件等。但是，单次调用商用软件的计算机模拟往往需要数十分钟、数小时甚至数天时间。而直接完成可靠性分析或优化设计，则可能需要数百次甚至数千次计算机模拟，这对实际工程来说是无法接受的。为此，缓解该瓶颈的一个办法就是使用代理模型(例如kriging模型、径向基函数插值、多项式回归、支持向量机回归、神经网络模型等)来构件输入变量与输出响应量的近似关系，建立高精度的机械结构输入-输出近似关系，对于结构的可靠性分析、优化设计具有重要意义。
3.近年来，kriging模型以其预测精度高且可以提供预测误差，在结构可靠性分析和优化设计领域应用越来越广泛。在结构优化设计领域，根据kriging模型的预测信息，可以选择最可能位于全局最优解附近的训练点，从而以尽量少的商用软件调用次数，完成最优解的找寻。在结构可靠性领域，根据kriging模型的预测信息，可以选择最可能位于极限状态函数附近的训练点，从而以尽量少的商用软件调用次数，完成失效概率的预测。kriging模型显著提高了结构可靠性分析和优化设计的效率。
4.但是，在完成高维度问题的函数近似时，kriging模型遭遇了“维度灾难”：训练kriging模型的计算时间大幅增长，达到或超过调用数值仿真模型的时间。建立高精度kriging模型的关键步骤在于对其似然函数进行全局优化，获取自相关参数θ的全局最优解，即最大似然估计(mle)。关于mle，需要说明的是：(1)一般情况下，自相关参数的个数等于随机变量的维度d；(2)计算似然函数，需要求解训练样本协方差矩阵的逆，其计算复杂度为o(n3)(n是训练点数目)；(3)似然函数是一个复杂的多模函数，存在多个局部最优解，因此，似然函数的寻优往往需要借助全局优化算法，例如遗传算法、差分进化算法等。上述因素导致，在处理高维问题时，自相关参数数目的增多，意味着全局优化搜寻空间的扩张以及似然函数计算次数的增加，从而导致建立kriging模型所需计算时间的增加。例如“专利公开号为：cn110008499a、名称为：一种基于bayesian kriging模型的机构优化设计方法”的中国发明专利，再如“专利公开号为：cn106202694a、名称为：基于组合预测方法的组合kriging模型构建方法”的中国发明专利，这些方法中kriging模型仍然存在大量待估计的自相关参数，因此，这些方法均不适用于高维问题。
5.为了减少自相关参数的个数，现有技术中已经有人提出了一种各项同性核函数，该核函数赋予输入变量相同的自相关参数。根据偏最小二乘法获取的主成分，研究人员构造了一种低维核函数，从而提出了适用于高维问题的kpls模型。其中，在构造核函数时，参数θ对kriging模型精度的影响至关重要。研究表明，参数θ的大小与变量对目标函数的影响
程度(近似)成比例。但是，上述方法构造的核函数，难以准确表征各个变量对目标函数的影响程度，kriging模型的预测精度也难以保证。为了构造更优质的核函数，研究人员对目标函数进行了相关性分析，获取了各个变量的影响系数，并将影响系数与缩放因子的乘积作为核心要素构造核函数。通过不断调整缩放因子，来使影响系数逼近最优的参数θ。考虑到sobol全局灵敏度分析的计算复杂性，研究人员分别引入了最大信息系数法(mic)和距离相关系数法(dic)来进行输入-输出相关性分析，取得了较好效果。两种方法分别记为kmic和kdic。
6.kriging模型由回归模型和高斯过程模型构成，其回归模型称为趋势函数(trend function)。目前，kriging模型建立时，趋势函数常常取为常数，趋势函数的选取问题常常被忽略。事实上，趋势函数对kriging模型的预测精度有较大影响，进而影响主动学习过程的收敛速度。早期，趋势函数一般从简单多项式中选取。趋势函数为常数的kriging模型称为ordinary kriging(ok)，为简单多项式的kriging模型称为universal kriging(uk)。为了选取最优的基函数，国内外多名学者提出了贝叶斯策略、优化策略、最小角回归策略等。但是，有关最优基函数选取的研究仍然集中在处理低维问题领域。
7.近年来，新的更加先进的回归技术被应用至kriging模型的趋势函数构造中。已有研究人员采用多项式混沌展开(pce)、多项式相关函数展开(pcfe)构造趋势函数。该文献采用径向基函数(rbf)构造趋势函数，从而提出了rbf-kriging模型。但是该技术存在以下问题：一方面，rbf-kriging是为低维优化问题提出的，将rbf-kriging拓展至高维问题的研究尚未出现。另一方面，rbf模型的精度常常受到基函数的形状参数的影响，rbf-kriging的精度也会相应受到影响。在单纯的rbf模型中，该形状参数可以通过交叉验证技术获取最优值。但是，在rbf与kriging的混合模型中，rbf的形状参数与kriging模型的超参数相互影响，如何根据当前训练点对rbf与kriging开展一体化参数校正，仍是需要攻克的问题。
8.因此，急需要提出一种融合径向基函数的、缩短高维问题训练时间的高精度kriging模型构建方法。


技术实现要素：

9.针对上述问题，本发明的目的在于提供融合径向基函数与克里金模型的高维代理模型构建方法，本发明采用的技术方案如下：
10.融合径向基函数与克里金模型的高维代理模型构建方法，其包括以下步骤：
11.采用拉丁超立方抽样法在原始的高维变量空间中抽取训练样本；
12.根据训练样本求得目标函数的真实值，并获得实验设计doe；
13.根据当前的实验设计doe，采用距离相关系数法求得变量的影响系数；
14.引入缩放因子，并构建单变核函数；
15.从训练样本中选取基函数参考点，并构建径向基函数；
16.将径向基函数对应的回归模型作为趋势函数，并将单变核函数作为高斯过程核函数，搭建得到kriging模型；
17.采用参数一体化校正法对高斯过程核函数对应的高斯过程超参数和径向基函数对应的径向基形状参数进行估计；
18.采用后验高斯过程模型进行预测，得到近似模型；
19.在原始的高维变量空间中抽取验证点，并进行近似模型的预测误差验证，得到高维代理模型。
20.进一步地，采用拉丁超立方抽样法在原始的高维变量空间中抽取训练样本，并组成初始训练点集x
(j)
(j＝1,
…
,n)；n表示训练样本的数量；
21.将训练样本代入目标函数g＝g(x)计算真实响应值g(x
(j)
)，形成当前doe[x
(j)
|g(x
(j)
)](j＝1,
…
,n)；其中，x＝(x1,
…
,xd)为d维输入变量。
[0022]
进一步地，根据当前的实验设计doe，采用距离相关系数法求得变量的影响系数，记为[w1,
…
,wd]，其中，wk表示第k个变量的影响系数。
[0023]
进一步地，所述单变核函数包括单变高斯核函数和单变样条线核函数；所述单变高斯核函数rh(x,x
′
)的表达式为：
[0024][0025]
其中，x＝(x1,
…
,xd)和x
′
＝(x1′
,
…
,x
′d)表示任意两个d维输入变量组成的向量；wk表示第k个变量的影响系数，λ表示引入的缩放因子，d表示维度；
[0026]
所述单变样条线核函数rh(x,x
′
)的表达为：
[0027][0028]
其中，s(m)表示关于变量m的样条函数，其表达式为：
[0029][0030]
进一步地，将径向基函数对应的回归模型作为趋势函数，并将单变核函数作为高斯过程核函数，搭建得到kriging模型；所述kriging模型的表达式为：
[0031][0032]
其中，表示径向基函数；c表示径向基函数的形状参数；zi(i＝1,
…
,nc)表示径向基函数的参考点，其从当前的实验设计中选取；β表示回归系数；∈(x)是高斯过程，且有e[(x),(x')]＝σ2rh(x,x
′
)，σ2是高斯过程方差。
[0033]
进一步地，采用参数一体化校正法对高斯过程核函数对应的高斯过程超参数和径向基函数对应的径向基形状参数进行估计，包括对参数{β,σ2,λ,c}进行估计。
[0034]
进一步地，采用最大似然估计法对参数进行估计，其表达为：
[0035][0036]
其中，n表示训练点数目；g
t
表示训练点位置功能函数值构成的向量；φ(c)表示训练点位置rbf值构成的矩阵；rh表示训练点之间的相关函数构成的矩阵，其由单变高斯核函
数rh(x,x
′
)求得；
[0037]
对公式(5)求偏导数，得到：
[0038][0039]
其中，表示回归系数β的估计值；表示高斯过程方差σ2的估计值；将公式(6)代入公式(5)，得到：
[0040][0041]
采用序列二次规划优化算法对l(λ,c)中的(λ,c)寻优，记其最优解为将最优解代入公式，即可回归系数β和高斯过程方差参数σ2的估计值。
[0042]
进一步地，，采用后验高斯过程模型进行预测，得到近似模型；
[0043]
预测均值μg(x)的表达式为：
[0044][0045]
其中r表示n
×
1的向量，其第i个元素ri＝rh(x,x(i))；表示形状参数为的径向基函数；
[0046]
预测方差的表达式为：
[0047][0048]
其中，t表示转置；表示形状参数为时训练点位置rbf值构成的矩阵。
[0049]
与现有技术相比，本发明具有以下有益效果：
[0050]
(1)本发明以径向基函数回归作为高维kriging模型的趋势函数，从而提高了高维kriging模型的预测精度。
[0051]
(2)本发明以单变核函数作为kriging模型的核函数，从而显著缩短高维问题中kriging模型的训练时间。
[0052]
(3)本发明采用参数一体化校正法对高斯过程核函数对应的高斯过程超参数和径向基函数对应的径向基形状参数进行估计，其好处在于，避免了径向基函数的形状参数取值不合适导致的精度下降。
[0053]
综上所述，本发明具有逻辑简单、有效解决高维问题等优点，在高维代理模型构建技术领域具有很高的实用价值和推广价值。
附图说明
[0054]
为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需使用的附图作简单介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对保护范围的限定，对于本领域技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
[0055]
图1为本发明的逻辑流程图。
具体实施方式
[0056]
为使本技术的目的、技术方案和优点更为清楚，下面结合附图和实施例对本发明作进一步说明，本发明的实施方式包括但不限于下列实施例。基于本技术中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本技术保护的范围。
[0057]
本实施例中，术语“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，a和/或b，可以表示：单独存在a，同时存在a和b，单独存在b这三种情况。
[0058]
本实施例的说明书和权利要求书中的术语“第一”和“第二”等是用于区别不同的对象，而不是用于描述对象的特定顺序。例如，第一目标对象和第二目标对象等是用于区别不同的目标对象，而不是用于描述目标对象的特定顺序。
[0059]
在本技术实施例中，“示例性的”或者“例如”等词用于表示作例子、例证或说明。本技术实施例中被描述为“示例性的”或者“例如”的任何实施例或设计方案不应被解释为比其它实施例或设计方案更优选或更具优势。确切而言，使用“示例性的”或者“例如”等词旨在以具体方式呈现相关概念。
[0060]
在本技术实施例的描述中，除非另有说明，“多个”的含义是指两个或两个以上。例如，多个处理单元是指两个或两个以上的处理单元；多个系统是指两个或两个以上的系统。
[0061]
如图1所示，本实施例提供了融合径向基函数与克里金模型的高维代理模型构建方法，本实施例高斯过程单变核函数构造和径向基趋势函数构造两部分进行介绍，以便理解本公开技术方法是如何将克里金模型与径向基函数联合起来并进行超参数一体化校正并对结构目标函数进行近似。具体步骤如下：
[0062]
第一步，采用在拉丁超立方抽样，在原始高维变量空间中抽取随机样本，并计算目标函数真实值。本方法在源空间中建立了kriging模型，而传统方法先通过降维技术，建立高维变量与低维隐变量之间的关系，而后在低维隐变量空间中建立kriging模型。其中，记x＝[x1,
…
,xd]是d维输入变量，g(x)为目标函数，在此，采用在拉丁超立方抽样，在原始高维变量空间中抽取随机样本，并计算目标函数真实值，记当前doe为[x
(j)
|g(x
(j)
)](j＝1,
…
,n)。
[0063]
第二步，根据当前doe，采用距离相关系数方法，获取变量的影响系数，记为[w1,
…
,wd]。该步骤的关键在于准确获取输入变量对输出变量的影响程度(或者称为灵敏度)。事实上，获取输入变量与输出变量的影响程度，一直是统计学领域研究的热点，常见的分析方法还包括sobol灵敏度、pearson相关系数等。sobol灵敏度的估计需要大量的样本，这与代理模型旨在尽量减少训练样本的目的是相悖的。pearson相关系数只适用于线性问题，对于变量之间相互影响、目标函数强非线性的问题，其容易给出错误的估计。dic能够以
较少的训练样本，获取较为准确的相关系数。因此，采用dic获取变量的相关系数。
[0064]
第三步，构造单变核函数。构造核函数时，自相关参数θ对kriging模型精度的影响至关重要，大量文献表明，参数θ的大小与变量对目标函数的影响程度(近似)成比例。因此，将dic获取的变量影响系数与缩放因子的乘积作为核心要素构造核函数，通过不断调整缩放因子，来使影响系数逼近最优的参数θ，该方案是合理的。kriging模型的核函数包括高斯核函数、样条线核函数。本发明采用的单变高斯核函数为
[0065][0066]
构造的单变样条线核函数为
[0067][0068]
且有
[0069]
第四步，从当前doe中选取nc≤n个训练点作为是参考点，记为当nc＝n时，rbf模型是插值模型；当nc＜n时，rbf模型为回归模型。本发明中，nc＝n/2，即一般的训练样本将作为参考点构造rbf。这一半的参考点可以通过k-means算法获取。由于拉丁超立方抽取的训练样本已经十分均匀，因此，本发明直接取训练样本的前面一半作为参考点。
[0070]
第五步，构造rbf回归模型的基函数
[0071][0072]
其中c是径向基函数的形状参数。径向基函数类型有多种，常见的有高斯径向基函数(gs-rbf)、薄板样条函数(th-rbf)、多二次函数(mq-rbf)等。详见表1。
[0073]
表1常见的径向基函数
[0074][0075]
第六步将rbf回归模型作为趋势函数，采用低维核函数高斯过程模型作为回归模型误差修正项，构造rbf-hd-kriging模型
[0076][0077]
且有e[(x),(x')]＝σ2rh(x,x
′
)。可以看到，rbf-hd-kriging是将趋势函数定义为径向基回归模型的高维kriging模型。趋势函数的精度对kriging模型整体的精度有重要影响。对于非线性功能函数，径向基回归模型的精度高于传统多项式函数，因此，将趋势函数定义为径向基函数，能够显著增强kriging模型对非线性目标函数的近似能力。
[0078]
第七步，采用参数一体化校正方法，对参数{β,σ2,λ,c}进行估计。
[0079]
(7.1)根据当前doe，推导似然函数
[0080][0081]
其中，g
t
为训练点位置功能函数值构成的向量，φ(c)为训练点位置rbf值构成的矩阵，rh为训练点之间的相关函数构成的矩阵(由核函数rh(x,x
′
)求得)。与传统kriging模型的似然函数相比，该似然函数与形状参数c有关。
[0082]
参数β的变量数目为nc。因此，上述似然函数的优化仍是高维问题。为避免求解该高维优化问题，特采用分阶段优化的策略：似然函数l(β,σ2,λ,c)要取得最大值，必有
[0083][0084]
进一步推导可得
[0085][0086]
(7.2)估计超参数(λ,c)。将上式代入似然函数，可得
[0087][0088]
(7.3)采用全局优化算法对(λ,c)寻优，代入和的方程即可获取β和σ2的估计值。这样就完成了rbf形状参数与高斯过程自相关参数的一体化校正。该优化问题有两个变量。这里采用序列二次规划算法，且搜寻的范围为[10-4
,100]2。
[0089]
第八步、利用下式对目标函数进行预测。获取了超参数的最优估计，即可开展预测。rbf-hd-kriging模型的预测均值为
[0090][0091]
其中r表示n
×
1的向量，其第i个元素ri＝rh(x,x(i))。
[0092]
预测方差的表达式为：
[0093][0094]
其中，
[0095]
第九步、开展rbf-hd-kriging模型的预测误差验证。近似模型精度评价指标常采用决定系数r2与标准化均方根误差nrmse来考核，其表达式分别为
[0096][0097][0098]
其中，yi和分别为第i(i＝1,
…
,n
val
)个验证点位置的响应值和预测值，为n
val
个验证点位置响应值的平均值。nrmse值越接近于0，r2值越接近于1，代表近似模型越精确。与训练点相同，验证点也采用拉丁超立方抽样方法获取。
[0099]
下面列举一实例：
[0100]
采用以下以若干算例，对本算法的有效性进行验证。如表2所示五个测试算例函数，其维度从20维值60维。训练点的数目为10倍维度，验证点数目均取5000。表3给出了采用所提方法rbf-hd-kriging所给出的测试结果，以及与其他方法的对比。本实施例采用的单变核函数为样条核函数，径向基函数为mq-rbf。可以看到，在相同条件下，所提方法的精度高于已有的方法，采用径向基函数作为趋势函数，并开展参数的一体化校正，对于提高高维度kriging模型的精度具有重要意义。
[0101]
表2测试算例函数
[0102][0103]
表3测试结果
[0104][0105]
上述实施例仅为本发明的优选实施例，并非对本发明保护范围的限制，但凡采用
本发明的设计原理，以及在此基础上进行非创造性劳动而作出的变化，均应属于本发明的保护范围之内。

技术特征：
1.融合径向基函数与克里金模型的高维代理模型构建方法，其特征在于，包括以下步骤：采用拉丁超立方抽样法在原始的高维变量空间中抽取训练样本；根据训练样本求得目标函数的真实值，并获得实验设计doe；根据当前的实验设计doe，采用距离相关系数法求得变量的影响系数；引入缩放因子，并构建单变核函数；从训练样本中选取基函数参考点，并构建径向基函数；将径向基函数对应的回归模型作为趋势函数，并将单变核函数作为高斯过程核函数，搭建得到kriging模型；采用参数一体化校正法对高斯过程核函数对应的高斯过程超参数和径向基函数对应的径向基形状参数进行估计；采用后验高斯过程模型进行预测，得到近似模型；在原始的高维变量空间中抽取验证点，并进行近似模型的预测误差验证，得到高维代理模型。2.根据权利要求1所述的融合径向基函数与克里金模型的高维代理模型构建方法，其特征在于，采用拉丁超立方抽样法在原始的高维变量空间中抽取训练样本，并组成初始训练点集x
(j)
(j＝1,
…
,n)；n表示训练样本的数量；将训练样本代入目标函数g＝g(x)计算真实响应值g(x
(j)
)，形成当前doe[x
(j)
|g(x
(j)
)](j＝1,
…
,n)；其中，x＝(x1,
…
,x
d
)为d维输入变量。3.根据权利要求2所述的融合径向基函数与克里金模型的高维代理模型构建方法，其特征在于，根据当前的实验设计doe，采用距离相关系数法求得变量的影响系数，记为[w1,
…
,w
d
]，其中，w
k
表示第k个变量的影响系数。4.根据权利要求3所述的融合径向基函数与克里金模型的高维代理模型构建方法，其特征在于，所述单变核函数包括单变高斯核函数和单变样条线核函数；所述单变高斯核函数r
h
(x,x
′
)的表达式为：其中，x＝(x1,
…
,x
d
)和x
′
＝(x
′1,
…
,x
′
d
)表示任意两个d维输入变量组成的向量；w
k
表示第k个变量的影响系数，λ表示引入的缩放因子，d表示维度；所述单变样条线核函数r
h
(x,x
′
)的表达为：其中，s(m)表示关于变量m的样条函数，其表达式为：5.根据权利要求4所述的融合径向基函数与克里金模型的高维代理模型构建方法，其
特征在于，将径向基函数对应的回归模型作为趋势函数，并将单变核函数作为高斯过程核函数，搭建得到kriging模型；所述kriging模型的表达式为：其中，表示径向基函数；c表示径向基函数的形状参数；z
i
(i＝1,
…
,n
c
)表示径向基函数的参考点，其从当前的实验设计中选取；β表示回归系数；∈(x)是高斯过程，且有e[(x),(x')]＝σ2r
h
(x,x
′
)，σ2是高斯过程方差。6.根据权利要求5所述的融合径向基函数与克里金模型的高维代理模型构建方法，其特征在于，采用参数一体化校正法对高斯过程核函数对应的高斯过程超参数和径向基函数对应的径向基形状参数进行估计，包括对参数{β,σ2,λ,c}进行估计。7.根据权利要求6所述的融合径向基函数与克里金模型的高维代理模型构建方法，其特征在于，采用最大似然估计法对参数进行估计，其表达为：其中，n表示训练点数目；g
t
表示训练点位置功能函数值构成的向量；φ(c)表示训练点位置rbf值构成的矩阵；r
h
表示训练点之间的相关函数构成的矩阵，其由单变高斯核函数r
h
(x,x
′
)求得；对公式(5)求偏导数，得到：其中，表示回归系数β的估计值；表示高斯过程方差σ2的估计值；将公式(6)代入公式(5)，得到：采用序列二次规划优化算法对l(λ,c)中的(λ,c)寻优，记其最优解为将最优解代入公式，即可回归系数β和高斯过程方差参数σ2的估计值。8.根据权利要求7所述的融合径向基函数与克里金模型的高维代理模型构建方法，其特征在于，采用后验高斯过程模型进行预测，得到近似模型；预测均值μ
g
(x)的表达式为：其中r表示n
×
1的向量，其第i个元素r
i
＝r
h
(x,x
(i)
)；表示形状参数为的径向基函数；
预测方差的表达式为：其中，t表示转置；表示形状参数为时训练点位置rbf值构成的矩阵。

技术总结
本发明公开了融合径向基函数与克里金模型的高维代理模型构建方法，包括：抽取训练样本；根据训练样本求得目标函数的真实值，并获得实验设计DoE；根据当前的实验设计DoE，采用距离相关系数法求得变量的影响系数；引入缩放因子，并构建单变核函数；从训练样本中选取基函数参考点，并构建径向基函数；将径向基函数对应的回归模型作为趋势函数，并将单变核函数作为高斯过程核函数，搭建得到Kriging模型；采用参数一体化校正法对高斯过程核函数对应的高斯过程超参数和径向基函数对应的径向基形状参数进行估计；采用后验高斯过程模型进行预测，得到近似模型；在原始的高维变量空间中抽取验证点，并进行近似模型的预测误差验证，得到高维代理模型。到高维代理模型。到高维代理模型。


技术研发人员：杨旭锋 汤炜 蒋文珂 邓杰 张懿
受保护的技术使用者：西南交通大学
技术研发日：2023.06.27
技术公布日：2023/9/9