CN105205262A - 一种基于二次响应面反演的随机模型修正方法 - Google Patents

Abstract

本发明提供一种基于二次响应面反演的随机模型修正方法，该方法针对工程系统参数识别或模型修正中存在的不确定性问题，通过单点逐次反演迭代，有效回避了常规随机模型修正方法中，因不确定性分析和灵敏度矩阵求解带来的误差。通过实例验证和方法对比说明，本方法可有效进行含明显不确定性结构的不确定参数识别，且识别效率和精度均很高。

Description

一种基于二次响应面反演的随机模型修正方法

技术领域

本发明涉及不确定性模型修正和系统参数识别领域，具体涉及一种基于二次响应面反演的随机模型修正方法。

背景技术

现行结构设计趋向精细化和大型化发展，以有限元分析为主的CAE分析技术在结构产品设计中扮演越来越重要的地位。在结构力学性能分析、失效预警、损伤识别等方面，准确而有效的有限元模型必不可少。然而由于结构本身的复杂性和建模过程中采取的不恰当简化策略等多方面原因，根据设计图纸建立的有限元模型并不能反应真实结构的力学性能。为了保证有限元模型和实验模型静动力响应一致性，可以利用模型修正技术进行有限元模型参数识别。

近年来，国内外学者在模型修正理论研究方面取得很大进展，发展了很多模型修正方法。但大多数模型修正方法都属于确定性模型修正方法，忽略了参数不确定性和测量响应不确定性。在实际问题中，不确定性是普遍存在的，如传感器配置或信号处理方面的不足；材料分散或设计制造误差；对复杂连接的不准确模拟等。确定性模型修正方法为寻找参数准确名义值努力，忽略了客观存在的不确定性，这极大限制了其在实际中的应用，可靠性大大降低。不确定性模型修正方法综合考虑系统存在的参数不确定性和响应测量不确定性，以对不确定参数进行可靠范围估计为目标，应用更广，尤其适用于含明显不确定性的复杂工程结构。

不确定性模型修正技术通过考虑实际结构存在的不确定性，并对不确定性量化分析，扩展了模型修正的适用范围，提高了修正可靠性。对于高频冲击、高度非线性大变形和随机耦合等存在明显不确定性的复杂结构，尤其适用。因此，开展相关理论研究并在实际结构上检验理论方法的有效性具有重大意义。

发明内容

本发明针对工程系统参数识别或模型修正中存在的不确定性问题，提供一种基于二次响应面反演的随机模型修正方法，该方法提高随机模型修正效率并保证了修正精度。该方法可有效进行含明显不确定性复杂结构不确定参数识别，且识别效率和精度均较高。

本发明采用如下技术方案实现：一种基于二次响应面反演的随机模型修正方法，其实施步骤如下：

步骤1：根据结构设计图纸建立初步有限元模型，并确定结构中存在不确定性的所有参数；结合实验模型，以结构响应输出(频率、位移等)构建优化目标函数，在不确定参数初始区间下进行优化修正，得到不确定参数名义值和与实验模型匹配度最高的初步修正模型；

步骤2：以步骤1中得到的不确定参数名义值为中心，构造包络涉及不确定性的初始超立方体；利用拉丁超立方抽样实验设计得到二次响应面构造所需原数据，并完成不确定参数初始超立方下的系统输入输出二次响应面函数拟合；

步骤3：根据步骤2中得到的二次响应面模型，利用最小二乘法建立基于二次响应面反演的不确定参数修正迭代算法，结合实验模型拉丁超立方抽样响应结果，进行不确定参数初始超立方修正；

步骤4：有效性判别：如果步骤3中所得修正后不确定参数超立方小于初始超立方体，则算法有效，转入下一步，否则，返回步骤2，重新确定初始超立方和二次响应面函数；

步骤5：计算迭代终止后不确定参数超立方中的二次响应面模型输出，并与实验模型响应拉丁超立方抽样结果对比，如果两者范围一致，则输出修正后不确定参数超立方体，否则，返回步骤2继续修正。

本发明的有益效果是在模型修正中综合考虑参数不确定性和响应测量不确定性，并利用非不确定量分析方法完成随机模型修正，提出一种基于二次响应面反演的随机模型修正方法。采用多项式响应面模拟系统输入输出关系并构建不确定参数反演迭代算法，不仅回避了运算复杂且精度较低的不确定性分析方法(概率方法、模糊方法和凸模型方法等)，提高了运算效率，同时，无需构建灵敏度矩阵，简化问题分析且回避了病态灵敏度矩阵问题，提高了修正精度。本发明可用于含明显不确定性的结构损伤识别、失效预警等方面，为结构静动力响应提供可靠估计，具有重要理论和现实意义。

