A robust identification method for nonlinear state-space systems based on thick-tailed noise

Abstract

The present invention provides a method for robust identification of nonlinear state space system based on thick-tailed noise, including: step 1: establishing a continuous-time dynamic model of a two-link manipulator system, and then converting it into a discrete state space model, collecting the vertical downward external force data of the end of the two-link manipulator system, and collecting the angular displacement of the two links relative to the horizontal plane to form an identification data set; step 2: introducing the system noise distribution and the prior distribution of the model parameters to form a probabilistic description of the two-link manipulator system identification problem; step 3: based on the identification data set and the probabilistic description of the two-link manipulator system identification problem, iteratively updating under the framework of the expectation maximization algorithm to obtain the parameters to be identified. The present invention can effectively solve the model parameter identification problem of such a system, and has greatly promoted the improvement of the nonlinear state space system identification theory of the two-link manipulator system and its industrial practical application.

Description

Technical Field

[0001] The present invention relates to the technical field of industrial system modeling and industrial system model parameter identification, and in particular to a thick-tail noise-based robust identification method for a nonlinear state space system of a double-link mechanical arm system. Background Art

[0002] In recent years, most factories at home and abroad have achieved automated or semi-automated production, which has greatly improved production efficiency and production costs. However, the double-link robot is nonlinear and time-varying. The double-link robot is a multi-input and multi-output system. The Lagrangian method can be used to establish a mathematical model for it. Since there are a lot of uncertain factors in the preliminary dynamic modeling in practical applications, the nonlinear state space model can not only well characterize the internal dynamics of the system, but also facilitate the design of the system controller. Therefore, it is necessary to study the nonlinear state space model of the double-link robot system.

[0003] Since industrial systems inevitably have problems such as sensor failure and external interference, there are often outliers in the collected process output data, and there may also be unknown disturbances that cannot be observed inside the system. If the role of outliers in the state update process and output measurement process is considered in the identification of the nonlinear state space model of the double-link robot arm system, the identification accuracy of the system model will often be reduced. Summary of the invention

[0004] Purpose of the invention: In order to solve the problem that in the prior art, when there are abnormal values ​​in the state update process and the output measurement process in the dual-link manipulator system, the system identification accuracy is reduced, the present invention provides a robust identification method for a nonlinear state space system based on thick-tailed state noise, comprising the following steps:

[0005] The following steps are involved:

[0006] Step 1: Establish a continuous-time dynamic model of the two-link manipulator system, convert the continuous-time dynamic model into a discrete state-space model after finite difference, collect the vertical downward external force data of the end of the two-link manipulator system, and collect the angular displacement of the two links relative to the horizontal plane to form an identification data set;

[0007] The discrete state space model includes a state transfer equation and a measurement equation;

[0008] Step 2: For the discrete state space model, introduce the system noise distribution and the prior distribution of model parameters to form a probabilistic description of the dual-link manipulator system identification problem;

[0009] Step 3: Based on the identification data set and the probabilistic description of the two-link robotic arm system identification problem, iterative updates are performed under the framework of the expectation maximization algorithm to obtain the parameters to be identified.

[0010] In step 1, the continuous-time dynamics model of the two-link manipulator system is expressed as:

[0011]

[0012] Among them, M, E are coefficient matrices; K(t) is a coefficient matrix, some elements of which change with time t; Φ is the angular displacement of the two joints, are the angular velocities of the two joints, are the two joint angular accelerations; and M, E, K(t) are specifically expressed as:

[0013]

[0014]

[0015] Where m and l represent the mass and length of the two connecting rods respectively; g represents the acceleration of gravity; in order to simplify the system, the damping of the system is assumed to be concentrated at the nodes, d1 and d2 represent the damping coefficients of the two nodes; the connecting rod is assumed to be rigid, the elasticity of joint 1 and joint 2 is represented by two constant rotation springs, k1 and k2 represent the spring stiffness; α 10 ,α 20 They respectively represent the steady-state positions of the angular displacements of the two connecting rods relative to the horizontal plane; F(t) represents an external force applied to the end of the second connecting rod, which changes with time t and maintains an angle of α3 with the horizontal plane.

[0016] In step 1, the state transfer equation is expressed as:

[0017] x(k)=x(k-1)+Δ t Ax(k-1)

[0018] Which set is the state quantity of the system, T represents the matrix transpose, x(k) and x(k-1) represent the state of the system at time k and time k-1 respectively; Δ t is the sampling time; A is the coefficient matrix, which is specifically expressed as:

[0019]

[0020]

[0021]

[0022] Among them, θ1, θ2 are the model parameters to be identified; b1, b2 are intermediate parameters.

[0023] In step 1, the measurement equation is expressed as:

[0024]

