CN115183969A - Method and system for estimating BWBN model parameters - Google Patents

Method and system for estimating BWBN model parameters Download PDF

Info

Publication number
CN115183969A
CN115183969A CN202210770940.3A CN202210770940A CN115183969A CN 115183969 A CN115183969 A CN 115183969A CN 202210770940 A CN202210770940 A CN 202210770940A CN 115183969 A CN115183969 A CN 115183969A
Authority
CN
China
Prior art keywords
identified
parameter
bwbn
model
optimization
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210770940.3A
Other languages
Chinese (zh)
Inventor
赵露
丁孙玮
涂田刚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Institute of Materials
Original Assignee
Shanghai Institute of Materials
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Institute of Materials filed Critical Shanghai Institute of Materials
Priority to CN202210770940.3A priority Critical patent/CN115183969A/en
Publication of CN115183969A publication Critical patent/CN115183969A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M7/00—Vibration-testing of structures; Shock-testing of structures
    • G01M7/02—Vibration-testing by means of a shake table
    • G01M7/022—Vibration control arrangements, e.g. for generating random vibrations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及一种BWBN模型的参数估计方法及系统，方法包括：获取实验的测量数据和物理参数，将测量数据以及物理参数代入BWBN模型并求解，确定待辨识参数；对待辨识参数进行初始化，确定取值范围和初始值，假定各个待辨识参数满足高斯分布，分别计算似然函数确定各个待辨识参数的初始先验分布；分别自每个待辨识参数的先验分布中抽取样本并基于样本更新各个待辨识参数的估算值，重复此步骤，直至满足预设置的终止条件；输出每个待辨识参数的估算值。与现有技术相比，本发明以测量数据为驱动，基于iTMCMC采样方法估计BWBN非线性迟滞模型的参数，具有较高通用性，可应用于各类迟滞阻尼振动近似满足BWBN模型的辨识问题中。

Figure 202210770940

The invention relates to a method and system for estimating parameters of a BWBN model. The method includes: acquiring measurement data and physical parameters of an experiment, substituting the measurement data and physical parameters into a BWBN model and solving, and determining parameters to be identified; initializing the parameters to be identified, and determining Value range and initial value, assuming that each parameter to be identified satisfies the Gaussian distribution, calculate the likelihood function to determine the initial prior distribution of each parameter to be identified; draw samples from the prior distribution of each parameter to be identified and update based on the sample For the estimated value of each parameter to be identified, this step is repeated until the preset termination condition is met; the estimated value of each parameter to be identified is output. Compared with the prior art, the present invention is driven by measurement data and estimates the parameters of the BWBN nonlinear hysteresis model based on the iTMCMC sampling method. .

Figure 202210770940

Description

一种BWBN模型的参数估计方法及系统A method and system for parameter estimation of BWBN model

技术领域technical field

本发明涉及振动控制领域，尤其是涉及一种BWBN模型的参数估计方法及系统。The invention relates to the field of vibration control, in particular to a parameter estimation method and system of a BWBN model.

背景技术Background technique

随着现代工业化和社会的发展，抑制或至少衰减可能影响系统或结构的不良振动，对建筑桥梁的安全性、加工的精度、设备的工作稳定性、交通工具乘坐的舒适性有重大意义。振动控制通常使用被动、半主动或主动系统来实现，其中每个系统都有相当大的滞回性能。Bouc-Wen-Baber-Noori(BWBN)模型可有效用于分析描述多系统滞回性能，捕获一系列滞回循环模式。With the development of modern industrialization and society, suppressing or at least attenuating undesirable vibrations that may affect the system or structure is of great significance to the safety of building bridges, the accuracy of processing, the working stability of equipment, and the comfort of vehicles. Vibration control is typically accomplished using passive, semi-active, or active systems, each of which has considerable hysteresis. The Bouc-Wen-Baber-Noori (BWBN) model can be effectively used to analyze and describe the hysteretic performance of multiple systems, capturing a series of hysteretic loop modes.

Bouc-Wen-Baber-Noori滞回模型在经典Bouc-Wen滞回模型基础上进一步发展，考虑强度、刚度退化和捏拢效应。许多工程结构在动荷载作用下会进入非弹性状态而表现出滞回特性。滞回特性也称为弹塑性。滞回一般来自材料的非线性特性、接触面的摩擦特性和结合面之间的接触变形等。在荷载作用下结构的恢复力与位移之间存在滞回关系，当荷载具有周期性时，由于结构进入弹塑性变形，塑性变形的位移不可全部恢复，使得曲线沿不同路径变化，形成滞回环。结构的实际滞回曲线十分复杂，难以应用于分析结构的非线性特性，因而需要建立既便于数学描述又能反映结构滞回特性的滞回模型来预测工程结构非线性动力响应。The Bouc-Wen-Baber-Noori hysteresis model is further developed on the basis of the classical Bouc-Wen hysteresis model, considering strength, stiffness degradation and pinching effect. Many engineering structures will enter an inelastic state and exhibit hysteretic properties under dynamic loads. Hysteretic properties are also known as elastoplasticity. Hysteresis generally comes from the nonlinear characteristics of the material, the friction characteristics of the contact surface and the contact deformation between the joint surfaces. There is a hysteretic relationship between the restoring force and the displacement of the structure under the action of the load. When the load is periodic, because the structure enters the elastic-plastic deformation, the displacement of the plastic deformation cannot be fully recovered, so that the curve changes along different paths, forming a hysteretic loop. The actual hysteresis curve of the structure is very complex, and it is difficult to apply it to analyze the nonlinear characteristics of the structure. Therefore, it is necessary to establish a hysteresis model that is convenient for mathematical description and can reflect the hysteresis characteristics of the structure to predict the nonlinear dynamic response of engineering structures.

BWBN模型改进了Bouc-Wen模型无法描述捏缩滞回现象的问题，同时带来的问题是参数维度较大，需要BWBN模型描述系统的滞回性能，就需要确定模型中的多项参数。现有的参数估计方法均需要用到不同的滤波器，传统MCMC方法容易落入单峰陷阱，且面对高纬度参数估计时有失效风险。而TMCMC适用于不确定参数很少的问题，如果不确定参数的数量很大，结果可能会有很大偏差，同时统计估计的效率较低。The BWBN model improves the problem that the Bouc-Wen model cannot describe the pinch-shrink hysteresis phenomenon. At the same time, the problem is that the parameter dimension is large. If the BWBN model is needed to describe the hysteretic performance of the system, it is necessary to determine several parameters in the model. The existing parameter estimation methods all need to use different filters, the traditional MCMC method is easy to fall into the single-peak trap, and there is a risk of failure in the face of high-latitude parameter estimation. However, TMCMC is suitable for problems with few uncertain parameters. If the number of uncertain parameters is large, the results may be greatly biased, and the efficiency of statistical estimation is low.

发明内容SUMMARY OF THE INVENTION

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种BWBN模型的参数估计方法及系统。The purpose of the present invention is to provide a method and system for estimating parameters of a BWBN model in order to overcome the above-mentioned defects in the prior art.

本发明的目的可以通过以下技术方案来实现：The object of the present invention can be realized through the following technical solutions:

根据本发明的第一方面，提供了一种BWBN模型的参数估计方法，包括以下步骤：According to a first aspect of the present invention, a parameter estimation method of a BWBN model is provided, comprising the following steps:

