基于区间响应面模型的随机模型修正方法

技术领域

本发明涉及有限元模型修正技术领域，特别是一种基于区间响应面模型的随机模型修正方法。 

背景技术

有限元模型修正技术已广泛应用于航空航天、土木、机械等工程领域。然而，传统的模型修正方法^[1-2]是建立在参数确定性的假设基础上的。实际工程结构中的不确定性因素却是普遍存在且不可避免的，比如材料离散性、测量误差、加工制造误差等从而导致结构参数的不确定性。若仍然采用确定性的理论和方法对有限元模型进行修正，必然导致修正结果的不可靠，与实际情况有较大出入，因此最近几年考虑了参数不确定性的随机模型修正方法开始得到关注。 

当前，概率方法、模糊方法和区间分析方法是结构不确定性的量化及传播分析的三种主要方法^[3]。其中，最常用的是把结构参数视为随机变量来对问题进行建模和分析，即采用概率方法。此时参数的联合概率密度分布函数应该是已知的。但是，当没有足够多的数据来验证上述分布函数的正确性时，概率方法的精度就难以得到保证^[4]。而采用模糊方法^[5]来描述不确定性时，需要参数的隶属度函数。但在很多情况下，确定隶属度函数甚至比概率分布函数更为困难，即研究人员往往不得不带有很大的主观性来选取相应的隶属函数，使得分析结 果的可靠性也值得怀疑。因此，为了反映客观实际、减少主观因素的干扰，在结构测试信息不够充分的情况下，将工程问题中的不确定性参数视为有界的区间变量，采用区间分析方法^[6-10]来修正有限元模型是十分必要的。 

目前，基于区间分析的反演问题可以采用三种方法：区间算法^[6]，顶点法^[7]和全局优化方法^[4,8]。其中，区间算法定义了一系列的运算（区间加、减、乘、除法），可以用来计算结构响应的区间。然而，区间运算无法考虑变量之间的相关性，直接采用区间运算求解结构响应时很容易导致响应区间的严重扩张，而且这种扩张的程度往往难以量化，使其难以应用于工程实际。顶点法的应用需要满足三大前提条件：（1）结构的整体质量矩阵和刚度矩阵为修正参数的线性函数；（2）上述矩阵可以分解为非负定的子结构质量和刚度矩阵；（3）输出响应为系统特征值。对于实际的工程结构而言，上述条件往往难以同时满足；同时当输出响应包含系统的特征向量，或结构的整体质量和刚度矩阵为修正参数的非线性函数时，顶点法容易失效。文献[9]基于全局优化理论，运用Kriging模型替代有限元模型，采用确定性的模型修正技术，优化得到与每组实验数据相对应的结构参数，进而得到结构参数的区间。上述过程实际是多次的确定性模型修正过程，具有一定的局限性：（1）该方法需要具有全局搜索能力的优化算法，而此类算法比较复杂，优化结果容易陷入局部最优值，由此得到的结构参数区间不一定准确；（2）当试验得到的统计数据较多时，该方法需要大量的重复计算，将耗费大量的计算成本；（3）对于实际的工程结构，研 究人员往往更关注结构响应的最大、最小值，即响应区间的上界和下界，因此在仅知道响应区间的情况下，该方法很可能会失效。 

发明内容

本发明的目的是提供一种基于区间响应面模型的随机模型修正方法，本发明不仅避免了区间运算的扩张性问题，同时兼顾了结构响应区间的快速计算，使得模型修正过程不必进行有限元分析计算和灵敏度矩阵的构建，从而节省了大量的计算时间和耗费，也尽可能避免了病态优化问题的出现。在没有大量统计数据的前提下，通过区间模型修正过程，实现了由结构响应区间辨识参数区间的逆过程，量化了结构参数的不确定性，使修正的有限元模型能够更准确地预测结构的静动力响应。 

本发明采用以下方式实现：一种基于区间响应面模型的随机模型修正方法，其特征在于包括以下步骤： 

步骤S01：基于试验设计和回归分析构建不包含交叉项的二阶多项式响应面模型； 

步骤S02：采用配方法将多项式响应面表达式转化为完全平方形式； 

步骤S03：将区间参数代入响应面表达式，使得确定性响应面模型变为区间响应面模型； 