[0025] Wherein, y(k) represents the load end data obtained by measurement.

[0026] In step 1, F(t) is expressed as:

[0027] F(t)=10sin(0.02πt).

[0028] In step 1, the discrete state space model is expressed as:

[0029] x k =f(x k-1 ,u k-1 ,θ f )+w k

[0030] y k =g(x k )+v k

[0031] Among them, x k =x(k),x k-1 =x(k-1),x k 、x k-1 They represent the state of the two-link manipulator system at time k and time k-1, namely the angular velocity and angular displacement of the two nodes; u k-1 =F((k-1)Δt),u k-1 represents the external force input of the two-link manipulator system at time k-1, that is, the vertical downward external force at the end of the rod; y k =y(k),y k represents the output data collected by the dual-link robotic arm system, i.e., the angular displacement of the two nodes; f(x k-1 ,u k-1 ,θ f )=x k-1 +Δ t Ax k-1 is the nonlinear mapping function of the state update process, θ f =[θ1,θ2] is the parameter to be identified in the nonlinear mapping function; g(x k ) is the output measurement function, w k ,v k are the state noise and measurement noise of the two-link robotic arm system, respectively.

[0032] Step 2 includes:

[0033] Step 2-1: Set w k ,v k Set to zero-mean Gaussian white noise, obeying the distribution N(w k |0,R w ),N(v k |0,R v );where R v is the scale parameter;

[0034] The state transition local noise change parameter is R w / γ k ,get:

[0035] w k |(R w ,γ k )~N(w k |0,R w / γ k )

[0036] Where R w is the scale parameter, N(w k |0,R w / γ k ) represents the noise sequence w k It obeys Gaussian distribution, with mean 0 and variance R w / γ k ; γ k is the weight parameter:

[0037] γ k |r w ~Gamma(γ k |r w / 2,r w / 2)

[0038] where r w represents the degree of freedom parameter, x k It obeys the following Gaussian distribution:

[0039] x k |(f(x k-1 ,u k-1 ,θ f ),R w ,γ k )~N(w k |(f(x k-1 ,u k-1 ,θ f ),R w / γ k ));

[0040] Step 2-2: x k The marginal distribution ofk ,γ k |f(x k-1 ,u k-1 ,θ f ),R w ,r w ) for γ k Integrating gives:

[0041]

[0042] where Γ(·) represents the Gamma function, p(γ k |r w ) is the probability density function that follows the Gamma distribution;

[0043] Step 2-3: Set the output measurement local noise change parameter to R v / λ k ,get:

[0044] v k |(R v ,λ k )~N(v k |0,R v / λ k )

[0045] Where N(v k |0,R v / λ k ) represents the noise sequence v k It obeys Gaussian distribution, with mean 0 and variance R v / λ k ; k is the weight parameter:

[0046] λ k |r v ~Gamma(λ k |r v / 2,r v / 2)

[0047] where v v represents the degree of freedom parameter, y k The Gaussian distribution is as follows:

[0048] y k |(g(x k ),R v ,λ k )~N(v k |g(x k ),R v / λ k );

[0049] Step 2-4:k The marginal distribution of k ,λ k |g(x k ),R v ,r v ) for δ k Integrating gives:

[0050]

[0051] Step 3 (The identification data set collected in step 1 is used as the observation data set C in step 3 obs Using, then obtaining unknown parameters in the system by the method of the present invention) includes:

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

[0053] Step 3-2, execute step E-step, calculate and derive the Q-function Q(Θ|Θ) of the two-link manipulator system s ):

[0054]

[0055] Where logp(C mis ,C obs |Θ) is the complete data set C={C mis ,C obs}'s log-likelihood function; Θ s is the parameter of the current nonlinear state space model; C obs = {u 1:N ,y 1:N} is the identification dataset of the two-link robotic arm system collected in step 1, u 1:N represents the external force input at the end of the link of the double-link robotic arm system, where u k =F(kΔt),Δ t is the sampling time, y 1:N represents the angular displacement of the two links of the double-link robotic arm system relative to the horizontal plane; C mis ={x 1:N ,γ 1:N ,λ 1:N} is a missing data set, that is, a data set that cannot be observed, where x 1:N represents the state quantity of the system, γ 1:N ,λ 1:N is the weight sequence, Θ={θ f ,R w ,r w ,R v ,r v} are the model parameters to be identified; represents log p(C mis ,C obs |Θ) conditional expectation;

[0056] Step 3-3, execute step M-step, and search for a new parameter estimate Θ by optimizing the search method. s+1 , so that the parameter Θ s+1 The following constraints are met:

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

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

[0059] Step 3-4, let s = s + 1, return to step 3-2, and repeat step E-step and step M-step until the estimated parameters converge. The convergence condition is:

[0060] In step 3-2, the step E-step includes:

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

[0062]

[0063] The log-likelihood function is respectively applied to the state quantity x of the system at time k k , weight sequence γ 1:N and the weight sequence λ 1:N Taking the conditional expectation, we get:

[0064]

[0065] Step 3-2-2, first calculate the state x of the system at time k k The posterior probability and Q(Θ|Θ s ) formula, and use the particle filter algorithm to establish the posterior probability distribution of the system state quantity, and obtain:

[0066]

[0067] where δ(·) is the Dirac function, represents the lth particle at time k, represents the weight of the lth particle at time k;

[0068] According to all known observation data, the state quantity x of the system at time k is obtained by particle smoothing. k The posterior probability distribution of :

[0069]

[0070] in,

[0071] Get the state x of the system at time k and time k-1 k ,x k-1 The joint posterior probability distribution of is:

[0072]

[0073] in,

[0074] The sum of the first integral term and the second integral term in the Q-function is:

[0075]

[0076] Where I1 is the sum of the first integral term and the second integral term in the Q-function, It is expressed as the particle weight after particle smoothing;

[0077] γ k and logγ k The expectation is expressed as <λ k > l and <logλ k > l , and there are:

[0078]

[0079]

[0080] in is the degree of freedom parameter of the measurement noise at the sth iteration, is the scale parameter of the measurement noise at the sth iteration;

[0081] The sum of the third and fourth integral terms in the Q-function is:

[0082]

[0083] Among them, I2 is the sum of the third integral term and the fourth integral term in the Q-function, It is represented as the joint weight of particles after particle smoothing; γ k and logγ k The expectation is expressed as <γ k > l and <logγ k > l , and there are:

[0084]

[0085]

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

[0087] Then the Q-function is expressed as:

[0088] Q(Θ|Θ s )=Q1(θ f ,R v )+Q2(R w )+Q3(r w )+Q4(r v )+C1

[0089] Where Q1(θ f ,R v ) represents the Q-function with respect to θ f ,R v Mathematical expression of Q2(R w ) represents the Q-function with respect to R w The mathematical expression of Q3(r w ) represents the Q-function with respect to r w Mathematical expression of Q4(r v ) represents the Q-function with respect to r v Mathematical expression of ; C1 represents the constant term, where:

[0090]

[0091]

[0092]

[0093]

[0094]

[0095] Where logp(u 1:N |Θ) is represented as the prior probability distribution of the input data, and logp(x1|γ1,Θ) is represented as the posterior probability distribution of the initial state.

[0096] Step 3-3, the step M-step includes:

[0097] Step 3-3-1, set Q1(θ f ,R w ) for R w Derivative, Q2(R v ) for R v Taking the derivative and setting it equal to zero, we get:

[0098]

[0099]

[0100] Step 3-3-2, Substituting into Q1(θ f ,R w )get:

[0101]

[0102] <λ k > l and <γ k > ij are the weights assigned to the square of the model prediction deviation corresponding to the state update process and the output measurement process data, respectively;

[0103] Using the damped Newton method to optimize the criterion function Q5(θ f ), and get the parameter θ f The optimized expression is:

[0104]

[0105] in, Denotes the parameter θ f The updated value for the s+1th time, μ represents the damping Newton constant, and I represents the identity matrix;

[0106] Step 3-3-3, Q3(r w ) w Derivative, Q4(r v ) v Taking the derivative and setting it equal to zero, we get:

[0107]

[0108]

[0109] Solve to get the parameter r w ,r v The estimates are

[0110] Step 3-3-4, through iterative cycle steps 3-3-1 to 3-3-3, obtain the model parameters of the system, and record the final iteration result as

[0111] The convergence conditions are:

[0112] Beneficial effect: Taking into account the existence of outliers in the state update process and output measurement process of the two-link robotic arm system, the method of the present invention proposes a robust identification method for nonlinear state-space systems based on thick-tailed noise under the framework of the generalized expectation-maximization algorithm. This method can effectively solve the model parameter identification problem of such systems, and has greatly promoted the improvement of the nonlinear state-space system identification theory of the two-link robotic arm system and its practical industrial application. BRIEF DESCRIPTION OF THE DRAWINGS

[0113] 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.

[0114] Figure 1 It is a double-link robotic arm system.

[0115] Figure 2 This is a comparison chart of system parameter estimation curves of the NSSM-Student method and the NSSM-Gaussian method when the proportion of abnormal values ​​in the state update process and output measurement process is 10%.

[0116] Figure 3 Flow chart of the method of the present invention. DETAILED DESCRIPTION

[0117] like Figure 1 , Figure 3 As shown, the present invention provides a robust identification method for nonlinear state space system based on thick-tail noise, comprising the following steps:

[0118] Step 1: Establish a continuous-time dynamic model of the double-link robotic arm system, convert the continuous-time dynamic model of the double-link robotic arm system into a discrete state-space model after finite difference, collect the vertical downward external force data of the end of the double-link robotic arm system, and collect the angular displacement of the two links relative to the horizontal plane to form an identification data set; the discrete state-space model includes a state transfer equation and a measurement equation;

[0119] Step 2: For the discrete state space model, introduce the system noise distribution and the prior distribution of model parameters to form a probabilistic description of the dual-link manipulator system identification problem;

[0120] Step 3: Based on the identification data set and the probabilistic description of the two-link robotic arm system identification problem, iterative updates are performed under the framework of the expectation maximization algorithm to obtain the parameters to be identified.

[0121] The continuous-time dynamics model of the two-link manipulator system is expressed as:

[0122]

[0123] Among them, M, E are coefficient matrices; K(t) is a coefficient matrix, some elements of which change with time t; Φ is the angular displacement of the two joints, are the angular velocities of the two joints, are the two joint angular accelerations. And M, E, K(t) are specifically expressed as:

[0124]

[0125]

[0126] Where m and l represent the mass and length of the two connecting rods respectively; g represents the acceleration of gravity; in order to simplify the system, the damping of the system is assumed to be concentrated at the nodes, d1 and d2 represent the damping coefficients of the two nodes; the connecting rod is assumed to be rigid. The elasticity of joints 1 and 2 is represented by two constant rotation springs, k1 and k2 represent the spring stiffness; α 10 ,α 20 They respectively represent the steady-state positions of the angular displacements of the two connecting rods relative to the horizontal plane; F(t) represents an external force applied to the end of the second connecting rod, which changes with time t and maintains an angle of α3 with the horizontal plane.

[0127] The state transfer equation of the dual-link manipulator system is expressed as:

[0128] x(k)=x(k-1)+Δ t Ax(k-1)

[0129] Which set is the state quantity of the system, T represents the matrix transpose, x(k) and x(k-1) represent the state of the system at time k and time k-1 respectively; Δ t is the sampling time; A is the coefficient matrix, which is specifically expressed as:

[0130]

[0131]

[0132]

[0133] Among them, θ1, θ2 are the model parameters to be identified; b1, b2 are intermediate parameters.

[0134] The measurement equation of the two-link robotic arm system is expressed as:

[0135]

[0136] Wherein, y(k) represents the load end data obtained by measurement.

[0137] The input of the two-link robot arm system is the vertical downward external force at the end:

[0138] F(t)=10sin(0.02πt)

[0139] The discrete state space model of the two-link manipulator system is expressed as:

[0140] x k =f(x k-1 ,u k-1 ,θ f )+w k

[0141] y k =g(x k )+v k

[0142] Among them, x k =x(k),x k-1 =x(k-1), which represent the state quantities of the two-link manipulator system at time k and time k-1, i.e., the angular velocity and angular displacement of the two nodes; u k-1 =F((k-1)Δt), which represents the external force input to the two-link manipulator system at time k-1, i.e., the vertical downward force at the end of the rod; y k =y(k), represents the output data collected by the dual-link robot system, i.e., the angular displacement of the two nodes; f(x k-1 ,u k-1 ,θ f )=x k-1 +Δ t Ax k-1 is a nonlinear mapping function, θ f =[θ1,θ2],θ f is the parameter to be identified in the nonlinear mapping function; w k ,v k are the state noise and measurement noise of the two-link robotic arm system, respectively.

[0143] w k ,v k Set to zero-mean Gaussian white noise, obeying the distribution N(w k |0,R w ),N(v k|0,R v );where R v is the scale parameter; step 2 includes:

[0144] Step 2-1: Change the state transfer local noise change parameter to R w / γ k ,get:

[0145] w k |(R w ,γ k )~N(w k |0,R w / γ k )

[0146] Where R w is the scale parameter, N(w k |0,R w / γ k ) represents the noise sequence w k It obeys Gaussian distribution, with mean 0 and variance R w / γ k ; γ k is the weight parameter:

[0147] γ k |r w ~Gamma(γ k |r w / 2,r w / 2)

[0148] where r w represents the degree of freedom parameter, x k It obeys the following Gaussian distribution:

[0149] x k |(f(x k-1 ,u k-1 ,θ f ),R w ,γ k )~N(w k |f(x k-1 ,u k-1 ,θ f ),R w / γ k );

[0150] Step 2-2: x k The marginal distribution of k ,γ k |f(x k-1 ,u k-1 ,θ f ),R w ,r w) for γ k Integrating gives:

[0151]

[0152] where Γ(·) represents the Gamma function, p(γ k |r w ) is the probability density function that follows the Gamma distribution;