获取实验的测量数据，所述测量数据包括力和位移，确定实验的物理参数，建立BWBN模型，将测量数据以及物理参数代入BWBN模型，求解BWBN模型，确定BWBN模型中的W个待辨识参数；Obtain the measurement data of the experiment, the measurement data includes force and displacement, determine the physical parameters of the experiment, establish a BWBN model, substitute the measurement data and the physical parameters into the BWBN model, solve the BWBN model, and determine the W parameters in the BWBN model to be identified;

分别对每个待辨识参数进行初始化，确定其最大值、最小值和初始值，假定各个待辨识参数满足高斯分布，分别计算各个待辨识参数的似然函数，得到各个待辨识参数的初始先验分布，将初始先验分布作为待辨识参数的先验分布，将初始值作为待辨识参数的估算值；Initialize each parameter to be identified, determine its maximum value, minimum value and initial value. Assuming that each parameter to be identified satisfies the Gaussian distribution, calculate the likelihood function of each parameter to be identified separately, and obtain the initial prior of each parameter to be identified. distribution, take the initial prior distribution as the prior distribution of the parameter to be identified, and take the initial value as the estimated value of the parameter to be identified;

分别自每个待辨识参数的先验分布中抽取样本并基于样本更新各个待辨识参数的估算值，完成一轮优化，重复此步骤，直至满足预设置的终止条件，停止迭代；Samples are drawn from the prior distribution of each parameter to be identified and the estimated value of each parameter to be identified is updated based on the sample, a round of optimization is completed, and this step is repeated until the preset termination condition is met, and the iteration is stopped;

输出每个待辨识参数的估算值。An estimate of each parameter to be identified is output.

进一步地，第j轮优化中自第s个待辨识参数的先验分布中抽取样本并基于样本更新第s个待辨识参数的估算值，其中，j＝1、2、…，s＝1、2、…、W，具体为：Further, in the j-th round of optimization, samples are drawn from the prior distribution of the s-th parameter to be identified, and the estimated value of the s-th parameter to be identified is updated based on the samples, where j=1, 2, ..., s=1, 2, ..., W, specifically:

从第s个待辨识参数的先验分布中抽取Nj,s个样本：

Figure BDA0003724195220000021
Nj,s为预设置的第j轮优化中第s个待辨识参数的样本抽取值；Draw N j,s samples from the prior distribution of the s-th parameter to be identified:
Figure BDA0003724195220000021
N j,s is the preset sample extraction value of the s-th parameter to be identified in the j-th round of optimization;

计算第j轮优化的渐变系数，计算公式如下：Calculate the gradient coefficient of the jth round of optimization, the calculation formula is as follows:

Figure BDA0003724195220000022
Figure BDA0003724195220000022

其中，

Figure BDA0003724195220000023
表示第j轮优化中第s个待辨识参数的样本变异系数；in,
Figure BDA0003724195220000023
Represents the sample variation coefficient of the sth parameter to be identified in the jth round of optimization;

分别计算每个样本的加权系数，计算公式如下：Calculate the weighting coefficient of each sample separately, the calculation formula is as follows:

Figure BDA0003724195220000024
Figure BDA0003724195220000024

Figure BDA0003724195220000025
表示第j轮优化中第s个待辨识参数的第k个样本的加权系数，D表示实验的测量数据；
Figure BDA0003724195220000025
represents the weighting coefficient of the k-th sample of the s-th parameter to be identified in the j-th round of optimization, and D represents the experimental measurement data;

计算Nj,s个样本的加权系数的平均值，计算公式如下：Calculate the average value of the weighting coefficients of N j, s samples, and the calculation formula is as follows:

Figure BDA0003724195220000031
Figure BDA0003724195220000031

其中，Sej,s表示第j轮优化中第s个待辨识参数的样本加权系数平均值；Among them, Se j,s represents the average value of the sample weighting coefficient of the s-th parameter to be identified in the j-th round of optimization;

计算协方差矩阵，计算公式如下：Calculate the covariance matrix, the calculation formula is as follows:

Figure BDA0003724195220000032
Figure BDA0003724195220000032

其中，∑j,s表示第j轮优化中第s个待辨识参数的Nj,s个样本的协方差矩阵，

Figure BDA0003724195220000033
表示第j轮优化中第s个待辨识参数的Nj,s个样本的均值，T表示矩阵转置，ξj表示预设置的第j轮优化中的自适应缩放因子；Among them, ∑ j,s represents the covariance matrix of N j,s samples of the sth parameter to be identified in the jth round of optimization,
Figure BDA0003724195220000033
represents the mean value of N j, s samples of the s-th parameter to be identified in the j-th round of optimization, T represents the matrix transposition, and ξ j represents the preset adaptive scaling factor in the j-th round of optimization;

从

Figure BDA0003724195220000034
中随机生成索引l，按照概率
Figure BDA0003724195220000035
从
Figure BDA0003724195220000036
中选择样本
Figure BDA0003724195220000037
将以
Figure BDA0003724195220000038
为中心的高斯分布和协方差矩阵∑j,s作为一个假设的μc的建议分布，得到μc作为初始样本分布开始的马尔可夫链；若
Figure BDA0003724195220000039
其中，r为从0到1的均匀分布中的样本，则令
Figure BDA00037241952200000310
更新第s个待辨识参数的先验分布和初始值，否则，重复此步骤。from
Figure BDA0003724195220000034
The index l is randomly generated in , according to the probability
Figure BDA0003724195220000035
from
Figure BDA0003724195220000036
select sample
Figure BDA0003724195220000037
will be
Figure BDA0003724195220000038
centered Gaussian distribution and covariance matrix ∑ j,s as a proposed distribution for a hypothetical μ c , obtaining μ c as a Markov chain starting with the initial sample distribution; if
Figure BDA0003724195220000039
where r is a sample in a uniform distribution from 0 to 1, then let
Figure BDA00037241952200000310
Update the prior distribution and initial value of the s-th parameter to be identified, otherwise, repeat this step.

进一步地，每轮优化开始时，确定自适应缩放因子ξj的值，其计算公式为：Further, at the beginning of each round of optimization, the value of the adaptive scaling factor ξ j is determined, and its calculation formula is:

Figure BDA00037241952200000311
Figure BDA00037241952200000311

Figure BDA00037241952200000312
Figure BDA00037241952200000312

其中，W表示待辨识参数的总数，pr表示样本接受率，tr表示目标接受率，Na表示需要调整的链条数，即第j轮优化中

Figure BDA00037241952200000313
的马尔可夫链的数量。Among them, W represents the total number of parameters to be identified, pr represents the sample acceptance rate, t r represents the target acceptance rate, and Na represents the number of chains to be adjusted, that is, in the jth round of optimization
Figure BDA00037241952200000313
the number of Markov chains.

进一步地，预设置的终止条件为第j轮优化的渐变系数qj等于1，认为每个待辨识参数均找到最优的估算值，停止迭代。Further, the preset termination condition is that the gradient coefficient q j of the jth round of optimization is equal to 1, and it is considered that each parameter to be identified has found the optimal estimated value, and the iteration is stopped.

进一步地，实验的物理参数包括质量和初始刚度，BWBN模型如下：Further, the physical parameters of the experiment include mass and initial stiffness, and the BWBN model is as follows:

Figure BDA00037241952200000314
Figure BDA00037241952200000314

