A robust identification method for ARX model based on skew measurement noise

Abstract

The present invention provides a robust identification method of an ARX model based on skew measurement noise, comprising: step 1: collecting identification data according to the continuous time dynamics equation of a mass spring damping system, and constructing an ARX model of a mass spring damping system outputted by a simulated position corresponding to the identification data; step 2: establishing an ARX robust probability model based on a generalized hyperbolic tilted Student's t distribution according to the ARX model; step 3: identifying the ARX robust probability model based on a generalized hyperbolic tilted Student's t distribution established in step 2 based on a generalized expectation maximization algorithm, and when the iteration converges, the final iteration result is used as the optimal estimated parameter. The method of the present invention can effectively improve the accuracy of parameter estimation of a mass spring damping system, and is of great significance to ensuring the accuracy of damping vibration reduction in actual industry.

Description

Technical Field

[0001] The invention belongs to the technical field of system modeling and parameter identification, and in particular relates to an ARX model robust identification method based on skew measurement noise. Background Art

[0002] Auto-Regressive with Extra Inputs (ARX) is an important time series analysis model that can describe many problems in practical applications and has a wide range of applications in industrial systems, vehicle driving, power systems, medical treatment, signal processing, and other fields. The obtained observation data is used to accurately estimate the model parameters of the system, thereby establishing an accurate ARX system to fit the actual system. Accurate modeling of these systems is a prerequisite for the design and implementation of control and optimization methods. Therefore, the relevant theories and methods based on the ARX system have received a lot of attention.

[0003] It is worth noting that in practical applications, due to problems such as sensor failure and external environmental interference, the collected observation data often contain skew noise or a certain proportion of outliers; these problems greatly affect the quality of the collected data set and put forward higher requirements for the identification algorithm of the actual process. Summary of the invention

[0004] Purpose of the invention: The technical problem to be solved by the present invention is to provide a robust identification method of the ARX model based on skew measurement noise in view of the shortcomings of the prior art, and to improve the accuracy of algorithm parameter identification in view of the problem that the output data of the existing industrial process is affected by skew noise and outliers, so as to better fit the actual system. The method of the present invention is first used for the identification problem of the mass spring damping system, which can effectively solve the problem of reduced modeling accuracy caused by skew noise and outliers.

[0005] The present invention comprises the following steps:

[0006] Step 1: Collect identification data according to the actual process of the mass-spring-damper system, and construct the ARX model of the mass-spring-damper system corresponding to the identification data:

[0007] Step 2: Based on the ARX model constructed in step 1, an ARX robust probability model based on generalized hyperbolic skewed Student's t distribution is established;

[0008] Step 3: Based on the generalized expectation maximization algorithm, the ARX robust probability model established in step 2 is identified. When the iteration converges, the final iterative result is used as the optimal estimated parameter.

[0009] Step 1 includes:

[0010] The dynamic equations of the actual variables of the mass-spring-damper system are written as:

[0011]

[0012] Where x(t) represents the position of the mass block, k s (t) is the stiffness of the spring, c d represents the damping coefficient of the fixed damper, F(t) is the external force applied to mass m;

[0013] In the actual system, the external force F(t) represents the input u(t) of the mass-spring-damper system, the mass block position x(t) represents the output y(t) of the mass-spring-damper system, and N pairs of input and output data of the mass-spring-damper system are collected at a time interval of T = 50 ms, where N represents the length of the data moment;

[0014] The following ARX model is constructed based on the input u(t) and output y(t) of the mass-spring-damper system:

[0015] A(z -1 )y k =B(z -1 ) k +e k

[0016]

[0017]

[0018] where u k represents the input of the mass-spring-damper system, y k represents the output of the mass-spring-damper system, e k represents the measurement noise of the mass spring damping system output, the subscript k represents the identification data moment, and the value ranges from 1 to N; A(z -1 ) and B(z -1 ) is the unit shift operator z -1 The polynomial of z -1 y k =y k-1 , the superscript of z represents the magnitude of the backward shift moment, for example and Respectively represent A(z -1 ) and B(z -1 ) is the parameter corresponding to the polynomial, the subscript represents the corresponding backshift operator z term, which is the unknown parameter to be identified; n a and n b They are A(z -1 ) and B(z -1 ) is a known order, which indicates the known number of polynomial parameters;

[0019] Rewrite the ARX model into a linear regression form, define the parameter vector θ and the input and output information vector for:

[0020]

[0021]

[0022] Where T represents the matrix transposition operator, and the input and output information vectors The subscripts k-1 to kn of the input and output a , k-1 to kn b Indicates the moment;

[0023] The final linear regression system model is:

[0024]