[0153] Step 2-3: Set the output measurement local noise change parameter to R v / λ k ,get:

[0154] v k |(R v ,λ k )~N(v k |0,R v / λ k )

[0155] Where N(v k |0,R v / λ k ) represents the noise sequence v k It obeys Gaussian distribution, with mean 0 and variance R v / λ k ; k is the weight parameter:

[0156] λ k |r v ~Gamma(λ k |r v / 2,r v / 2)

[0157] where v v represents the degree of freedom parameter, y k The Gaussian distribution is as follows:

[0158] y k |(g(x k ),R v ,λ k )~N(v k |g(x k ),R v / λ k );

[0159] Step 2-4: k The marginal distribution of k ,λ k |g(x k ),R v ,r v ) for δk Integrating gives:

[0160]

[0161] Step 3 includes:

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

[0163] Step 3-2, execute step E-step, calculate and derive the Q-function Q(Θ|Θ) of the two-link manipulator system s ):

[0164]

[0165] Where logp(C mis ,C obs |Θ) is the complete data set C={C mis ,C obs}'s log-likelihood function; Θ s is the parameter of the current nonlinear state space model; C obs = {u 1:N ,y 1:N} is the identification dataset of the two-link robotic arm system collected in step 1, u 1:N represents the external force input at the end of the link of the double-link robotic arm system, where u k =F(kΔt),Δ t is the sampling time, y 1:N represents the angular displacement of the two links of the double-link robotic arm system relative to the horizontal plane; C mis ={x 1:N ,γ 1:N ,λ 1:N} is a missing data set, that is, a data set that cannot be observed, where x 1:N represents the state quantity of the system, γ 1:N ,λ 1:N is the weight sequence, Θ={θ f ,R w ,r w ,R v ,r v} are the model parameters to be identified; represents logp(C mis ,C obs |Θ) conditional expectation;

[0166] Step 3-3, execute step M-step, and search for a new parameter estimate Θ by optimizing the search method. s+1 , so that the parameter Θ s+1 The following constraints are met:

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

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

[0169] Step 3-4, let s = s + 1, return to step 3-2, and repeat step E-step and step M-step until the estimated parameters converge. The convergence condition is:

[0170] In step 3-2, the step E-step includes:

[0171] 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:

[0172]

[0173] The log-likelihood function is respectively applied to the state quantity x of the system at time k k , weight sequence γ 1:N and the weight sequence λ 1:N Taking the conditional expectation, we get:

[0174]

[0175] Step 3-2-2, first calculate the state x of the system at time k k The posterior probability and Q(Θ|Θ s ) formula, and use the particle filter algorithm to establish the posterior probability distribution of the system state quantity, and obtain:

[0176]

[0177] where δ(·) is the Dirac function, represents the lth particle at time k, represents the weight of the lth particle at time k;

[0178] According to all known observation data, the state quantity x of the system at time k is obtained by particle smoothing. k The posterior probability distribution of :

[0179]

[0180] in,

[0181] Get the state x of the system at time k and time k-1 k ,x k-1 The joint posterior probability distribution of is:

[0182]

[0183] in,

[0184] The first and second integral terms in the Q-function are:

[0185]

[0186] Where I1 is the sum of the first integral term and the second integral term in the Q-function, It is expressed as the particle weight after particle smoothing;

[0187] γ k and logγ k The expectation is expressed as <λ k > l and <logλ k > l , and there are:

[0188]

[0189]

[0190] in is the degree of freedom parameter of the measurement noise at the sth iteration, is the scale parameter of the measurement noise at the sth iteration,

[0191] The third and fourth integral terms in the Q-function are:

[0192]

[0193] Among them, I2 is the sum of the third integral term and the fourth integral term in the Q-function, It is represented as the joint weight of particles after particle smoothing; γ k and logγ k The expectation is expressed as <γ k > l and <logγ k > l , and there are:

[0194]

[0195]

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

[0197] Then the Q-function is expressed as:

[0198] Q(Θ|Θ s )=Q1(θ f ,R v )+Q2(R w )+Q3(r w )+Q4(r v )+C1

[0199] Where Q1(θ f ,R v ) represents the Q-function with respect to θ f ,R v Mathematical expression of Q2(R w ) represents the Q-function with respect to R w The mathematical expression of Q3(r w ) represents the Q-function with respect to r w Mathematical expression of Q4(r v ) represents the Q-function with respect to r v Mathematical expression of ; C1 represents the constant term, where:

[0200]

[0201]

[0202]

[0203]

[0204]

[0205] Where logp(u 1:N |Θ) is represented as the prior probability distribution of the input data, and logp(x1|γ1,Θ) is represented as the posterior probability distribution of the initial state.