Figure BDA00037241952200000315
Figure BDA00037241952200000315

v(t)＝1.0+δvε(t)v(t)=1.0+δ v ε(t)

η(t)＝1.0+δηε(t)η(t)=1.0+δ η ε(t)

其中，m表示质量，x表示位移，c表示阻尼系数，ke表示刚度，z表示阻尼位移，F(t)表示力，h(z)体现捏缩效应，v(t)表示强度退化，η(t)表示刚度退化，ε(t)表示滞回耗散能量，n、a、β、γ、δv、δη为待辨识参数，

Figure BDA0003724195220000041
和
Figure BDA0003724195220000042
表示x的一阶导数和二阶导数，
Figure BDA0003724195220000043
表示z的一阶导数。where m represents mass, x represents displacement, c represents damping coefficient, ke represents stiffness, z represents damping displacement, F(t) represents force, h(z) represents pinch effect, v(t) represents strength degradation, η (t) represents stiffness degradation, ε(t) represents hysteretic dissipation energy, n, a, β, γ, δ v , δ η are parameters to be identified,
Figure BDA0003724195220000041
and
Figure BDA0003724195220000042
represents the first and second derivatives of x,
Figure BDA0003724195220000043
represents the first derivative of z.

进一步地，使用四阶龙格库塔求解BWBN恢复力测量模型中的微分方程，如下：Further, use the fourth-order Runge-Kutta to solve the differential equation in the BWBN restoring force measurement model, as follows:

Figure BDA0003724195220000044
Figure BDA0003724195220000044

K1＝f(xn,yn)K 1 =f(x n ,y n )

Figure BDA0003724195220000045
Figure BDA0003724195220000045

Figure BDA0003724195220000046
Figure BDA0003724195220000046

K4＝f(xn+h,yn+hK3)K 4 =f(x n +h,y n +hK 3 )

其中，f()表示BWBN模型的微分方程，xn为测量数据中的位移，yn＝F(t)为测量数据中的力，h为龙格库塔时间步长，K1，K2，K3，K4为龙格库塔中的传递参数。where f() represents the differential equation of the BWBN model, x n is the displacement in the measurement data, y n =F(t) is the force in the measurement data, h is the Runge-Kutta time step, K 1 , K 2 , K 3 , K 4 are the transfer parameters in Runge-Kutta.

进一步地，似然函数的建立原理为：Further, the establishment principle of the likelihood function is:

Figure BDA0003724195220000047
Figure BDA0003724195220000047

Figure BDA0003724195220000048
Figure BDA0003724195220000048

ln P3＝ln P1+ln P2 ln P 3 =ln P 1 +ln P 2

其中，σ1表示位移预测误差的未知方差，σ2表示总耗散能量预测误差的未知方差，

Figure BDA0003724195220000049
表示x的平均值，N()表示高斯分布。where σ 1 represents the unknown variance of the displacement prediction error, σ 2 represents the unknown variance of the total dissipated energy prediction error,
Figure BDA0003724195220000049
represents the mean of x, and N() represents a Gaussian distribution.

进一步地，计算各个待辨识参数的似然函数时，需要对预测误差耗散能量进行归一化处理，计算公式为：Further, when calculating the likelihood function of each parameter to be identified, it is necessary to normalize the dissipated energy of the prediction error, and the calculation formula is:

Figure BDA00037241952200000410
Figure BDA00037241952200000410

其中，L＝ln P3，θg为待辨识参数组成的相邻，下标g为测量数据中力和位移的索引，Nj,s为预设置的第j轮优化中第s个待辨识参数的样本抽取值。Among them, L=ln P 3 , θ g is the neighbor composed of the parameters to be identified, the subscript g is the index of the force and displacement in the measured data, N j,s is the preset jth round of optimization to be identified in the sth The sample draw value for the parameter.

进一步地，计算似然函数之前，需要从测量数据的滞回曲线中确定总耗散能量。Further, before calculating the likelihood function, the total dissipated energy needs to be determined from the hysteresis curve of the measured data.

根据本发明的第二方面，提供了一种BWBN模型的参数估计系统，包括：According to a second aspect of the present invention, a parameter estimation system of a BWBN model is provided, comprising:

数据获取模块，用于获取实验的测量数据和物理参数，所述测量数据包括力和位移；a data acquisition module for acquiring experimental measurement data and physical parameters, the measurement data including force and displacement;

BWBN模型模块，用于建立BWBN模型，将测量数据以及物理参数代入BWBN模型，求解BWBN模型，确定BWBN模型中的W个待辨识参数；The BWBN model module is used to establish the BWBN model, substitute the measurement data and physical parameters into the BWBN model, solve the BWBN model, and determine the W parameters to be identified in the BWBN model;

初始化模块，用于对每个待辨识参数进行初始化，确定其最大值、最小值和初始值，假定各个待辨识参数满足高斯分布，分别计算各个待辨识参数的似然函数，得到各个待辨识参数的初始先验分布，将初始先验分布作为待辨识参数的先验分布，将初始值作为待辨识参数的估算值；The initialization module is used to initialize each parameter to be identified, determine its maximum value, minimum value and initial value, and assume that each parameter to be identified satisfies the Gaussian distribution, calculate the likelihood function of each parameter to be identified, and obtain each parameter to be identified. The initial prior distribution of , takes the initial prior distribution as the prior distribution of the parameter to be identified, and takes the initial value as the estimated value of the parameter to be identified;

迭代优化模块，用于执行至少一轮优化操作，并在满足预设置的终止条件，停止迭代，每轮优化操作中分别自每个待辨识参数的先验分布中抽取样本并基于样本更新各个待辨识参数的估算值；The iterative optimization module is used to perform at least one round of optimization operations, and stop the iteration when the preset termination conditions are met. In each round of optimization operations, samples are drawn from the prior distribution of each parameter to be identified and updated based on the samples. Estimated values of identification parameters;

输出模块，用于输出每个待辨识参数的估算值。The output module is used to output the estimated value of each parameter to be identified.

与现有技术相比，本发明具有以下有益效果：Compared with the prior art, the present invention has the following beneficial effects:

本申请引入了可变的自适应调整参数，在每个需调整的马尔可夫链后均进行样本加权系数的调整，从而调整目标分布，提高了算法的性能，减少估计值中的偏差，让样本平均值和方差值估计偏差小，加快搜寻速度，利用过渡马尔科夫链蒙特卡洛方法本身的优势，可有效解决多峰或极平坦单峰问题的采样，同时应用BWBN模型所具有捏缩挤压功能，让系统具有更高的通用性。The application introduces variable adaptive adjustment parameters, and adjusts the sample weighting coefficient after each Markov chain to be adjusted, so as to adjust the target distribution, improve the performance of the algorithm, reduce the deviation in the estimated value, and make the The estimated deviation of the sample mean and variance value is small, the search speed is accelerated, and the advantages of the transitional Markov chain Monte Carlo method can be used to effectively solve the sampling of multi-peak or extremely flat single-peak problems. The shrinking and squeezing function makes the system more versatile.

附图说明Description of drawings

图1为本发明的流程图；Fig. 1 is the flow chart of the present invention;

图2为实施例中归一化耗能图；2 is a normalized energy consumption diagram in an embodiment;

图3为实施例中各个待辨识参数的概率分布图；Fig. 3 is the probability distribution diagram of each parameter to be identified in the embodiment;