步骤S04：在区间响应面模型上进行区间运算得到预测的结构响应区间，并结合实测响应区间建立目标函数； 

步骤S05：构建优化反演问题来识别参数的区间分布。 

在本发明一实施例中，所述步骤S01的实现方式如下： 

步骤S011：基于试验设计方法中的中心复合设计法寻求拟合响应面所需的设计点，用于建立不包含交叉项的二次多项式响应面： 

$y={\beta }_0+\underset{i=1}{\overset{m}{\Sigma }}{\beta }_i{x}_i+\underset{i=1}{\overset{m}{\Sigma }}{\beta }_{\mathrm{ii}}{x}_i^2+ϵ---\left(1\right)$

式中：y代表结构响应，x_i表示第i个参数，m为参数个数，β₀,β_i,β_ii分别为常数项、一次项和二次项的待定系数，ε为拟合误差； 

步骤S012：在确定设计点后，通过数值计算获取每个设计点所对应的响应值，以得到一个样本； 

步骤S013：基于所有样本，通过最小二乘回归分析可估计多项式的待定系数；拟合后的响应面模型表示为： 

$y={\hat{\beta }}_0+\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_i{x}_i+\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_{\mathrm{ii}}{x}_i^2---\left(2\right)$

式中为常数项、一次项、二次项的回归系数。 

在本发明一实施例中，所述步骤S02和S03的实现方式如下： 

在得到如式（2）所示的二次多项式响应面模型以后，采用配方法，将响应面模型转化成完全平方形式，使设计参数在函数表达式中仅出现一次，转化后的形式如下所示： 

$y={\hat{\beta }}_0+\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_i{x}_i+\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_{\mathrm{ii}}{x}_i^2$           （3） 

$=\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_{\mathrm{ii}}{(x_i+\frac{{\hat{\beta }}_i}{2{\hat{\beta }}_{\mathrm{ii}}})}^2+{\hat{\beta }}_0-\underset{i=1}{\overset{m}{\Sigma }}\frac{{\hat{\beta }}_i^2}{4{\hat{\beta }}_{\mathrm{ii}}}$

而后将参数x_i采用区间数x^I表示，即得到区间响应面模型： 

$y^I=\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_{\mathrm{ii}}{(x_i^I+\frac{{\hat{\beta }}_i}{2{\hat{\beta }}_{\mathrm{ii}}})}^2+{\hat{\beta }}_0-\underset{i=1}{\overset{m}{\Sigma }}\frac{{\hat{\beta }}_i^2}{4{\hat{\beta }}_{\mathrm{ii}}}---\left(4\right).$

在本发明一实施例中，所述步骤S04和S05的实现方式如下： 

区间模型修正过程可以转化为以下的优化反问题： 

$\min \left[R\left(x^I\right)\right],R\left(x^I\right)=\underset{i=1}{\overset{P}{\Sigma }}{\left(\frac{\underset{‾}{y_{\mathrm{ia}}}-\underset{‾}{y_{\mathrm{ie}}}}{\underset{‾}{y_{\mathrm{ie}}}}\right)}^2+\underset{i=1}{\overset{P}{\Sigma }}{\left(\frac{\overline{y_{\mathrm{ia}}}-\overline{y_{\mathrm{ie}}}}{\overline{y_{\mathrm{ie}}}}\right)}^2---\left(5\right)$

st.VLB≤x≤VUB 

式中x^I为结构区间参数；和分别为结构第i阶理论和实测响应的下界、上界；P为结构响应的个数；st.表示约束条件；VLB，VUB是设计空间范围；R为残差,这里作为目标函数；然后通过寻找残差R的最小值，求解参数的区间分布。 

工程结构的材料、几何参数、边界条件和连接条件中不可避免地存在不确定性，因此在模型修正过程中考虑上述参数的不确定性是十分必要的，可以使模型修正方法更符合实际。对于复杂的工程结构，因为试验成本很高，所以通过大量的试验来获取足够多的统计数据以描述结构响应和参数的概率密度分布函数往往是不现实的。此时采用基于概率统计的随机模型修正方法可能会有较大的计算误差，甚至是不可行的。因此，在测量数据较少的情况下，基于区间分析的随机模型修正方法具有较好的优势。本发明的优点在于（1）忽略了交叉项（参数间相互效应）的区间响应面模型可以避免区间运算可能导致的区间扩张问题，从而保证参数区间的计算精度；（2）修正过程无需构建灵敏度矩阵，大大简化了优化问题并避免了病态收敛过程的出现；（3）修正过程直接在区间响应面模型上进行优化迭代，无需进行复杂的数值分析计算，大大提高了区间模型修正的效率；（4）可以适用 于参数不确定性（区间分布范围）程度较高的情况；（5）可以同时适用于试验数据较多和较少的情况；（6）模型修正过程减少了主观因素的干扰，所得结果更真实可靠。 