[0206] Step 3-3, the step M-step includes:

[0207] Step 3-3-1, set Q1(θ f ,R w ) for R w Derivative, Q2(R v ) for R v Taking the derivative and setting it equal to zero, we get:

[0208]

[0209]

[0210] Step 3-3-2, Substituting into Q1(θ f ,R w )get:

[0211]

[0212] <λ k > l and <γ k > ij are the weights assigned to the square of the model prediction deviation corresponding to the state update process and the output measurement process data, respectively;

[0213] Using the damped Newton method to optimize the criterion function Q5(θ f ), and get the parameter θ f The optimized expression is:

[0214]

[0215] in, Denotes the parameter θ f The updated value for the s+1th time, μ represents the damping Newton constant, and I represents the identity matrix;

[0216] Step 3-3-3, Q3(r w ) w Derivative, Q4(r v ) v Taking the derivative and setting it equal to zero, we get:

[0217]

[0218]

[0219] Solve to get the parameter r w ,r v The estimates are

[0220] Step 3-3-4, through iterative cycle steps 3-3-1 to 3-3-3, obtain the model parameters of the system, and record the final iteration result as

[0221] The convergence conditions are:

[0222] Example

[0223] In this embodiment, the robust identification method for nonlinear state-space systems based on thick-tailed noise is specifically implemented according to the following steps:

[0224] The input of the algorithm is C obs = {u 1:N ,y 1:N}, the output of the algorithm includes: model parameters Scale parameter Degree of freedom parameter

[0225] Step 1: Construct observation data sets separately: C obs ;

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

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

[0228] Using the formula given in E-step, calculate the smoothed probability density function p(x k |y 1:N ,Θ s ),p(x k ,x k-1 |y 1:N ,Θ s ) and the expectations at each moment <γ k > ij , <logγ k > ij ,<λ k > l , <logλ k > l ;

[0229] Step 4: Update model parameter estimates:

[0230] Update model parameters θ f ,θ g , scale parameter R w ,R v , the degree of freedom parameter r w ,r v ;

[0231] Step 5: Set s=s+1 until the estimated parameters converge.

[0232] <1> Collect the input external force data of the double-link robotic arm system and the angular displacement data of the two nodes of the output section and perform preprocessing;

[0233] <2> Set the simulation parameters:

[0234] In order to verify that the method of the present invention is effective in solving problems such as abnormal values ​​in the state updating process and the output measurement process, it is assumed that the output data contains 10% abnormal values.

[0235] <3> Simulation verification:

[0236] For the convenience of description, the method of the present invention is abbreviated as: NSSM-Student. As a comparison, the same simulation model is simulated using the traditional method NSSM-Gaussian. The comparison of simulation results can well prove the effectiveness of the method of the present invention. In order to better verify the method of the present invention, the system model parameters are estimated for the NSSM-Student and NSSM-Gaussian methods when the proportion of abnormal values ​​in the state update process and the output measurement process is 10%. The specific simulation results are shown in Figure 2 . Figure 2 This is a comparison of the system parameter estimation curves of the NSSM-Student and NSSM-Gaussian methods when the proportion of abnormal values ​​in the state update process and the output measurement process is 10%. The solid curve represents the curve of the parameter estimation of the NSSM-Student method, the dotted curve represents the curve of the parameter estimation of the NSSM-Gaussian method, and the dotted straight line represents the true value of the system parameters. The closer it is to the dotted straight line, the better the parameter estimation effect.

[0237] Relative parameter estimation errors of NSSM-Student and NSSM-Gaussian methods when the proportion of outliers is 10% during the state update process and the output measurement process The root mean square error RMSE of the angular displacement of the two nodes is shown in Table 1:

[0238] Table 1

[0239]

[0240] Summary of simulation results:

[0241] When the proportion of outliers increases, the standard deviation of each parameter continues to decrease. The smaller the standard deviation, the better the estimation effect;

[0242] Under the condition of consistent outlier ratio, the parameter estimation effect of the NSSM-Student method of the present invention is significantly better than that of the NSSM-Gaussian method used as a comparison;

[0243] According to the simulation results in Table 1, it can be seen that compared with the method of the present invention, the RPEE value and RMSE value of the NSSM-Gaussian method are larger, indicating that the method of the present invention is highly robust to the presence of abnormal values ​​in the state update process and the output measurement process.