图4为实施例中的滞回曲线图。FIG. 4 is a hysteresis curve diagram in the embodiment.

具体实施方式Detailed ways

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. This embodiment is implemented on the premise of the technical solution of the present invention, and provides a detailed implementation manner and a specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

实施例1：Example 1:

针对现有技术中BWBN模型的参数估计方法的缺陷，申请人指出，iTMCMC算法自适应调整建议分布，并在每个MCMC步骤后调整样本权重，所以在平均值和标准偏差偏差以及有效样本数采样方面优于原始TMCMC方法，可更好地应用于多维参数估计问题，故对iTMCMC算法进行改进，使用iTMCMC算法进行BWBN模型的参数估计。In view of the defects of the parameter estimation method of the BWBN model in the prior art, the applicant pointed out that the iTMCMC algorithm adaptively adjusts the proposal distribution, and adjusts the sample weight after each MCMC step, so the mean and standard deviation deviation and the number of valid samples are sampled It is better than the original TMCMC method and can be better applied to the multi-dimensional parameter estimation problem. Therefore, the iTMCMC algorithm is improved, and the iTMCMC algorithm is used to estimate the parameters of the BWBN model.

根据本发明的第一方面，提供了一种BWBN模型的参数估计方法，如图1所示，包括以下步骤：According to the first aspect of the present invention, a method for estimating parameters of a BWBN model is provided, as shown in FIG. 1 , including the following steps:

S1、获取实验的测量数据，测量数据包括力F(t)和位移x，确定实验的物理参数，建立BWBN模型，将测量数据以及物理参数代入BWBN模型，求解BWBN模型，确定BWBN模型中的W个待辨识参数；S1. Obtain the measurement data of the experiment, including force F(t) and displacement x, determine the physical parameters of the experiment, establish a BWBN model, substitute the measurement data and physical parameters into the BWBN model, solve the BWBN model, and determine the W in the BWBN model parameters to be identified;

其中，本实施例中测量数据即滞回力实验中采集的数据，也可以使用数值模拟数据，需要注意的是，测量数据或数值模拟数据表示的滞回图像需符合或近似BWBN模型，这样才能使用BWBN模型对测量数据进行描述。Among them, the measured data in this embodiment is the data collected in the hysteretic force experiment, and numerical simulation data can also be used. It should be noted that the hysteretic image represented by the measured data or numerical simulation data must conform to or approximate the BWBN model, so that it can be used The BWBN model describes the measurement data.

BWBN恢复力测量模型是一种较成熟的数学模型，对其原理和建立过程不再赘述，相关从业人员可以理解。建立的BWBN恢复力测量模型如下：The BWBN resilience measurement model is a relatively mature mathematical model, its principle and establishment process will not be repeated, and relevant practitioners can understand it. The established BWBN resilience measurement model is as follows:

Figure BDA0003724195220000061
Figure BDA0003724195220000061

Figure BDA0003724195220000062
Figure BDA0003724195220000062

v(t)＝1.0+δvε(t)v(t)=1.0+δ v ε(t)

η(t)＝1.0+δηε(t)η(t)=1.0+δ η ε(t)

其中，m表示质量，x表示位移，c表示阻尼系数，ke表示刚度，z表示阻尼位移，F(t)表示力，h(z)体现捏缩效应，v(t)表示强度退化，η(t)表示刚度退化，ε(t)表示滞回耗散能量，n、a、β、γ、δv、δη为待辨识参数，

Figure BDA0003724195220000063
和
Figure BDA0003724195220000064
表示x的一阶导数和二阶导数，
Figure BDA0003724195220000065
表示z的一阶导数。where m represents mass, x represents displacement, c represents damping coefficient, ke represents stiffness, z represents damping displacement, F(t) represents force, h(z) represents pinch effect, v(t) represents strength degradation, η (t) represents stiffness degradation, ε(t) represents hysteretic dissipation energy, n, a, β, γ, δ v , δ η are parameters to be identified,
Figure BDA0003724195220000063
and
Figure BDA0003724195220000064
represents the first and second derivatives of x,
Figure BDA0003724195220000065
represents the first derivative of z.

Figure BDA0003724195220000066
即BWBN模型中的微分方程，实验中的物理参数包括质量和初始刚度等，不再赘述。
Figure BDA0003724195220000066
That is, the differential equation in the BWBN model, the physical parameters in the experiment include mass and initial stiffness, etc., which will not be repeated.

将已知值和测量数据代入BWBN模型后，可由四阶龙格库塔求解，如下：After substituting known values and measurement data into the BWBN model, the fourth-order Runge-Kutta can be solved as follows:

Figure BDA0003724195220000071
Figure BDA0003724195220000071

K1＝f(xn,yn)K 1 =f(x n ,y n )

Figure BDA0003724195220000072
Figure BDA0003724195220000072

Figure BDA0003724195220000073
Figure BDA0003724195220000073

K4＝f(xn+h,yn+hK3)K 4 =f(x n +h,y n +hK 3 )

其中，f()表示BWBN模型的微分方程，xn为测量数据中的位移，yn＝F(t)为测量数据中的力，h为龙格库塔时间步长，K1，K2，K3，K4为龙格库塔中的传递参数。where f() represents the differential equation of the BWBN model, x n is the displacement in the measurement data, y n =F(t) is the force in the measurement data, h is the Runge-Kutta time step, K 1 , K 2 , K 3 , K 4 are the transfer parameters in Runge-Kutta.

S2、分别对每个待辨识参数进行初始化，确定其最大值、最小值和初始值，假定各个待辨识参数满足高斯分布，分别计算各个待辨识参数的似然函数，得到各个待辨识参数的初始先验分布，将初始先验分布作为待辨识参数的先验分布，将初始值作为待辨识参数的估算值；S2, respectively initialize each parameter to be identified, determine its maximum value, minimum value and initial value, assume that each parameter to be identified satisfies the Gaussian distribution, calculate the likelihood function of each parameter to be identified separately, and obtain the initial value of each parameter to be identified Prior distribution, the initial prior distribution is taken as the prior distribution of the parameter to be identified, and the initial value is taken as the estimated value of the parameter to be identified;

四阶龙格库塔求解微分方程后，相当于得到每个待辨识参数的表示，具体的，在BWBN模型中有17个待定参数，并引入2个参数描述输入噪音，因此总计19个待辨识参数。其他应用场景下，根据需要，可以设置不同数量的待辨识参数，如7个、21个等等。按照经验，设置其取值范围和初始值，如下表所示：After the fourth-order Runge-Kutta solves the differential equation, it is equivalent to obtaining the representation of each parameter to be identified. Specifically, there are 17 parameters to be determined in the BWBN model, and 2 parameters are introduced to describe the input noise, so a total of 19 parameters to be identified parameter. In other application scenarios, different numbers of parameters to be identified, such as 7, 21, etc., can be set as required. According to experience, set its value range and initial value, as shown in the following table:

表1待辨识参数的初始化Table 1 Initialization of parameters to be identified

Figure BDA0003724195220000074
Figure BDA0003724195220000074

Figure BDA0003724195220000081
Figure BDA0003724195220000081