附图说明

图1是本发明方法流程示意图。 

图2是本发明实施例区间模型修正流程图。 

具体实施方式

下面结合附图及实施例对本发明做进一步说明。 

如图1所示，一种基于区间响应面模型的随机模型修正方法，其特征在于包括以下步骤： 

步骤S01：基于试验设计和回归分析构建不包含交叉项的二阶多项式响应面模型； 

步骤S02：采用配方法将多项式响应面表达式转化为完全平方形式； 

步骤S03：将区间参数代入响应面表达式，使得确定性响应面模型变为区间响应面模型； 

步骤S04：在区间响应面模型上进行区间运算得到预测的结构响应区间，并结合实测响应区间建立目标函数； 

步骤S05：构建优化反演问题来识别参数的区间分布。 

在区间运算中，参数变量是作为独立变量进行处理的，容易导致计算结果的区间扩张。从结构响应与参数的函数关系式来看，若参数在函数表达式中多次出现，区间算法很可能会导致响应区间的扩张。 因此，为了避免区间运算的扩张性，本发明提出了“区间响应面模型”的概念，即采用不包含交叉项的二次多项式响应面来近似表示结构响应与不确定性参数之间的关系，再通过配方法将响应面模型转化为完全平方形式，并将参数用区间数表示，从而得到区间响应面模型。由于结构参数在区间响应面模型中仅出现一次，因此，采用区间运算求解响应区间时将不会引起区间扩张。本发明不仅避免了区间运算的扩张性问题，同时兼顾了结构响应区间的快速计算，使得模型修正过程不必进行有限元分析计算和构建灵敏度矩阵，节省了大量的计算时间和耗费，并尽可能避免了病态优化问题的出现。在没有大量统计数据的前提下，通过区间模型修正过程，实现了由结构响应区间辨识参数区间的逆过程，量化了结构参数的不确定性，使修正的有限元模型能够更准确地预测结构的静动力响应。 

具体的，请参见图2，图2是本实施例区间模型修正流程图，本实施例所采用的技术方案主要包括二次多项式响应面的建立，区间响应模型的构造，响应的区间运算，以及基于区间响应面模型的随机模型修正过程。具体步骤如下： 

1、建立二次多项式响应面 

响应面模型^[11]实质上是一种多项式数学表达式，可以用于建立系统输入（参数）与输出（响应）之间的关系，该表达式通过拟合设计样本得到。本发明基于试验设计方法中的中心复合设计法寻求拟合响应面所需的设计点（design points），建立了不包含交叉项（相互效应）的二次多项式响应面： 

$y={\beta }_0+\underset{i=1}{\overset{m}{\Sigma }}{\beta }_i{x}_i+\underset{i=1}{\overset{m}{\Sigma }}{\beta }_{\mathrm{ii}}{x}_i^2+ϵ---\left(1\right)$

式中：y代表结构响应，x_i表示第i个参数，m为参数个数，β₀,β_i,β_ii分别为常数项、一次项和二次项的待定系数，ε为拟合误差； 

在确定设计点后，可以通过数值计算获取每个设计点所对应的响应值，以得到一个样本。然后基于所有样本，通过最小二乘回归分析可估计多项式的待定系数。拟合后的响应面模型可以表示为： 

$y={\hat{\beta }}_0+\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_i{x}_i+\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_{\mathrm{ii}}{x}_i^2---\left(2\right)$

式中为常数项、一次项、二次项的回归系数。 

2、区间响应面模型的构造 

在得到如式（2）所示的二次多项式响应面模型以后，即可采用配方法，将响应面模型转化成完全平方形式，使设计参数在函数表达式中仅出现一次，以尽可能避免区间运算过程的区间扩张问题。转化后的形式如下所示： 