附图说明

图1是本发明一种基于二次响应面反演的随机模型修正方法实现流程图；

图2是三自由度弹簧质量系统模型；

图3是三自由度弹簧质量系统一阶特征值二次响应面；

图4是三自由度弹簧质量系统二阶特征值二次响应面；

图5是三自由度弹簧质量系统三阶特征值二次响应面；

图6是三自由度弹簧质量系统修正前后eig1～eig2二维域分布；

图7是三自由度弹簧质量系统修正前后eig1～eig3二维域分布；

图8是三自由度弹簧质量系统修正前后eig1～|modal11|二维域分布；

图9是三自由度弹簧质量系统修正前后eig3～|modal11|二维域分布。

具体实施方式

下面结合附图以及具体实施例进一步说明本发明。

如图1所示，本发明为一种基于二次响应面反演的随机模型修正方法，包括以下步骤：

步骤1：根据结构设计图纸建立初步有限元模型，并确定结构中存在不确定性的所有参数。结合实验模型，以结构响应输出(频率、位移等)构建优化目标函数，在不确定参数初始区间下进行优化修正，得到不确定参数名义值和与实验模型匹配度最高的初步修正模型；

步骤2：以步骤1中得到的不确定参数名义值为中心，构造包络涉及不确定性的初始超立方体。利用拉丁超立方抽样实验设计得到二次响应面构造所需原数据，并完成不确定参数初始超立方下的系统输入输出二次响应面函数拟合；

步骤3：根据步骤2中得到的二次响应面模型，利用最小二乘法建立基于二次响应面反演的不确定参数修正迭代算法，结合实验模型拉丁超立方抽样响应结果，进行不确定参数初始超立方修正；

步骤4：有效性判别：如果步骤3中所得修正后不确定参数超立方小于初始超立方体，则算法有效，转入下一步，否则，返回步骤2，重新确定初始超立方和二次响应面函数；

步骤5：计算迭代终止后不确定参数超立方中的二次响应面模型输出，并与实验模型响应拉丁超立方抽样结果对比，如果两者范围一致，则输出修正后不确定参数超立方体，否则，返回步骤2继续修正。

在步骤1中，确定性模型修正包含以下步骤：

步骤1.1：根据结构设计图纸建立初步有限元模型，对结构中可能存在不确定性的参数或连接，如材料弹性模量、构件几何尺寸、复杂连接阻尼系数等，利用不确定参数变量标记并用区间变量描述其可能变化范围。

步骤1.2：以结构前n阶固有频和模态振型为系统输出，构建如下目标函数：

其中，w_fre,i、w_shap,ij分别为第i阶频率和第i阶振型的第j自由度位移权重系数；f_a,i(X)为第i阶频率理论计算值；f_m,i为第i阶频率实验测量值；为第i阶振型的第j自由度位移理论计算值；为第i阶振型的第j个自由度位移实验测量值。

步骤1.3：利用优化算法求取使得目标函数最小的不确定参数名义值。

在步骤2中，拟合二次响应面包含以下步骤：

步骤2.1：综合考虑结构可能涉及的不确定性范畴，以不确定性参数名义值为中心，构建不确定参数初始超立方体。

步骤2.2：在超立方体范围内，利用拉丁超立方实验设计确定实验方案，并计算结构响应输出，构建结构输入输出二次多项式响应面函数。

$\hat{Y}\left(X\right)=a+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_i{X}_i+\underset{i=1}{\overset{n-1}{\Σ}}\underset{j=i}{\overset{n}{\Σ}}{d}_{ij}{X}_i{X}_j+\mathrm{\ϵ}$

其中，a,b_i,d_ij为响应面系数，为响应面输出，X_i,X_j为输入参数，ε为残差。

在步骤3中，基于二次响应面反演进行不确定参数修正包含以下步骤：

步骤3.1：首先将上式写为如下形式：

$\mathrm{\ϵ}=\hat{Y}\left(X\right)-a-\underset{i=1}{\overset{n}{\Σ}}{b}_i{X}_i-\underset{i=1}{\overset{n-1}{\Σ}}\underset{j=i}{\overset{n}{\Σ}}{d}_{ij}{X}_i{X}_j$

其矩阵形式为： $\mathrm{\ϵ}=\hat{Y}-C-BX-DF\left(X\right)$