计算出系统固有频率，并假设未知的待辨识参数满足高斯分布，计算出似然函数，得到其初始先验分布；似然函数具体建立过程如下，建立以位移均值为均值，总耗能预测误差为方差的高斯分布N()，建立如下的似然函数：Calculate the natural frequency of the system, and assume that the unknown parameter to be identified satisfies the Gaussian distribution, calculate the likelihood function, and obtain its initial prior distribution; the specific establishment process of the likelihood function is as follows, the establishment takes the displacement mean as the mean value, and the total energy consumption prediction error For the Gaussian distribution of variance N(), the following likelihood function is established:

Figure BDA0003724195220000082
Figure BDA0003724195220000082

Figure BDA0003724195220000083
Figure BDA0003724195220000083

ln P3＝ln P1+ln P2 ln P 3 =ln P 1 +ln P 2

其中，σ1表示位移预测误差的未知方差，σ2表示总耗散能预测误差的未知方差，

Figure BDA0003724195220000084
表示x的平均值，N()表示高斯分布；where σ 1 represents the unknown variance of the displacement prediction error, σ 2 represents the unknown variance of the total dissipated energy prediction error,
Figure BDA0003724195220000084
Represents the mean value of x, and N() represents a Gaussian distribution;

再偏导求解参数值，从而通过计算各个待辨识参数的似然函数，得到各个待辨识参数的初始先验分布。Then the partial derivative is used to solve the parameter value, so that the initial prior distribution of each parameter to be identified can be obtained by calculating the likelihood function of each parameter to be identified.

在计算各个待辨识参数的似然函数时，需要对预测误差耗散能量进行归一化处理，图2给出了归一化耗能图的示例，计算公式为：When calculating the likelihood function of each parameter to be identified, the prediction error dissipation energy needs to be normalized. Figure 2 shows an example of the normalized energy dissipation map. The calculation formula is:

Figure BDA0003724195220000085
Figure BDA0003724195220000085

其中，L＝ln P3，θg为待辨识参数组成的相邻，下标g为测量数据中力和位移的索引，Nj,s为预设置的第j轮优化中第s个待辨识参数的样本抽取值。Among them, L=ln P 3 , θ g is the neighbor composed of the parameters to be identified, the subscript g is the index of the force and displacement in the measured data, N j,s is the preset jth round of optimization to be identified in the sth The sample draw value for the parameter.

在计算似然函数之前，需要从测量数据的滞回曲线中确定总耗散能量，其计算方法为本领域公知常识，不再赘述。Before calculating the likelihood function, it is necessary to determine the total dissipated energy from the hysteresis curve of the measurement data, and the calculation method thereof is common knowledge in the art, and will not be repeated here.

S3、分别自每个待辨识参数的先验分布中抽取样本并基于样本更新各个待辨识参数的估算值，完成一轮优化，重复此步骤，直至满足预设置的终止条件，停止迭代；S3, extract samples from the prior distribution of each parameter to be identified and update the estimated value of each parameter to be identified based on the sample, complete a round of optimization, and repeat this step until the preset termination condition is met, and the iteration is stopped;

以第j轮优化中自第s个待辨识参数的先验分布中抽取样本并基于样本更新第s个待辨识参数的估算值为例，对迭代优化过程进行说明，其中，j＝1、2、…，s＝1、2、…、W，具体为：Taking samples from the prior distribution of the s-th parameter to be identified in the j-th round of optimization and updating the estimated value of the s-th parameter to be identified based on the samples as an example, the iterative optimization process is described, where j=1, 2 , ..., s = 1, 2, ..., W, specifically:

(1)从第s个待辨识参数的先验分布中抽取Nj,s个样本：

Figure BDA0003724195220000091
Nj,s为预设置的第j轮优化中第s个待辨识参数的样本抽取值，每个待辨识参数在每轮优化中的样本抽取值可以各不相同，本实施例中，为了便于计算，统一规定Nj,s＝1000；(1) Extract N j,s samples from the prior distribution of the s-th parameter to be identified:
Figure BDA0003724195220000091
N j,s is the preset sample extraction value of the sth parameter to be identified in the jth round of optimization, and the sample extraction value of each parameter to be identified in each round of optimization may be different. In this embodiment, for convenience Calculate, uniformly stipulate N j,s =1000;

(2)计算第j轮优化的渐变系数，预设置的终止条件为第j轮优化的渐变系数qj等于1，认为每个待辨识参数均找到最优的估算值，因此，在此步骤中，如果计算出的渐变系数为0，则停止迭代，完成优化，执行步骤S4。此外，可以设置最大迭代次数作为终止条件。渐变系数的计算公式如下：(2) Calculate the gradient coefficient of the jth round of optimization, the preset termination condition is that the gradient coefficient q j of the jth round of optimization is equal to 1, and it is considered that each parameter to be identified finds the optimal estimated value. Therefore, in this step , if the calculated gradient coefficient is 0, stop the iteration, complete the optimization, and execute step S4. Also, a maximum number of iterations can be set as a termination condition. The formula for calculating the gradient coefficient is as follows:

q0＝0，qj＝min(|CV(τj,s)-1|，s＝1、2、...、W)q 0 =0, q j =min(|CV(τ j,s )-1|, s=1, 2, . . . , W)

其中，CV(τj,s)表示第j轮优化中第s个待辨识参数的样本变异系数，min(|CV(τj,s)-1|，s＝1、2、...、W)即表示计算第j轮中19个待辨识参数的样本变异系数，并从中选取与1的差值最小的样本变异系数作为第j轮优化的渐变系数，样本变异系数是本领域常用技术手段，其计算公式和原理不再赘述，相关从业人员可以理解；Among them, CV(τ j,s ) represents the sample variation coefficient of the sth parameter to be identified in the jth round of optimization, min(|CV(τ j,s )-1|, s=1, 2,..., W) means to calculate the sample variation coefficient of the 19 parameters to be identified in the jth round, and select the sample variation coefficient with the smallest difference from 1 as the gradient coefficient of the jth round optimization, and the sample variation coefficient is a common technical means in this field. , its calculation formula and principle will not be repeated, and relevant practitioners can understand;

(3)分别计算每个样本的加权系数，计算公式如下：(3) Calculate the weighting coefficient of each sample separately, and the calculation formula is as follows:

Figure BDA0003724195220000092
Figure BDA0003724195220000092

Figure BDA0003724195220000093
表示第j轮优化中第s个待辨识参数的第k个样本的加权系数，D表示实验的测量数据；
Figure BDA0003724195220000093
represents the weighting coefficient of the k-th sample of the s-th parameter to be identified in the j-th round of optimization, and D represents the experimental measurement data;

(4)计算Nj,s个样本的加权系数的平均值，计算公式如下：(4) Calculate the average value of the weighting coefficients of N j, s samples, and the calculation formula is as follows:

Figure BDA0003724195220000094
Figure BDA0003724195220000094

其中，Sej,s表示第j轮优化中第s个待辨识参数的样本加权系数平均值；Among them, Se j,s represents the average value of the sample weighting coefficient of the s-th parameter to be identified in the j-th round of optimization;

(5)计算协方差矩阵，计算公式如下：(5) Calculate the covariance matrix, the calculation formula is as follows:

Figure BDA0003724195220000095
Figure BDA0003724195220000095

其中，∑j,s表示第j轮优化中第s个待辨识参数的Nj,s个样本的协方差矩阵，