$y={\hat{\beta }}_0+\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_i{x}_i+\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_{\mathrm{ii}}{x}_i^2$              （3） 

$=\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_{\mathrm{ii}}{(x_i+\frac{{\hat{\beta }}_i}{2{\hat{\beta }}_{\mathrm{ii}}})}^2+{\hat{\beta }}_0-\underset{i=1}{\overset{m}{\Sigma }}\frac{{\hat{\beta }}_i^2}{4{\hat{\beta }}_{\mathrm{ii}}}$

而后将参数x_i采用区间数x^I表示，即得到区间响应面模型： 

$y^I=\underset{i=1}{\overset{m}{\Sigma }}{\hat{\beta }}_{\mathrm{ii}}{(x_i^I+\frac{{\hat{\beta }}_i}{2{\hat{\beta }}_{\mathrm{ii}}})}^2+{\hat{\beta }}_0-\underset{i=1}{\overset{m}{\Sigma }}\frac{{\hat{\beta }}_i^2}{4{\hat{\beta }}_{\mathrm{ii}}}---\left(4\right)$

3、结构响应的区间运算 

得到区间响应面模型后，便可直接基于区间运算法则进行区间运算，得到结构响应的区间分布。 

4、区间模型修正过程 

最后，区间模型修正过程可以转化为以下的优化反问题： 

$\min \left[R\left(x^I\right)\right],R\left(x^I\right)=\underset{i=1}{\overset{P}{\Sigma }}{\left(\frac{\underset{‾}{y_{\mathrm{ia}}}-\underset{‾}{y_{\mathrm{ie}}}}{\underset{‾}{y_{\mathrm{ie}}}}\right)}^2+\underset{i=1}{\overset{P}{\Sigma }}{\left(\frac{\overline{y_{\mathrm{ia}}}-\overline{y_{\mathrm{ie}}}}{\overline{y_{\mathrm{ie}}}}\right)}^2---\left(5\right)$

st.VLB≤x≤VUB 

式中x^I为结构区间参数；和分别为结构第i阶理论和实测响应的下界、上界；P为结构响应的个数；st.表示约束条件；VLB，VUB是设计空间范围；R为残差，这里即目标函数。 

然后通过寻找残差R的最小值，求解参数的区间分布。 

以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。 

参考文献 

[1]M.I.Friswell,J.E.Mottershead,Finite Element Model Updating in Structural Dynamics,Kluwer Academic Press,Dordrecht,1995. 

[2]李辉,丁桦.结构动力模型修正方法研究进展[J].《力学进展》,2005,35(2):170‐180. 

[3]D.Moens,D.Vandepitte,A survey of non‐probabilistic uncertainty treatment in finite element analysis,Computer Methods in Applied Mechanics and Engineering194(12–16)(2005)1527–1555. 

[4]I.Elishakoff.Essay on uncertainties in elastic and viscoelastic structures:from A.M.Freudenthal’s criticisms to modern convex modeling[J].Computers and Structures,1995,56(6):871‐895. 

[5]Li Chen,S.S.Rao.Fuzzy finite‐element approach for the vibration analysis of imprecisely‐defined systems.Finite Elements in Analysis and Design27(1997)69‐83. 

[6]R.Moore,R.Kearfott,M.Cloud,Introduction to Interval Analysis,Society for Industrial and Applied Mathematics,Philadelphia,2009. 

[7]Qiu Zhi‐ping,Wang Xiao‐jun,M.I Friswell.Eigenvalue bounds of structures with uncertain‐but‐bounded parameters[J].Journal of Sound and Vibration,2005,282(1):297‐312. 

[8]E.Hansen and G.W.Walster,Global Optimization Using Interval Analysis.Marcel Dekker,New York,2003. 

[9]H.H.Khodaparast,J.E.Mottershead,K.J.Badcock.Interval model updating with irreducible uncertainty using the Kriging predictor,Mech.Syst.Sig.Process.25(2011),1204–1226. 

[10]Y.S Erdogan,P.G Bakir,Inverse propagation of uncertainties in finite element model updating through use of fuzzy arithmetic,Engineering applications of artificial intelligence  26(2013)357‐367. 

D.C.Montgomery,Design and analysis of experiments,3rd ed,New York:John Wiley&Sons,Inc;2006。