其中，ε＝(ε₁,ε₂,...,ε_K)^T为残差向量；为响应测量向量；C＝(a₁,a₂,...,a_K)^T为二次响应面常数向量；B＝[b_ki](k＝1,...,K；i＝1,...,n)为二次响应面一次项系数矩阵；为二次响应面二次项系数矩阵； $F\left(X\right)={(X_1^2,X_{2n}^2,...,X_n^2,X_1{X}_2,X_1{X}_3,...,X_{n-1}{X}_n)}^T$ 为待修正参数二次齐次矩阵；X＝(X₁,X₂,...,X_n)^T为待修正系数向量。

通过系统输出修正不确定参数区间转化为如下优化问题：

$\underset{X}{min}\left({\mathrm{\ϵ}}^T\mathrm{\ϵ}\right)$

$s.t.\underset{\‾}{X}\≤X\≤\stackrel{\‾}{X}$

根据最小二乘法，

$\begin{array}{c}{\mathrm{\ϵ}}^T\mathrm{\ϵ}={\left(\hat{Y}-C-BX-DF\left(X\right)\right)}^T\left(\hat{Y}-C-BX-DF\left(X\right)\right)\\ ={\left(\hat{Y}-C\right)}^T\left(\hat{Y}-C\right)-{\left(\hat{Y}-C\right)}^{T}BX-{\left(\hat{Y}-C\right)}^{T}DF\left(X\right)-X^T{B}^T\left(\hat{Y}-C\right)+\\ X^T{B}^{T}BX+X^T{B}^{T}DF\left(X\right)-F^T\left(X\right)D^T\left(\hat{Y}-C\right)+F^T\left(X\right)D^T\left(X\right)D^{T}BX+\\ F^T\left(X\right)D^{T}DF\left(X\right)\end{array}$

$\begin{array}{c}\frac{\∂\left({\mathrm{\ϵ}}^T\mathrm{\ϵ}\right)}{\∂X}=-B^T\left(\hat{Y}-C\right)-\frac{\∂F^T\left(X\right)}{\∂X}{D}^T\left(\hat{Y}-C\right)+B^{T}BX+B^{T}DF\left(X\right)+\\ \frac{\∂F^T\left(X\right)}{\∂X}{D}^{T}BX+\frac{\∂F^T\left(X\right)}{\∂X}{D}^{T}DF\left(X\right)=0\end{array}$

其中，

根据上式构造基于二次响应面反演的不确定参数修正迭代公式：

$X{|}_{l+1}={\left(B^{T}B+W+\frac{\∂F^T\left(X\right)}{\∂X}{D}^{T}B\right)}^{-1}{|}_l\{\left(B^T+\frac{\∂F^T\left(X\right)}{\∂X}{D}^T\right)\[\left(\hat{Y}-C\right)-DF\left(X\right)\]+WX\}$

其中，W＝λI为加权矩阵。

步骤3.2：将实验模型响应拉丁超立方抽样结果依次代入上式，进行不确定参数反演修正。

在步骤4中，算法有效性判断包含以下步骤：

将修正后不确定参数超立方体与初始超立方体比较。如果修正后超立方体大于初始超立方体，则退回步骤2，重新构造不确定参数初始超立方体及二次响应面；如果修正后参数超立方体小于不确定参数初始超立方体，则进行下一步。

在步骤5中，修正结果有效性判断包含以下步骤：

判断迭代收敛所得参数超立方体范围内，二次响应面输出范围与实验模型拉丁超立方抽样是否一致。如果两者一致，则输出最终修正所得不确定参数超立方体，即为修正结果；如果两者响应分布超出误差容限，则退回步骤2重新进行修正。

应用实例如下：

以图2所示三自由度弹簧质量系统说明本方法。系统参数取值为：

m₁＝1kg,m₂＝4kg,m₃＝1kg

，

k₁＝k₃＝0,k₆＝1N/m

k₂,k₄,k₅为不确定参数，且处于如下区间：

k₂＝[7.58.5]N/m,k₄＝[1.82.2]N/m,k₅＝[1.82.2]N/m

不确定参数初始区间假设如下：

k₂＝[6.59.5]N/m,k₄＝[1.62.4]N/m,k₅＝[1.62.4]N/m

系统输出定义为结构前三阶特征值(固有频率平方)，即：

z＝[eig1eig2eig3]^T。

其中，z为系统输出，eig1为第一阶特征值，eig2为第二阶特征值，eig3为第三阶特征值。

1)系统输出二次响应面拟合：