Figure BDA0003724195220000096
表示第j轮优化中第s个待辨识参数的Nj,s个样本的均值，T表示矩阵转置，ξj表示预设置的第j轮优化中的自适应缩放因子；Among them, ∑ j,s represents the covariance matrix of N j,s samples of the sth parameter to be identified in the jth round of optimization,
Figure BDA0003724195220000096
represents the mean value of N j, s samples of the s-th parameter to be identified in the j-th round of optimization, T represents the matrix transposition, and ξ j represents the preset adaptive scaling factor in the j-th round of optimization;

(6)从

Figure BDA0003724195220000101
中随机生成索引l，按照概率
Figure BDA0003724195220000102
从
Figure BDA0003724195220000103
中选择样本
Figure BDA0003724195220000104
将以
Figure BDA0003724195220000105
为中心的高斯分布和协方差矩阵∑j,s作为一个假设的μc的建议分布，得到μc作为初始样本分布开始的马尔可夫链；若
Figure BDA0003724195220000106
其中，r为从0到1的均匀分布中的样本，则令
Figure BDA0003724195220000107
更新第s个待辨识参数的先验分布和初始值，否则，重复此步骤。(6) From
Figure BDA0003724195220000101
The index l is randomly generated in , according to the probability
Figure BDA0003724195220000102
from
Figure BDA0003724195220000103
select sample
Figure BDA0003724195220000104
will be
Figure BDA0003724195220000105
centered Gaussian distribution and covariance matrix ∑ j,s as a proposed distribution for a hypothetical μ c , obtaining μ c as a Markov chain starting with the initial sample distribution; if
Figure BDA0003724195220000106
where r is a sample in a uniform distribution from 0 to 1, then let
Figure BDA0003724195220000107
Update the prior distribution and initial value of the s-th parameter to be identified, otherwise, repeat this step.

生成从0到1的均匀分布，主要是用于评判此次采样的概率分布是否小于分布，如果均匀分布的概率大于采样的样本概率，即

Figure BDA0003724195220000108
说明此次采样的链条是错误的，可以放弃，重新生成索引l，如果
Figure BDA0003724195220000109
说明可以接受此次采样，继而使用μc更新第s个待辨识参数的先验分布和初始值。Generating a uniform distribution from 0 to 1 is mainly used to judge whether the probability distribution of this sampling is smaller than the distribution. If the probability of the uniform distribution is greater than the sample probability of the sampling, that is
Figure BDA0003724195220000108
It means that the sampling chain is wrong, you can give up and regenerate the index l, if
Figure BDA0003724195220000109
It shows that this sampling can be accepted, and then use μ c to update the prior distribution and initial value of the s-th parameter to be identified.

完成上述步骤(1)-(6)后，即完成了第j轮优化中自第s个待辨识参数的先验分布和估算值更新，再令s+1，再次执行步骤(1)-(6)，完成下一个待辨识参数在第j轮的优化。完成全部19个待辨识参数的优化后，由于没有满足预设置的终止条件，令j+1，开始下一轮的迭代优化。After completing the above steps (1)-(6), the prior distribution and estimated value update of the s-th parameter to be identified in the j-th round of optimization is completed, and then let s+1, and execute steps (1)-( 6), to complete the optimization of the next parameter to be identified in the jth round. After completing the optimization of all 19 parameters to be identified, since the preset termination conditions are not met, let j+1, and start the next round of iterative optimization.

在步骤(5)中，自适应缩放因子ξj是在每轮优化开始时确定的，其计算公式为：In step (5), the adaptive scaling factor ξ j is determined at the beginning of each round of optimization, and its calculation formula is:

Figure BDA00037241952200001010
Figure BDA00037241952200001010

Figure BDA00037241952200001011
Figure BDA00037241952200001011

其中，W表示待辨识参数的总数，pr表示样本接受率，tr表示目标接受率，Na表示需要调整的链条数，即第j轮优化中

Figure BDA00037241952200001012
的马尔可夫链的数量。Among them, W represents the total number of parameters to be identified, pr represents the sample acceptance rate, t r represents the target acceptance rate, and Na represents the number of chains to be adjusted, that is, in the jth round of optimization
Figure BDA00037241952200001012
the number of Markov chains.

S4、输出每个待辨识参数的估算值，以及概率密度函数pdf，其中，本实施例中待辨识参数的概率分布如图3所示。S4 , output the estimated value of each parameter to be identified, and the probability density function pdf, wherein the probability distribution of the parameter to be identified in this embodiment is shown in FIG. 3 .

本实施例中，经过多轮优化后，得到19个待辨识参数的估算值和方差值，如下表所示：In this embodiment, after multiple rounds of optimization, the estimated values and variance values of 19 parameters to be identified are obtained, as shown in the following table:

表2待辨识参数的估算值Table 2 Estimated values of parameters to be identified

Figure BDA00037241952200001013
Figure BDA00037241952200001013

Figure BDA0003724195220000111
Figure BDA0003724195220000111

本实施例中，得到BWBN模型中各个参数的估算值后，就可以将其代回BWBN模型，再次将实验的测量数据和物理参数代入BWBN模型，求解出实验的滞回曲线，模拟滞回曲线图如图4所示。In this embodiment, after obtaining the estimated value of each parameter in the BWBN model, it can be substituted back into the BWBN model, and the measured data and physical parameters of the experiment are substituted into the BWBN model again, the hysteresis curve of the experiment is solved, and the hysteresis curve is simulated The diagram is shown in Figure 4.

根据本发明的第二方面，提供了一种BWBN模型的参数估计系统，包括：According to a second aspect of the present invention, a parameter estimation system of a BWBN model is provided, comprising:

数据获取模块，用于获取实验的测量数据和物理参数，测量数据包括力和位移；The data acquisition module is used to acquire the measurement data and physical parameters of the experiment, and the measurement data includes force and displacement;

BWBN模型模块，用于建立BWBN模型，将测量数据以及物理参数代入BWBN模型，求解BWBN模型，确定BWBN模型中的W个待辨识参数；The BWBN model module is used to establish the BWBN model, substitute the measurement data and physical parameters into the BWBN model, solve the BWBN model, and determine the W parameters to be identified in the BWBN model;

初始化模块，用于对每个待辨识参数进行初始化，确定其最大值、最小值和初始值，假定各个待辨识参数满足高斯分布，分别计算各个待辨识参数的似然函数，得到各个待辨识参数的初始先验分布，将初始先验分布作为待辨识参数的先验分布，将初始值作为待辨识参数的估算值；The initialization module is used to initialize each parameter to be identified, determine its maximum value, minimum value and initial value, and assume that each parameter to be identified satisfies the Gaussian distribution, calculate the likelihood function of each parameter to be identified, and obtain each parameter to be identified. The initial prior distribution of , takes the initial prior distribution as the prior distribution of the parameter to be identified, and takes the initial value as the estimated value of the parameter to be identified;

迭代优化模块，用于执行至少一轮优化操作，并在满足预设置的终止条件，停止迭代，每轮优化操作中分别自每个待辨识参数的先验分布中抽取样本并基于样本更新各个待辨识参数的估算值；The iterative optimization module is used to perform at least one round of optimization operations, and stop the iteration when the preset termination conditions are met. In each round of optimization operations, samples are drawn from the prior distribution of each parameter to be identified and updated based on the samples. Estimated values of identification parameters;

输出模块，用于输出每个待辨识参数的估算值。The output module is used to output the estimated value of each parameter to be identified.