[0025] In step 2, an autoregressive ARX robust probability model with external input based on generalized hyperbolic skewed Student's t distribution is established, specifically including:

[0026] Step 2-1: Measurement noise e at the output at time k k Obey the generalized hyperbolic skew Student's t distribution:

[0027] p(e k |β,Λ,υ)=GHSkewt(e k |0,β,Λ,υ)

[0028] Where p(e k |β,Λ,υ) represents the measurement noise e under β,Λ,υ parameters k The posterior probability of , GHSkewt represents the generalized hyperbolic skew Student's t distribution, β, Λ, υ represent the skew parameter, scale parameter and degree of freedom parameter of GHSkewt distribution respectively;

[0029] Step 2-2: Introduce weight variable w k Decompose the noise distribution in step 2-1:

[0030]

[0031] Where p(e k |w k ,β,Λ,υ) indicates that k ,β,Λ,υ parameter measurement noise e k The posterior probability, p(w k |υ) represents the weight variable w k The posterior probability under the scale parameter υ follows the following distribution:

[0032] p(w k |υ)=IG(w k |υ / 2,υ / 2)

[0033] Among them, IG(w k |υ / 2,υ / 2) represents the inverse gamma distribution, υ represents the degree of freedom parameter;

[0034] When introducing the weight variable w k Down:

[0035] p(e k |w k ,β,Λ,υ)=N(e k |βw k ,Λw k )

[0036] Where N(e k |βw k ,Λw k ) represents Gaussian distribution, βw k represents the mean of the Gaussian distribution, Λw k represents the variance of the Gaussian distribution;

[0037] Step 2-3: Combine the system model in the form of linear regression Get the system output y k The generalized hyperbolic skew Student's t distribution is:

[0038]

[0039] in is a sub-Gaussian distribution, represents the mean of the sub-Gaussian, Λw k represents the variance of each sub-Gaussian;

[0040] The parameter set for constructing the ARX model of the mass-spring-damper system is Θ = {θ, β, Λ, υ}.

[0041] Step 3 includes:

[0042] Step 3-1, initialization: Initialize the parameter Θ to be estimated 1 , the number of iterations s is set to 1;

[0043] Step 3-2, execute step E-step: calculate and derive the Q-function Q(Θ|Θ s ):

[0044]

[0045] Where logp(C mis ,C obs |Θ) is the complete data set C={Cmis ,C obs}, Θ s Estimated parameter set for the model of the current linear system, C obs ={u1: N ,y1: N} is the observed data set, u 1:N Represents the input data of the system at each moment, y 1:N Represents the output data of the system at each moment, C mis ={w1: N} is a missing data set, w 1:N is the weight sequence, Θ = {θ, β, Λ, υ} is the model parameter set; represents logp(C mis ,C obs |Θ) conditional expectation;

[0046] Step 3-3, execute step M-step: by maximizing the Q-function, obtain a new parameter estimate Θ s+1 , so that the parameter Θ s+1 Satisfy the constraints:

[0047] Q(Θ s+1 |Θ s )≥Q(Θ s |Θ s )

[0048] Let: Θ s =Θ s+1 ;

[0049] Step 3-4, let s=s+1, return to step 3-2, and repeat step E-step and step M-step until the parameters to be estimated converge.

[0050] In step 3-2, step E-step includes:

[0051] Step 3-2-1, calculate the complete data set C = {C mis ,C obs The log-likelihood function logp(C mis ,C obs |Θ), using the probability chain rule, the log-likelihood function is written as:

[0052]

[0053]

[0054] Where C1=p(u 1:N |Θ) is a constant term;

[0055] The log-likelihood function is applied to the weight sequence w 1:N Taking the conditional expectation, we get:

[0056]

[0057] Where C1=p(u 1:N |Θ) is a constant term, Θ s ={θ s ,β s ,Λ s ,υ s} represents the result of the sth parameter iteration, θ s ,β s ,Λ s ,υ s They represent the sth iteration results of model parameters, skew parameters, scale parameters and degrees of freedom parameters, respectively, for the weight sequence w 1:N Find the integral term of the conditional expectation:

[0058]

[0059]

[0060]

[0061] Step 3-2-2, calculate the observable data and the sth parameter set iteration result C obs ,Θ s The weight variable w k The posterior probability p(w k |C obs ,Θ s ):

[0062]

[0063] Gaussian distribution for the last row Inverse gamma distribution IG(w k |υ s / 2,υ s / 2), generalized hyperbolic skewed Student's t distribution Take the logarithm, and when the parameter result of the sth iteration is known, divide the weight variable w at each moment by k The remaining parameters are constants:

[0064]

[0065]

[0066] Among them, const1, const2, and const3 are constants, which are related to the weight variable w at each moment. k Not relevant;

[0067] For the posterior probability p(w k |C obs ,Θ s ) and substitute the logarithm into the combined like terms to obtain:

[0068]

[0069]

[0070] After de-logarithmic processing, we get:

[0071]

[0072] Where const is the product of the logarithms of three constants const1, const2, and const3; according to the probability density function about w k The functional relationship can be obtained as follows: GIG(w k |λ,β,γ); GIG(w k |λ,β,γ) represents the generalized inverse Gaussian distribution, which obeys the parameters Weight variable w k , logw k Expected to be:

[0073]

[0074]

[0075]

[0076] Where E represents the random variable w k , logw k The expected calculation of <w k >、 <log w k > K a (b) represents the modified Bessel function of the second kind with order a and independent variable b;

[0077] After substituting the parameters λ, δ, and γ, we can calculate <w k >、 <log w k >:

[0078]

[0079]

[0080]

[0081] Step 3-2-3, according to the calculation <w k >、 , <log w k >, bring into Q-function calculation:

[0082]

[0083] Where:

[0084]

[0085] in It is independent of the parameter set Θ and the function P(υ) is only related to the degree of freedom parameter υ.

[0086] In step 3-2, in step 3-3, step M-step includes:

[0087] Step 3-3-1: To maximize the Q-function, find the partial derivative formula of the Q-function with respect to the parameters θ, β, Λ, and υ, and set the partial derivative formula to 0:

[0088]

[0089]

[0090]

[0091]

[0092] Where, ψ(x)=dlog(Γ(x)) / dx is the Digamma function;

[0093] The estimated parameters of θ, β, Λ, and υ in the parameter set Θ are θ s+1 , β s+1 ,Λ s+1 ,υ s+1 ;

[0094] Step 3-3-2: In step M-step, the model parameters of the system are estimated by iterating step 3-3-1, and the final iteration result is recorded as Θ * ={θ * , β * ,Λ * ,υ *}, θ * , β * ,Λ * ,υ * They represent the optimal model parameters, optimal skew parameters of statistical noise, optimal scale parameters and optimal degrees of freedom parameters respectively;

[0095] The convergence conditions are:

[0096] The present invention also provides a storage medium storing a computer program or instruction. When the computer program or instruction is executed, the ARX model robust identification method based on skew measurement noise is implemented.

[0097] Beneficial effect: Considering that the collected output data set is interfered by skew noise and a certain proportion of outliers, the present invention proposes a robust identification method of the ARX model based on skew measurement noise under the framework of the generalized expectation-maximization algorithm, which can effectively solve the model parameter identification problem in such cases, so as to better fit the actual system, and has been verified on the mass-spring damping system, which has greatly promoted the improvement of the ARX model system identification theory and its practical application. BRIEF DESCRIPTION OF THE DRAWINGS

[0098] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments, and the above and / or other advantages of the present invention will become more clear.

[0099] Figure 1 Flow chart of the method of the present invention.

[0100] Figure 2 This is the schematic diagram of the mass-spring-damper system.

[0101] Figure 3 Identification data collected for a mass-spring-damper system.

[0102] Figure 4 The simulated output curves of the method of the present invention and the comparison algorithm are shown in Figure 1 when the output data of the example of the present invention is disturbed by the GHSkewt distribution noise with scale, skew and degree of freedom parameters of (0.001, -0.2, 60).

[0103] Figure 5 It is the simulated output error curve of the method of the present invention and the comparison algorithm when the output data of the example of the present invention is disturbed by the GHSkewt distribution noise with scale, skew and degree of freedom parameters of (0.001, -0.2, 60). DETAILED DESCRIPTION

[0104] like Figure 1 As shown, the present invention provides an ARX model robust identification method based on skew measurement noise, comprising the following steps:

[0105] Step 1: Collect identification data according to the actual process of the mass-spring-damper system, and construct an ARX model of the mass-spring-damper system corresponding to the identification data;

[0106] Step 2: Based on the ARX model constructed in step 1, an ARX robust probability model based on generalized hyperbolic skewed Student's t distribution is established;

[0107] Step 3: Based on the generalized expectation maximization algorithm, the ARX robust probability model based on the generalized hyperbolic skewed Student's t distribution established in step 2 is identified. When the iteration converges, the optimal model parameter θ is obtained. * and the optimal skew parameter β for statistical noise * , optimal degrees of freedom parameter υ * and the optimal scale parameter Λ * .

[0108] Wherein, step 1 comprises:

[0109] The dynamic equations of the actual variables of the mass-spring-damper system are written as:

[0110]

[0111] where x(t) represents the position of the mass block, k s (t) is the stiffness of the spring, c d represents the damping coefficient of the fixed damper, and F(t) is the external force applied to the mass m. In the actual system, the external force F(t) represents the input u(t) of the mass-spring damping system, the mass block position x(t) represents the output y(t) of the mass-spring damping system, and N pairs of input and output data of the mass-spring damping system are collected at a time interval of T = 50ms, where N represents the length of the data moment;

[0112] The following ARX model is constructed based on the external force input u(t) and position output y(t) of the mass-spring-damper system:

[0113] A(z -1 )y k =B(z -1 ) k +e k

[0114]

[0115]

[0116] where u k represents the input of the mass-spring-damper system, y k represents the output of the mass-spring-damper system, e k represents the measurement noise of the mass spring damping system output, the subscript k represents the identification data time, from 1 to N; A(z -1 ) and B(z -1 ) is the unit shift operator z -1 The polynomial of z-1 y k =y k-1 , the superscript of z represents the magnitude of the backward shift moment, for example and A(z -1 ) and B(z -1 ) is the parameter corresponding to the polynomial, the subscript represents the corresponding backshift operator z term, which is the unknown parameter to be identified; n a and n b They are A(z -1 ) and B(z -1 ), indicating the known number of parameters of the polynomial.

[0117] Rewrite the above system into a linear regression form, define the parameter vector θ and the input and output information vector for:

[0118] θ=[a1,…,a na ,b1,…,b nb ] T

[0119]

[0120] Where T represents the matrix transposition operator, and the input and output information vectors The subscripts k-1 to kn of the input and output a , k-1 to kn b Indicates the moment;

[0121] The final linear regression form of the ARX model is:

[0122]

[0123] in, is the model parameter to be identified, e k is the measured noise of the output.

[0124] Step 2 includes:

[0125] In many existing recognition methods, it is usually assumed that the model noise e k Obey Gaussian distribution. These methods are effective in some cases, and when the system noise satisfies Gaussian estimation, these methods can provide unbiased parameter estimates. However, when the process data contains not only Gaussian noise but also outliers, the performance of these methods may degrade. Therefore, in order to solve the outlier problem, a noise model based on heavy-tailed distribution can be used to deal with the situation where the data set contains outliers. However, when the collected observation data is interfered by skewed noise, the above methods of noise modeling based on Gaussian distribution or heavy-tailed distribution often cannot obtain accurate parameter estimates, and thus cannot fit the actual system well.

[0126] Step 2-1: Noise e of the system model in linear regression form k Obey the generalized hyperbolic skew Student's t distribution:

[0127] p(e k |β,Λ,υ)=GHSkewt(e k |0,β,Λ,υ)

[0128] Where p(e k |β,Λ,υ) represents the measurement noise e under β,Λ,υ parameters k The posterior probability of GHSkewt is , GHSkewt represents the generalized hyperbolic skew Student's t distribution, β, Λ, υ represent the skew parameter, scale parameter and degree of freedom parameter of GHSkewt distribution respectively; in the above formula, the location parameter of this distribution is set to 0;

[0129] Step 2-2: Introduce weight variable w k Decompose the noise distribution in step 2-1:

[0130]

[0131] Where p(e k |w k ,β,Λ,υ) indicates that k ,β,Λ,υ parameter measurement noise e k The posterior probability, p(w k |υ) represents the weight variable w k The posterior probability under the scale parameter υ follows the following distribution:

[0132] p(w k |υ)=IG(w k |υ / 2,υ / 2)

[0133] Among them, IG(w k |υ / 2,υ / 2) represents the inverse gamma distribution, υ represents the degree of freedom parameter;

[0134] When introducing the weight variable w k Down:

[0135] p(e k |w k ,β,Λ,υ)=N(e k |βw k ,Λw k )

[0136] Where N(e k |βw k ,Λw k ) represents Gaussian distribution, βw krepresents the mean of the Gaussian distribution, Λw k represents the variance of the Gaussian distribution;

[0137] Step 2-3: Combine the system model in the form of linear regression Get the system output y k The generalized hyperbolic skew Student's t distribution, that is, the marginal distribution, is:

[0138]

[0139] in can be regarded as a sub-Gaussian distribution, represents the mean of the sub-Gaussian, Λw k represents the variance of each sub-Gaussian; IG(w k |υ / 2,υ / 2) represents the inverse gamma distribution, w k The weight of each sub-Gaussian distribution can be expressed;

[0140] By y k It can be seen from the marginal distribution expression that the GHSkewt distribution can be regarded as a convex combination of infinite sub-Gaussian distributions with different variances, and the introduction of skew and degree of freedom parameters in this distribution can flexibly describe the skew noise, which makes the GHSkewt distribution effectively deal with the negative impact of skew noise on the system.

[0141] The parameter set for constructing the ARX model of the mass-spring-damper system is Θ = {θ, β, Λ, υ}.

[0142] Step 3 includes:

[0143] Step 3-1, initialization: Initialize the parameter Θ to be estimated 1 , the number of iterations s is set to 1;

[0144] Step 3-2, execute step E-step: calculate and derive the Q-function Q(Θ|Θ s ):

[0145]

[0146] Where logp(C mis ,C obs |Θ) is the complete data set C={C mis ,C obs}, Θ s is the model parameter of the current linear system, C obs = {u 1:N ,y 1:N} is the observed data set, u 1:N Represents the input data of the system at each moment, y 1:NRepresents the input data of the system at each moment, C mis ={w 1:N} is a missing data set, w 1:N is the weight sequence, Θ = {θ, β, Λ, υ} is the model parameter; represents log p(C mis ,C obs |Θ) conditional expectation;

[0147] Step 3-3, M-step: By maximizing the Q-function, we get a new parameter estimate Θ s+1 , so that the parameter Θ s+1 Satisfy the constraints:

[0148] Q(Θ s+1 |Θ s )≥Q(Θ s |Θ s )

[0149] Let: Θ s =Θ s+1 ;

[0150] Step 3-4, let s = s + 1, return to step 3-2, and repeat the E-step and M-step until the parameters to be estimated converge.

[0151] The convergence conditions are:

[0152] In step 3-2, step E-step includes:

[0153] Step 3-2-1, calculate the complete data set C = {C mis ,C obs The log-likelihood function logp(C mis ,C obs |Θ), using the probability chain rule, the log-likelihood function is written as:

[0154]

[0155]

[0156] Where C1=p(u 1:N |Θ) is a constant term;

[0157] The log-likelihood function is applied to the weight sequence w 1:N To obtain the conditional expectation, we can get:

[0158]

[0159] Where C1=p(u 1:N |Θ) is a constant term, Θs ={θ s ,β s ,Λ s ,υ s} represents the sth iteration result of model parameters, skew parameters, scale parameters and degrees of freedom parameters on the weight sequence w 1:N The integral term for conditional expectation can be expressed as:

[0160]

[0161]

[0162]

[0163] Step 3-2-2, calculate the observable data and the sth parameter set iteration result C obs ,Θ s The weight variable w k The posterior probability p(w k |C obs ,Θ s ):

[0164]

[0165] Gaussian distribution of the last line of the above equation Inverse gamma distribution IG(w k |υ s / 2,υ s / 2), generalized hyperbolic skewed Student's t distribution Take the logarithm, and when the parameter result of the sth iteration is known, divide the weight variable w at each moment by k The remaining parameters are constants:

[0166]

[0167]

[0168]

[0169] Among them, const1, const2, and const3 are constants, which are related to the weight variable w at each moment. k Not relevant;

[0170] According to the above formula, the posterior probability p(w k |C obs ,Θ s ) and substitute it into the combined like terms to obtain:

[0171]

[0172]

[0173] De-logarithmize the above formula to get:

[0174]

[0175] Where const is the product of the logarithms of three constants const1, const2, and const3; according to the probability density function about w k The functional relationship of GIG(w k |λ,β,γ); GIG(w k |λ,β,γ) represents the generalized inverse Gaussian distribution, which obeys the parameters Weight variable w k , logw k Expected to be:

[0176]

[0177]

[0178]

[0179] Where E represents the random variable w k , logw k The expected calculation of <w k >、 <log w k >, K a (b) represents the second kind modified Bessel function with order a and independent variable b, and its specific calculated value is as shown in the above formula;

[0180] Substituting the parameters λ, δ, and γ into the above formula, we can get <w k >、 <log w k >:

[0181]

[0182]

[0183]

[0184] Step 3-2-3, according to the calculation <w k >、 <log w k >Substitute into Q-function calculation:

[0185]

[0186]

[0187] in It is independent of the parameter set Θ, and the function P(υ) is only related to the degree of freedom parameter υ;

[0188] In step 3-3, step M-step includes:

[0189] Step 3-3-1: To maximize the Q-function, find the partial derivative formula of the Q-function with respect to the parameters θ, β, Λ, and υ, and set the partial derivative formula to 0:

[0190]

[0191]

[0192]

[0193]

[0194] Among them, ψ(x)=dlog(Γ(x)) / dx is the digamma function.

[0195] Solving the above nonlinear equations yields the estimated θ of the parameters θ, β, Λ, and υ in the parameter set Θ. s+1 , β s+1 ,Λ s+1 ,υ s+1 ;

[0196] Step 3-3-2: In the M-step, the model parameters of the system are estimated by iterating through step 3-1, and the final iteration results are recorded as: Θ * ={θ * , β * ,Λ * ,υ *}, respectively represent the optimal model parameters θ * and the optimal skew parameter β for statistical noise * , optimal degrees of freedom parameter υ * and the optimal scale parameter Λ * ;

[0197] The convergence conditions are:

[0198] At this point, the entire process of parameter identification of the algorithm of the present invention is completed.

[0199] Embodiment 1:

[0200] In this embodiment, the ARX model robust identification method based on skew measurement noise is specifically implemented according to the following steps:

[0201] Step 1: Construct observation data sets separately: C obs and missing data set: C mis ;

[0202] Step 2: Initialize parameter Θ 1 , and set s=1;

[0203] Step 3: Update the posterior probability distribution and expectation:

[0204] Using the formula given in E-step, calculate the weight sequence w 1:N The posterior probability p(w k |C obs ,Θ s ) and the expectation at each moment k <w k >、 <log w k >, reconstruct the Q-function;

[0205] Step 4: Update model parameter estimates:

[0206] Using the formula given in M-step, update the model parameters θ s+1 and the noise parameter β s+1 ,Λ s+1 ,υ s+1 .

[0207] Step 5: Set s = s + 1 until the iteration converges and obtain the optimal model parameter and noise parameter set θ * ={θ * , β * ,Λ * ,υ *}.

[0208] The effectiveness of the present invention is verified by the following examples:

[0209] Physical simulation example:

[0210] Use Figure 2 The effectiveness of the method of the present invention is verified by the principle of the mass spring damping system. The system mainly includes a mass block, a spring and a damper, and its dynamic equation is described as:

[0211]

[0212] where x(t) represents the position of the mass block, k s (t) is the stiffness of the spring, c dRepresents the damping coefficient of the fixed damper, and F(t) is the external force applied to mass m. In the actual system, the external force F(t) represents the input u(t) of the mass spring damping system, the mass block position x(t) represents the output y(t) of the mass spring damping system, and the input and output data of the mass spring damping system at 1000 moments are collected at a time interval of T = 50ms, as shown in Figure 3 ;

[0213] Two methods using Laplace distribution and Student's t distribution modeling are selected as comparison algorithms, respectively named comparison algorithm 1 and comparison algorithm 2. There is skew noise in the output measurement of the spring damping system. The dynamic output of the mass spring damping system is predicted using the method of the present invention and the comparison method 2.

[0214] Based on the method of the present invention, the simulation prediction results of the dynamic position output of the spring damping system are shown in Table 1, where the root mean square error (RMSE) and the coefficient of determination (R 2 ), defined as follows:

[0215]

[0216]

[0217] Where N is the length of the identification data, and Represent the true output at the nth sampling moment, the estimated output, and the average of the true output at all moments. As can be seen from Table 1, the output estimation of the method of the present invention has the smallest root mean square error (RMSE), and the coefficient of determination (R 2 ) is closer to 1, which indicates that the present invention has the best model output estimation;

[0218] Figure 4 The simulated output curve of the mass-spring damping system is plotted, the black line represents the real output of the system, the * line represents the estimated output of the method of the present invention, and the × line and the + line represent two comparative methods. It is not difficult to find that the * line almost completely overlaps with the black line, while the comparative algorithm cannot, which proves that the algorithm of the present invention can accurately capture the dynamic characteristics of the mass-spring damping system. Figure 5 The error discrete points of the simulated output at each moment are plotted, the black dotted line represents the same as the actual output, the * line represents the simulated output error discrete points of the method of the present invention, the × line and the + line represent the simulated output error discrete points of the comparison algorithm, respectively. It can be clearly seen from the comparison results that the simulated output error curve of the method of the present invention is closer to 0 and the curve fluctuation is small, while the comparison algorithm has a large deviation and large fluctuation, which further verifies that the method of the present invention can better capture the dynamic characteristics of the mass spring damping system.

[0219] Table 1

[0220]

[0221]

[0222] Simulation results explanation:

[0223] The present invention can accurately capture the dynamic characteristics of the mass spring damping system under the interference of skew noise and has better identification accuracy, while the comparison algorithm does not have this ability.

[0224] In a specific implementation, the present application provides a computer storage medium and a corresponding data processing unit, wherein the computer storage medium can store a computer program, and when the computer program is executed by the data processing unit, the invention content of the robust identification method of an ARX model based on skew measurement noise provided by the present invention and some or all steps in each embodiment can be executed. The storage medium can be a disk, an optical disk, a read-only memory (ROM) or a random access memory (RAM), etc.

[0225] Those skilled in the art can clearly understand that the technical solutions in the embodiments of the present invention can be implemented by means of computer programs and their corresponding general hardware platforms. Based on this understanding, the technical solutions in the embodiments of the present invention are essentially or partly contributed to the prior art can be embodied in the form of a computer program, i.e., a software product, which can be stored in a storage medium and includes several instructions for enabling a device including a data processing unit (which can be a personal computer, a server, a single-chip microcomputer, a MUU or a network device, etc.) to execute the methods described in various embodiments of the present invention or certain parts of the embodiments.

[0226] The present invention provides a robust identification method for an ARX model based on skew measurement noise. There are many methods and approaches to implement the technical solution. The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the protection scope of the present invention. All components not specified in this embodiment can be implemented by existing technologies.

Claims

1. A robust identification method for ARX model based on skew measurement noise, characterized in that: The following steps are involved: Step 1: Collect identification data according to the actual process of the mass-spring-damper system, and construct the ARX model of the mass-spring-damper system corresponding to the identification data: Step 2: Based on the ARX model constructed in step 1, an ARX robust probability model based on generalized hyperbolic skewed Student's t distribution is established; Step 3: Based on the generalized expectation maximization algorithm, the ARX robust probability model established in step 2 is identified. When the iteration converges, the final iteration result is used as the optimal estimation parameter; Step 1 includes: The dynamic equations of the actual variables of the mass-spring-damper system are written as: Where x(t) represents the position of the mass block, k s (t) is the stiffness of the spring, c d represents the damping coefficient of the fixed damper, F(t) is the external force applied to mass m; In the actual system, the external force F(t) represents the input u(t) of the mass-spring-damper system, the mass block position x(t) represents the output y(t) of the mass-spring-damper system, and N pairs of input and output data of the mass-spring-damper system are collected at time intervals T, where N represents the length of the data moment; The following ARX model is constructed based on the input u(t) and output y(t) of the mass-spring-damper system: A(z -1 )y k =B(z -1 )u k +e k where u k represents the input of the mass-spring-damper system, y k represents the output of the mass-spring-damper system, e k represents the measurement noise of the mass spring damping system output, the subscript k represents the identification data moment, and the value ranges from 1 to N; A(z -1 ) and B(z -1 ) is the unit shift operator z -1 The polynomial of z -1 y k =y k-1 , the superscript z represents the magnitude of the backward movement moment, a1, a2, …, a na and Respectively represent A(z -1 ) and B(z -1 ) is the parameter corresponding to the polynomial, the subscript represents the corresponding backshift operator z term, which is the unknown parameter to be identified; n a and n b They are A(z -1 ) and B(z -1 ) is a known order, which indicates the known number of polynomial parameters; In step 1, the ARX model is rewritten as a linear regression model, defining the parameter vector θ and the input and output information vector for: θ=[a1,…,a na ,b1,…,b nb ] T Where T represents the matrix transposition operator, and the input and output information vectors The subscripts k-1 to kn of the input and output a , k-1 to kn b Indicates the moment; The final linear regression system model is: In step 2, an autoregressive ARX robust probability model with external input based on generalized hyperbolic skewed Student's t distribution is established, specifically including: Step 2-1: Measurement noise e at the output at time k k Obey the generalized hyperbolic skew Student's t distribution: p(e k |β,Λ,υ)=GHSkewt(e k |0,b,L,u) Where p(e k |β,Λ,υ) represents the measurement noise e under β,Λ,υ parameters k The posterior probability of , GHSkewt represents the generalized hyperbolic skew Student's t distribution, β, Λ, υ represent the skew parameter, scale parameter and degree of freedom parameter of GHSkewt distribution respectively; Step 2-2: Introduce weight variable w k Decompose the noise distribution in step 2-1: Where p(e k |w k ,β,Λ,υ) indicates that k ,β,Λ,υ parameter measurement noise e k The posterior probability, p(w k |υ) represents the weight variable w k The posterior probability under the scale parameter υ follows the following distribution: p(in k |υ)=IG(in k |υ / 2,υ / 2) Among them, IG(w k |υ / 2,υ / 2) represents the inverse gamma distribution, υ represents the degree of freedom parameter; When introducing the weight variable w k Down: p(e k |w k ,β,Λ,υ)=N(e k |bw k ,Λw k ) Where N(e k |βw k ,Λw k ) represents Gaussian distribution, βw k represents the mean of the Gaussian distribution, Λw k represents the variance of the Gaussian distribution; Step 2-3: Combine the system model in the form of linear regression Get the system output y k The generalized hyperbolic skew Student's t distribution is: in is a sub-Gaussian distribution, represents the mean of the sub-Gaussian, Λw k represents the variance of each sub-Gaussian; The parameter set for constructing the ARX model of the mass-spring-damper system is Θ = {θ, β, Λ, υ}.

2. The method according to claim 1, characterized in that Step 3 includes: Step 3-1, initialization: Initialize the parameter Θ to be estimated 1 , the number of iterations s is set to 1; Step 3-2, execute step E-step: calculate and derive the Q-function Q(Θ|Θ s ): Where logp(C mis ,C obs |Θ) is the complete data set C={C mis ,C obs }, Θ s Estimated parameter set for the model of the current linear system, Cobs = {u1: N ,y1: N } is the observed data set, u 1:N Represents the input data of the system at each moment, y 1:N Represents the output data of the system at each moment, Cmis={w1: N } is a missing data set, w 1:N is the weight sequence, Θ = {θ, β, Λ, υ} is the model parameter set; represents log p(C mis ,C obs |Θ) conditional expectation; Step 3-3, execute step M-step: by maximizing the Q-function, obtain a new parameter estimate Θ s+1 , so that the parameter Θ s+1 Satisfy the constraints: Q(Θ s+1 |I s )≥Q(Θ s |I s ) Let: Θ s =Θ s+1 ; Step 3-4, let s=s+1, return to step 3-2, and repeat step E-step and step M-step until the parameters to be estimated converge.

3. The method according to claim 2, characterized in that In step 3-2, step E-step includes: Step 3-2-1, calculate the log-likelihood function logp(C mis ,C obs |Θ), using the probability chain rule, the log-likelihood function is written as: Where C1=p(u 1:N |Θ) is a constant term; The log-likelihood function is applied to the weight sequence w 1:N Taking the conditional expectation, we get: Where C1=p(u 1:N |Θ) is a constant term, Θ s ={θ s ,β s ,Λ s ,υ s } represents the result of the sth parameter iteration, θ s ,β s ,Λ s ,υ s They represent the sth iteration results of model parameters, skew parameters, scale parameters and degrees of freedom parameters, respectively, for the weight sequence w 1:N Find the integral term of the conditional expectation: Step 3-2-2, calculate the observable data and the sth parameter set iteration result C obs ,Θ s The weight variable w k The posterior probability p(w k |C obs ,Θ s ): Gaussian distribution for the last row Inverse gamma distribution IG(w k |υ s / 2,υ s / 2), generalized hyperbolic skewed Student's t distribution Take the logarithm, and when the parameter result of the sth iteration is known, divide the weight variable w at each moment by k The remaining parameters are constants: Among them, const1, const2, and const3 are constants, which are related to the weight variable w at each moment. k Not relevant; For the posterior probability p(w k |C obs ,Θ s ) and substitute the logarithm into the combined like terms to obtain: After de-logarithmic processing, we get: Where const is the product of the logarithms of three constants const1, const2, and const3; according to the probability density function about w k The functional relationship can be obtained as follows: GIG(w k |λ,β,γ); GIG(w k |λ,β,γ) represents the generalized inverse Gaussian distribution, which obeys the parameters Weight variable w k , logw k Expected to be: Where E represents the random variable w k , logw k The expected calculation of <w k >、 <w k -1 >、 <log w k > K a (b) represents the modified Bessel function of the second kind with order a and independent variable b; After substituting the parameters λ, δ, and γ, we can calculate <w k >、 <log w k >: Step 3-2-3, according to the calculation <w k >、 <log w k >, bring into Q-function calculation: Where: in It is independent of the parameter set Θ and the function P(υ) is only related to the degree of freedom parameter υ.

4. The method according to claim 3, characterized in that In step 3-2, in step 3-3, step M-step includes: Step 3-3-1: To maximize the Q-function, find the partial derivative formula of the Q-function with respect to the parameters θ, β, Λ, and υ, and set the partial derivative formula to 0: Where, ψ(x)=dlog(Γ(x)) / dx is the Digamma function; The estimated parameters of θ, β, Λ, and υ in the parameter set Θ are θ s+1 , β s+1 ,Λ s+1 ,υ s+1 ; Step 3-3-2: In step M-step, the model parameters of the system are estimated by iterating step 3-3-1, and the final iteration result is recorded as Θ * ={θ * , β * ,Λ * ,υ * }, θ * , β * ,Λ * ,υ * They represent the optimal model parameters, optimal skew parameters of statistical noise, optimal scale parameters and optimal degrees of freedom parameters respectively; The convergence conditions are:

5. A storage medium, characterized in that: A computer program or instruction is stored, and when the computer program or instruction is executed, the method according to any one of claims 1 to 4 is implemented.

Images