在参数初始区间内，二次响应面待定形式如下：

$\hat{y}\left(x\right)=a+\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i+\underset{i=1}{\overset{3}{\Σ}}\underset{j=i}{\overset{3}{\Σ}}{d}_{ij}{x}_i{x}_j$

根据最小二乘法，得到响应面计算结果如下：

一阶特征值：

$\begin{array}{c}eig1=-0.005k_2^2-0.0258k_4^2-0.0196k_5^2+0.0178k_2{k}_4+0.0278k_2{k}_5-\\ 0.0268k_4{k}_5+0.0813k_2+0.0800k_4+0.0023k_5+0.0951\end{array}$

二阶特征值：

$\begin{array}{c}eig2=0.0111k_2^2-0.8887k_4^2-1.3331k_5^2+0.0250k_2{k}_4+0.0244k_2{k}_5+\\ 2.5825k_4{k}_5-0.1488k_2-1.1498k_4+0.5179k_5+3.3411\end{array}$

三阶特征值：

$\begin{array}{c}eig3=-0.0062k_2^2+0.9145k_4^2+1.3827k_5^2-0.0428k_2{k}_4-0.0522k_2{k}_5-\\ 2.5557k_4{k}_5+0.3175k_2+2.3198k_4+0.4798k_5-1.1862\end{array}$

2)不确定反向传播模型修正：

基于二次响应面推导不确定参数修正迭代公式和拉丁超立方抽样测量响应，修正三自由度弹簧质量系统不确定参数。修正结果列于表1。

表1三自由度弹簧质量系统不同方法修正结果

与基于Kriging模型和径向基模型修正结果对比，本方法所得结果区间与真实区间吻合度最高。

Claims (6)

1.一种基于二次响应面反演的随机模型修正方法，其特征在于：包括以下步骤：

步骤1：根据结构设计图纸建立初步有限元模型，并确定结构中存在不确定性的所有参数，结合实验模型，以结构响应输出构建优化目标函数，在不确定参数初始区间下进行优化修正，得到不确定参数名义值和与实验模型匹配度最高的初步修正模型；

步骤2：以步骤1中得到的不确定参数名义值为中心，构造包络涉及不确定性的初始超立方体；利用拉丁超立方抽样实验设计得到二次响应面构造所需原数据，并完成不确定参数初始超立方下的系统输入输出二次响应面函数拟合；

步骤3：根据步骤2中得到的二次响应面模型，利用最小二乘法建立基于二次响应面反演的不确定参数修正迭代算法，结合实验模型拉丁超立方抽样响应测量结果，进行不确定参数初始超立方修正；

步骤4：有效性判别：如果步骤3中所得修正后不确定参数超立方小于初始超立方体，则算法有效，转入下一步，否则，返回步骤2，重新确定初始超立方和二次响应面函数；

步骤5：计算迭代终止后不确定参数超立方中的二次响应面模型输出，并与实验模型响应拉丁超立方抽样结果对比，如果两者范围一致，则输出修正后不确定参数超立方体，否则，返回步骤2继续修正。

2.根据权利要求1所述的基于二次响应面反演的随机模型修正方法，其特征在于：在步骤1中，包含以下步骤：

步骤1.1：根据结构设计图纸建立初步有限元模型，对结构中可能存在不确定性的参数或连接，利用不确定参数变量标记并用区间变量描述其可能变化范围；

步骤1.2：以结构前n阶固有频和模态振型为系统输出，构建如下目标函数：

其中，w_fre,i、w_shap,ij分别为第i阶频率和第i阶振型的第j自由度位移权重系数；f_a,i(X)为第i阶频率理论计算值；f_m,i为第i阶频率实验测量值；为第i阶振型的第j自由度位移理论计算值；为第i阶振型的第j个自由度位移实验测量值；

步骤1.3：利用优化算法求取使得目标函数最小的不确定参数名义值。

3.根据权利要求2所述的基于二次响应面反演的随机模型修正方法，其特征在于，在步骤2中，包含以下步骤：

步骤2.1：综合考虑结构可能涉及的不确定性范畴，以不确定性参数名义值为中心，构建不确定参数初始超立方体；

步骤2.2：在超立方体范围内，利用拉丁超立方实验设计确定实验方案，并计算结构响应输出，构建结构输入输出二次多项式响应面函数：

$\hat{Y}\left(X\right)=a+\underset{i=1}{\overset{n}{\Σ}}{b}_i{X}_i+\underset{i=1}{\overset{n-1}{\Σ}}\underset{j=i}{\overset{n}{\Σ}}{d}_{ij}{X}_i{X}_j+\mathrm{\ϵ}.$