iTMCMC采样方法基于贝叶斯模型更新技术，应用iTMCMC采样方法估计BWBN非线性迟滞模型的参数，可以估计除非敏感滞后振幅参数A以外的多维参数向量，返回参数后验分布、估算值和方差。相较于MCMC采样方法，本申请调整每个MCMC采样步骤后的样本权重，以降低模型证据估计的平均偏差，在MCMC步骤中应用老化期，以改善后验近似，自适应选择MCMC算法的建议目标分布区间，以实现接近最优的接受率。因此，本申请在克服了已有算法的缺点同时保留了MCMC采样的优点。The iTMCMC sampling method is based on the Bayesian model update technology. The iTMCMC sampling method is used to estimate the parameters of the BWBN nonlinear hysteresis model, which can estimate multi-dimensional parameter vectors other than the non-sensitive lag amplitude parameter A, and return the parameter posterior distribution, estimated value and variance. Compared with the MCMC sampling method, the present application adjusts the sample weight after each MCMC sampling step to reduce the average deviation of the model evidence estimation, applies an aging period in the MCMC step to improve the posterior approximation, and adaptively selects the MCMC algorithm suggestions Target distribution interval to achieve near-optimal acceptance rates. Therefore, the present application retains the advantages of MCMC sampling while overcoming the shortcomings of the existing algorithms.

应用本申请可带来以下优点：The application of this application can bring the following advantages:

采样的独立样本有效数量多，采样效率高，平均值与方差值估计偏差较小，可有效解决多峰或极平坦单峰问题的采样，适用范围广，可很好模拟有捏缩滞回现象的滞回图像，实现简便，不需要滤波器实现有噪声影响的参数估计，可提供BWBN模型参数的概率分布，对模型选择的评估有重要意义。The effective number of independent samples to be sampled is large, the sampling efficiency is high, and the estimated deviation between the average value and the variance value is small, which can effectively solve the sampling of multi-peak or extremely flat single-peak problems. The hysteresis image of the phenomenon is easy to implement and does not require a filter to achieve parameter estimation with noise effects, and can provide the probability distribution of the parameters of the BWBN model, which is of great significance to the evaluation of model selection.

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。The preferred embodiments of the present invention have been described in detail above. It should be understood that those skilled in the art can make numerous modifications and changes according to the concept of the present invention without creative efforts. Therefore, all technical solutions that can be obtained by those skilled in the art through logical analysis, reasoning or limited experiments on the basis of the prior art according to the concept of the present invention shall fall within the protection scope determined by the claims.

Claims (10)

1. A method for estimating parameters of a BWBN model is characterized by comprising the following steps:
acquiring measurement data of an experiment, wherein the measurement data comprises force and displacement, determining physical parameters of the experiment, establishing a BWBN model, substituting the measurement data and the physical parameters into the BWBN model, solving the BWBN model, and determining W parameters to be identified in the BWBN model;
initializing each parameter to be identified respectively, determining the maximum value, the minimum value and the initial value of each parameter to be identified, supposing that each parameter to be identified meets Gaussian distribution, calculating the likelihood function of each parameter to be identified respectively to obtain initial prior distribution of each parameter to be identified, taking the initial prior distribution as the prior distribution of the parameter to be identified, and taking the initial value as the estimated value of the parameter to be identified;
respectively extracting samples from the prior distribution of each parameter to be identified, updating the estimated value of each parameter to be identified based on the samples, completing a round of optimization, repeating the step until a preset termination condition is met, and stopping iteration;
and outputting the estimated value of each parameter to be identified.
2. The method of claim 1, wherein in the j-th optimization, samples are extracted from the prior distribution of the s-th parameter to be identified, and the estimated value of the s-th parameter to be identified is updated based on the samples, wherein j =1, 2, \8230, s =1, 2, \8230, and W is specifically:
extracting N from the prior distribution of the s-th parameter to be identified j,s One sample:
Figure FDA0003724195210000011
N j,s extracting a value for the sample of the s-th parameter to be identified in the preset j-th optimization;
and (4) calculating the optimized gradient coefficient of the j-th round, wherein the calculation formula is as follows:
Figure FDA0003724195210000012
wherein,
Figure FDA0003724195210000013
representing the sample variation coefficient of the s-th parameter to be identified in the j-th optimization;
and respectively calculating the weighting coefficient of each sample, wherein the calculation formula is as follows:
Figure FDA0003724195210000014
Figure FDA0003724195210000015
representing the weighting coefficient of the kth sample of the s-th parameter to be identified in the j-th optimization, and D represents the measurement data of the experiment;
calculating N j，s The average value of the weighting coefficients of the samples is calculated according to the following formula:
Figure FDA0003724195210000016
wherein, se j,s Representing the sample weighting coefficient average value of the s-th parameter to be identified in the j-th optimization;
and (3) calculating a covariance matrix, wherein the calculation formula is as follows:
Figure FDA0003724195210000021
therein, sigma j,s N for representing s to-be-identified parameter in j round optimization j,s The covariance matrix of the individual samples is,
Figure FDA0003724195210000022
n for representing s to-be-identified parameter in j optimization j,s Mean value of samples, T denotes matrix transpose, ξ j Representing the preset adaptive scaling factor in the jth round of optimization;
from
Figure FDA0003724195210000023
Generating index l at random according to probability
Figure FDA0003724195210000024
From
Figure FDA0003724195210000025
In selecting samples
Figure FDA0003724195210000026
Will be provided with
Figure FDA0003724195210000027
Centered gaussian distribution and covariance matrix sigma j,s As a hypothesis c Of the proposed distribution, to obtain μ c A Markov chain that starts as an initial sample distribution; if it is
Figure FDA0003724195210000028
Where r is a sample in a uniform distribution from 0 to 1, then let
Figure FDA0003724195210000029
And updating the prior distribution and the initial value of the s-th parameter to be identified, otherwise, repeating the step.
3. The method of claim 2, wherein the adaptive scaling factor ξ is determined at the beginning of each round of optimization j The calculation formula is as follows:
Figure FDA00037241952100000210
Figure FDA00037241952100000211
wherein W represents the total number of parameters to be identified, p r Denotes sample acceptance rate, t r Indicating the target acceptance rate, na indicating the number of chains to be adjusted, i.e. in the j-th optimization
Figure FDA00037241952100000212
Number of markov chains.
4. The method of claim 2 for estimating the parameters of a BWBN model, whereinCharacterized in that the preset termination condition is the gradient coefficient q optimized in the j' th round j And if the value is equal to 1, considering that each parameter to be identified finds the optimal estimated value, and stopping iteration.
5. The method of claim 1, wherein the physical parameters of the experiment include mass and initial stiffness, and the BWBN model is as follows:
Figure FDA00037241952100000213
Figure FDA00037241952100000214
v(t)＝1.0+δ v ε(t)
η(t)＝1.0+δ η ε(t)
wherein m represents mass, x represents displacement, c represents damping coefficient, k e Representing stiffness, z damping displacement, F (t) force, h (z) pinch effect, v (t) strength degradation, η (t) stiffness degradation, ε (t) hysteresis dissipation energy, n, a, β, γ, δ v 、δ η As the parameter to be identified, the identification information is obtained,
Figure FDA0003724195210000031
and
Figure FDA0003724195210000032
the first and second derivatives of x are indicated,
Figure FDA0003724195210000033
representing the first derivative of z.
6. The parameter estimation method of the BWBN model according to claim 5, wherein the differential equation in the BWBN restoring force measurement model is solved using a fourth order longgure tower as follows:
Figure FDA0003724195210000034
K 1 ＝f(x n ,y n )
Figure FDA0003724195210000035
Figure FDA0003724195210000036
K 4 ＝f(x n +h,y n +hK 3 )
wherein f () represents a differential equation of the BWBN model, x n For measuring displacements in the data, y n = F (t) force in the measurement data, h is the Longgasta time step, K 1 ，K 2 ，K 3 ，K 4 Are transfer parameters in the dragon lattice tower.
7. The method of claim 5 wherein the likelihood function is built based on the following principles:
Figure FDA0003724195210000037
Figure FDA0003724195210000038
lnP 3 ＝lnP 1 +lnP 2
wherein σ 1 Representing the unknown variance, σ, of the prediction error of the displacement 2 An unknown variance representing the total dissipated energy prediction error,
Figure FDA0003724195210000039
represents the average value of x, and N () represents a gaussian distribution.
8. The method of claim 7, wherein the likelihood function of each parameter to be identified is calculated by normalizing the predicted error dissipation energy, and the calculation formula is:
Figure FDA00037241952100000310
wherein, L = lnP 3 ，θ g For the neighbourhood of the composition of the parameters to be identified, the subscript g is the index of the forces and displacements in the measured data, N j,s And extracting a sample value of the s-th parameter to be identified in the preset j-th optimization.
9. The method of claim 8, wherein the total dissipated energy is determined from a hysteresis curve of the measured data before the likelihood function is calculated.
10. A parameter estimation system of a BWBN model, characterized in that the parameter estimation method based on the BWBN model according to any one of claims 1-9 comprises:
the data acquisition module is used for acquiring measurement data and physical parameters of an experiment, wherein the measurement data comprises force and displacement;
the BWBN model module is used for establishing a BWBN model, substituting the measurement data and the physical parameters into the BWBN model, solving the BWBN model and determining W parameters to be identified in the BWBN model;
the initialization module is used for initializing each parameter to be identified, determining the maximum value, the minimum value and the initial value of the parameter to be identified, respectively calculating the likelihood function of each parameter to be identified on the assumption that each parameter to be identified meets Gaussian distribution to obtain initial prior distribution of each parameter to be identified, taking the initial prior distribution as the prior distribution of the parameter to be identified, and taking the initial value as the estimated value of the parameter to be identified;
the iteration optimization module is used for executing at least one round of optimization operation, stopping iteration when a preset termination condition is met, respectively extracting samples from the prior distribution of each parameter to be identified in each round of optimization operation, and updating the estimation value of each parameter to be identified based on the samples;
and the output module is used for outputting the estimation value of each parameter to be identified.
CN202210770940.3A 2022-06-30 2022-06-30 Method and system for estimating BWBN model parameters Pending CN115183969A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210770940.3A CN115183969A (en) 2022-06-30 2022-06-30 Method and system for estimating BWBN model parameters

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210770940.3A CN115183969A (en) 2022-06-30 2022-06-30 Method and system for estimating BWBN model parameters