[0244] The present invention provides a method for robust identification of nonlinear state space systems based on thick-tailed noise. There are many methods and ways to implement the technical solution. The above is only a preferred implementation 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 nonlinear state-space systems based on thick-tailed noise, characterized in that: The following steps are involved: Step 1: Establish a continuous-time dynamic model of the two-link manipulator system, convert the continuous-time dynamic model into a discrete state-space model after finite difference, collect the vertical downward external force data of the end of the two-link manipulator system, and collect the angular displacement of the two links relative to the horizontal plane to form an identification data set; The discrete state space model includes a state transfer equation and a measurement equation; Step 2: For the discrete state space model, introduce the system noise distribution and the prior distribution of model parameters to form a probabilistic description of the dual-link manipulator system identification problem; Step 3: Based on the identification data set and the probabilistic description of the dual-link robot arm system identification problem, iteratively update under the framework of the expectation maximization algorithm to obtain the parameters to be identified; In step 1, the continuous-time dynamics model of the two-link manipulator system is expressed as: Among them, M, E are coefficient matrices; K(t) is a coefficient matrix, some elements of which change with time t; Φ is the angular displacement of the two joints, are the angular velocities of the two joints, are the two joint angular accelerations; and M, E, K(t) are specifically expressed as: Where m and l represent the mass and length of the two connecting rods respectively; g represents the acceleration of gravity; in order to simplify the system, the damping of the system is assumed to be concentrated at the nodes, d1 and d2 represent the damping coefficients of the two nodes; the connecting rod is assumed to be rigid, the elasticity of joint 1 and joint 2 is represented by two constant rotation springs, k1 and k2 represent the spring stiffness; α 10 ,α 20 They respectively represent the steady-state positions of the angular displacements of the two connecting rods relative to the horizontal plane; F(t) represents an external force applied to the end of the second connecting rod, which changes with time t and maintains an angle of α3 with the horizontal plane.

2. The method according to claim 1, characterized in that: In step 1, the state transfer equation is expressed as: x(k)=x(k-1)+Δ t Ax(k-1) Among them, set is the state quantity of the system, T represents the matrix transpose, x(k) and x(k-1) represent the state of the system at time k and time k-1 respectively; Δ t is the sampling time; A is the coefficient matrix, which is specifically expressed as: Among them, θ1, θ2 are the model parameters to be identified; b1, b2 are intermediate parameters.

3. The method according to claim 2, characterized in that In step 1, the measurement equation is expressed as: Wherein, y(k) represents the load end data obtained by measurement.

4. The method according to claim 3, characterized in that In step 1, F(t) is expressed as: F(t)=10sin(0.02πt).

5. The method according to claim 4, characterized in that In step 1, the discrete state space model is expressed as: x k =f(x k-1 ,u k-1 ,θ f )+w k y k =g(x k )+v k Among them, x k =x(k),x k-1 =x(k-1),x k 、x k-1 They represent the state of the two-link manipulator system at time k and time k-1, namely the angular velocity and angular displacement of the two nodes; u k-1 =F((k-1)Δt),u k-1 represents the external force input of the two-link manipulator system at time k-1, that is, the vertical downward external force at the end of the rod; y k =y(k),y k represents the output data collected by the dual-link robotic arm system, i.e., the angular displacement of the two nodes; f(x k-1 ,u k-1 ,θ f )=x k-1 +Δ t Ax k-1 is the nonlinear mapping function of the state update process, θ f =[θ1,θ2] is the parameter to be identified in the nonlinear mapping function; g(x k ) is the output measurement function, w k ,v k are the state noise and measurement noise of the two-link robotic arm system, respectively.

6. The method according to claim 5, characterized in that Step 2 includes: Step 2-1: Set w k ,v k Set to zero-mean Gaussian white noise, obeying the distribution N(w k |0,R w ),N(v k |0,R v );where R v is the scale parameter; The state transition local noise change parameter is R w / γ k ,get: w k |(R w ,c k )~N(w k |0,R w / c k ) Where R w is the scale parameter, N(w k |0,R w / γ k ) represents the noise sequence w k It obeys Gaussian distribution, with mean 0 and variance R w / γ k ; γ k is the weight parameter: γ k |r w ~Gamma(γ k |r w / 2,r w / 2) where r w represents the degree of freedom parameter, x k It obeys the following Gaussian distribution: x k |(f(x k-1 ,u k-1 ,i f ),R w ,c k )~N(w k |(f(x k-1 ,u k-1 ,i f ),R w / c k )); Step 2-2: x k The marginal distribution of k ,γ k |f(x k-1 ,u k-1 ,θ f ),R w ,r w ) for γ k Integrating gives: where Γ(·) represents the Gamma function, p(γ k |r w ) is the probability density function that follows the Gamma distribution; Step 2-3: Set the output measurement local noise change parameter to R v / λ k ,get: v k |(R v ,l k )~N(v k |0,R v / l k ) Where N(v k |0,R v / λ k ) represents the noise sequence v k It obeys Gaussian distribution, with mean 0 and variance R v / λ k ; k is the weight parameter: l k |r v ~Gamma(λ k |r v / 2,r v / 2) where v v represents the degree of freedom parameter, y k The Gaussian distribution is as follows: y k |(g(x k ),R v ,l k )~N(v k |g(x k ),R v / l k ); Step 2-4: k The marginal distribution of k ,λ k |g(x k ),R v ,r v ) for δ k Integrating gives:

7. The method according to claim 6, 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(Θ|Θ) of the two-link manipulator system s ): Where logp(C mis ,C obs |Θ) is the complete data set C={C mis ,C obs }'s log-likelihood function; Θ s is the parameter of the current nonlinear state space model; C obs = {u 1:N ,y 1:N } is the identification dataset of the two-link robotic arm system collected in step 1, u 1:N represents the external force input at the end of the link of the double-link robotic arm system, where u k =F(kΔt),Δ t is the sampling time, y 1:N represents the angular displacement of the two links of the double-link robotic arm system relative to the horizontal plane; C mis ={x 1:N ,γ 1:N ,λ 1:N } is a missing data set, that is, a data set that cannot be observed, where x 1:N represents the state quantity of the system, γ 1:N ,λ 1:N is the weight sequence, Θ={θ f ,R w ,r w ,R v ,r v } are the model parameters to be identified; represents logp(C mis ,C obs |Θ) conditional expectation; Step 3-3, execute step M-step, and search for a new parameter estimate Θ by optimizing the search method. s+1 , so that the parameter Θ s+1 The following constraints are met: 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 estimated parameters converge. The convergence condition is:

8. The method according to claim 7, characterized in that In step 3-2, the step E-step includes: 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: The log-likelihood function is respectively applied to the state quantity x of the system at time k k , weight sequence γ 1:N and the weight sequence λ 1:N Taking the conditional expectation, we get: Step 3-2-2, first calculate the state x of the system at time k k The posterior probability and Q(Θ|Θ s ) formula, and use the particle filter algorithm to establish the posterior probability distribution of the system state quantity, and obtain: where δ(·) is the Dirac function, represents the lth particle at time k, represents the weight of the lth particle at time k; According to all known observation data, the state quantity x of the system at time k is obtained by particle smoothing. k The posterior probability distribution of : in, Get the state x of the system at time k and time k-1 k ,x k-1 The joint posterior probability distribution of is: in, The sum of the first integral term and the second integral term in the Q-function is: Where I1 is the sum of the first integral term and the second integral term in the Q-function, It is expressed as the particle weight after particle smoothing; γ k and logγ k The expectation is expressed as <λ k > l and <logλ k > l , and there are: in is the degree of freedom parameter of the measurement noise at the sth iteration, is the scale parameter of the measurement noise at the sth iteration; The sum of the third and fourth integral terms in the Q-function is: Among them, I2 is the sum of the third integral term and the fourth integral term in the Q-function, It is represented as the joint weight of particles after particle smoothing; γ k and logγ k The expectation is expressed as <γ k > l and <logγ k > l , and there are: Where, ψ(x)=dlog(Γ(x)) / dx is the Digamma function; Then the Q-function is expressed as: Q(Θ|Θ s )=Q1(θ f ,R v )+Q2(R w )+Q3(r w )+Q4(r v )+C1 Where Q1(θ f ,R v ) represents the Q-function with respect to θ f ,R v Mathematical expression of Q2(R w ) represents the Q-function with respect to R w The mathematical expression of Q3(r w ) represents the Q-function with respect to r w Mathematical expression of Q4(r v ) represents the Q-function with respect to r v Mathematical expression of ; C1 represents the constant term, where: Where logp(u 1:N |Θ) is represented as the prior probability distribution of the input data, and logp(x1|γ1,Θ) is represented as the posterior probability distribution of the initial state.

9. The method according to claim 8, characterized in that Step 3-3, the step M-step includes: Step 3-3-1, set Q1(θ f ,R w ) for R w Derivative, Q2(R v ) for R v Taking the derivative and setting it equal to zero, we get: Step 3-3-2, Substituting into Q1(θ f ,R w )get: <λ k > l and <γ k > ij are the weights assigned to the square of the model prediction deviation corresponding to the state update process and the output measurement process data, respectively; Using the damped Newton method to optimize the criterion function Q5(θ f ), and get the parameter θ f The optimized expression is: in, Denotes the parameter θ f The updated value for the s+1th time, μ represents the damping Newton constant, and I represents the identity matrix; Step 3-3-3, Q3(r w ) w Derivative, Q4(r v ) v Taking the derivative and setting it equal to zero, we get: Solve to get the parameter r w ,r v The estimates are Step 3-3-4, through iterative cycle steps 3-3-1 to 3-3-3, obtain the model parameters of the system, and record the final iteration result as The convergence conditions are:

Images