其中，a,b_i,d_ij为响应面系数，为响应面输出，X_i,X_j为输入参数，ε为残差。

4.根据权利要求3所述的基于二次响应面反演的随机模型修正方法，其特征在于，在步骤3中，包含以下步骤：

步骤3.1：首先将步骤2.2中式写为如下形式：

$\mathrm{\ϵ}=\hat{Y}\left(X\right)-a-\underset{i=1}{\overset{n}{\Σ}}{b}_i{X}_i-\underset{i=1}{\overset{n-1}{\Σ}}\underset{j=i}{\overset{n}{\Σ}}{d}_{ij}{X}_i{X}_j$

其矩阵形式为：

$\mathrm{\ϵ}=\hat{Y}-C-BX-DF\left(X\right)$

其中，ε＝(ε₁,ε₂,...,ε_K)^T为残差向量；为响应测量向量；C＝(a₁,a₂,...,a_K)^T为二次响应面常数向量；B＝[b_ki](k＝1,...,K；i＝1,...,n)为二次响应面一次项系数矩阵；为二次响应面二次项系数矩阵； $F\left(X\right)={(X_1^2,X_2^2,...,X_n^2,X_1{X}_2,X_1{X}_3,...,X_{n-1}{X}_n)}^T$ 为待修正参数二次齐次矩阵；X＝(X₁,X₂,...,X_n)^T为待修正系数向量；

通过系统输出修正不确定参数区间转化为如下优化问题：

$\begin{array}{c}\underset{x}{\min }\left({\mathrm{\ϵ}}^T\mathrm{\ϵ}\right)\\ s.t.\underset{\‾}{X}\≤X\≤\stackrel{\‾}{X}\end{array}$

根据最小二乘法，

$\begin{array}{c}{\mathrm{\ϵ}}^T\mathrm{\ϵ}={(\hat{Y}-C-BX-DF(X\left)\right)}^T(\hat{Y}-C-BX-DF(X))\\ ={(\hat{Y}-C)}^T(\hat{Y}-C)-{(\hat{Y}-C)}^{T}BX-{(\hat{Y}-C)}^{T}DF\left(X\right)-X^T{B}^T(\hat{Y}-C)+\\ X^T{B}^{T}BX+X^T{B}^{T}DF\left(X\right)-F^T\left(X\right)D^T(\hat{Y}-C)+F^T\left(X\right)D^{T}BX+\\ F^T\left(X\right)D^{T}DF\left(X\right)\end{array}$

$\begin{array}{c}\frac{\∂\left({\mathrm{\ϵ}}^T\mathrm{\ϵ}\right)}{\∂X}=-B^T(\hat{Y}-C)-\frac{\∂F^T\left(X\right)}{\∂X}{D}^T(\hat{Y}-C)+B^{T}BX+B^{T}DF\left(X\right)+\\ \frac{\∂F^T\left(X\right)}{\∂X}{D}^{T}BX+\frac{\∂F^T\left(X\right)}{\∂X}{D}^{T}DF\left(X\right)=0\end{array}$

其中，

根据上式构造基于二次响应面反演的不确定参数修正迭代公式：

$X{|}_{l+1}={(B^{T}B+W+\frac{\∂F^T\left(X\right)}{\∂X}{D}^{T}B)}^{-1}{|}_l\{(B^T+\frac{\∂F^T\left(X\right)}{\∂X}{D}^T)\[(\hat{Y}-C)-DF\left(X\right)\]+WX\}$

其中，W＝λI为加权矩阵；

步骤3.2：将实验模型响应拉丁超立方抽样结果依次代入上式，进行不确定参数反演修正。

5.根据权利要求4所述的基于二次响应面反演的随机模型修正方法，其特征在于，在步骤4中，操作如下：

将修正后不确定参数超立方体与初始超立方体比较，如果修正后超立方体大于初始超立方体，则退回步骤2，重新构造不确定参数初始超立方体及二次响应面；如果修正后参数超立方体小于不确定参数初始超立方体，则进行下一步。

6.根据权利要求5所述的基于二次响应面反演的随机模型修正方法，其特征在于，在步骤5中，操作如下：

判断迭代收敛所得参数超立方体范围内，二次响应面输出范围与实验模型拉丁超立方抽样是否一致；如果两者一致，则输出最终修正所得不确定参数超立方体，即为修正结果；如果两者响应分布超出误差容限，则退回步骤2重新进行修正。