Publications (1)

Publication Number Publication Date
CN115183969A true CN115183969A (en) 2022-10-14

Family

ID=83514573

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210770940.3A Pending CN115183969A (en) 2022-06-30 2022-06-30 Method and system for estimating BWBN model parameters

Country Status (1)

Country Link
CN (1) CN115183969A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117216846A (en) * 2023-09-12 2023-12-12 华南理工大学 Reinforced concrete member hysteresis curve prediction method, system, equipment and medium

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106909738A (en) * 2017-02-24 2017-06-30 北京工商大学 A kind of model parameter identification method

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106909738A (en) * 2017-02-24 2017-06-30 北京工商大学 A kind of model parameter identification method

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
GILBERTO A. ORTIZ 等: "Identification of Bouc–Wen type models using the Transitional Markov Chain Monte Carlo method", COMPUTERS AND STRUCTURES, vol. 146, 6 November 2014 (2014-11-06), pages 252 - 269 *
ZHAO LU 等: "Simulation Study of New Buckling Restraint Bracing Based on Bayesian Optimization Algorithm to Identify the Parameters of Modified Bouc–Wen Model", ADVANCES IN CIVIL ENGINEERING, vol. 2023, 3 January 2023 (2023-01-03), pages 1 - 12 *
赵露: "新型摩擦阻尼器的设计与性能研究", 工程科技Ⅱ, 18 May 2023 (2023-05-18), pages 1 - 65 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117216846A (en) * 2023-09-12 2023-12-12 华南理工大学 Reinforced concrete member hysteresis curve prediction method, system, equipment and medium
CN117216846B (en) * 2023-09-12 2024-04-19 华南理工大学 Method, system, equipment and medium for predicting hysteresis curve of reinforced concrete components

Similar Documents

Publication Publication Date Title
CN111913803B (en) Service load fine granularity prediction method based on AKX hybrid model
CN113391211B (en) A method for predicting the remaining life of lithium batteries under the condition of small samples
CN108090621B (en) Short-term wind speed prediction method and system based on staged overall optimization
KR20170056687A (en) Actually-measured marine environment data assimilation method based on sequence recursive filtering three-dimensional variation
CN108445759B (en) A Random Fault Detection Method for Networked Systems Under Sensor Saturation Constraints
CN113240170A (en) Air quality prediction method based on seasonal cyclic neural network
CN116500454A (en) A method, system, device, and medium for estimating the state of health of a lithium-ion battery based on a multi-feature input time-series model
CN108304685A (en) A kind of non-linear degradation equipment method for predicting residual useful life and system
CN110677297A (en) Combined network flow prediction method based on autoregressive moving average model and extreme learning machine
CN113852432A (en) RCS-GRU model-based spectrum prediction sensing method
CN108734287A (en) Compression method and device, terminal, the storage medium of deep neural network model
CN112086100A (en) Quantization error entropy based urban noise identification method of multilayer random neural network
CN111460692B (en) Method and system for predicting remaining life of equipment considering mutual influence of degradation rates
CN102063524A (en) Performance reliability simulation method based on improved self-adaption selective sampling
CN108734268A (en) Compression method and device, terminal, the storage medium of deep neural network model
CN110471768B (en) FastPCA-ARIMA-based load prediction method
CN118114118A (en) A typical weapon equipment fault diagnosis method based on CNDT
CN115183969A (en) Method and system for estimating BWBN model parameters
CN116227324B (en) Fractional order memristor neural network estimation method under variance limitation
CN118189870A (en) A bridge displacement prediction and early warning method and system
CN113837443A (en) Transformer substation line load prediction method based on depth BilSTM
CN109978138A (en) The structural reliability methods of sampling based on deeply study
CN111325308A (en) A Nonlinear System Identification Method
CN119202604A (en) A method and system for rapid prediction of ship cabin noise based on convolutional neural network
CN118535851A (en) A sparse regular graph pooling network sensor optimization configuration method for diesel engine vibration monitoring system

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information

Address after: 200437 No. 99, Handan Road, Shanghai, Hongkou District

Applicant after: Shanghai Material Research Institute Co.,Ltd.

Address before: 200437 No. 99, Handan Road, Shanghai, Hongkou District

Applicant before: SHANGHAI Research Institute OF MATERIALS

CB02 Change of applicant information