Robot load dynamics parameter identification method based on current dynamics

Abstract

The invention discloses a robot load dynamics parameter identification method based on current dynamics, and relates to the field of robot dynamics parameter calibration. The problem that an existing load kinetic parameter identification method has errors in joint torque constant acquisition, joint torque prediction errors are increased, and load kinetic parameter identification errors are accumulated is solved. According to the method, a new model for load kinetic parameter identification is constructed; the identification method based on weighted iteration is used for estimating kinetic parameters of a current level, and is beneficial for obtaining a weighted least square solution and eliminating abnormal values. In addition, in the kinetic parameter identification process of the current level, a zero-speed continuous nonlinear joint combination friction model is introduced to improve the overall identification precision. The method is mainly used for robot current level kinetic parameter identification under the condition that the joint torque is not used or is difficult to obtain.

Description

Technical Field

[0001] The invention relates to the field of robot dynamics parameter calibration. Background Art

[0002] As the load capacity of robots increases, the proportion of load-induced current or torque in actuators increases accordingly, increasing their importance in robotics. For example, accurate load dynamics parameters can help accurately predict joint current or torque, enhance robot motion control performance, and improve human-robot safe interaction capabilities.

[0003] When the load is fixedly connected to the end of the robot, there are two main identification methods. One is to use force or torque sensors installed on the robot joints to identify the dynamic characteristics of the load. However, the high cost and large size and weight make this method not feasible in some cases. The load is fixedly connected to the end effector, which can be regarded as part of the robot's end link. Therefore, another method is to identify the load dynamic parameters by identifying the dynamic parameters of the robot itself. This method is widely used, specifically:

[0004] Swevers et al. provided a classic load identification model in the literature. For loads connected to the robot joint torque and links, this method distinguishes the load from the end link, thereby reducing the identification error. Khalil et al. proposed two load identification methods. One is to let the robot track the same reference trajectory twice: first without load, and then with load. The difference between the two joint torques is attributed to the load; the other method uses weighted least squares to simultaneously identify the robot's dynamic parameters and load parameters, but does not consider friction changes. Gaz et al. proposed a static identification method based on the linear correlation between the changes in the robot's static parameters and the changes in the load's static parameters. Bahloul et al. tried to use nonlinear optimization techniques to identify the dynamic parameters of the load, applying basic physical constraints such as ensuring that the load mass is positive. Xu et al. proposed an identification method based on double weighting technology to improve the estimation accuracy of the payload dynamics.

[0005] The above algorithms provide a good theoretical basis and valuable reference for load dynamic parameter identification. However, these methods all rely on joint torque signals to identify load dynamic parameters. In fact, most industrial robots are not equipped with joint torque sensors, and joint torque can only be predicted by multiplying the actuator current with the joint torque constant, but the joint torque constant provided by the manufacturer usually has an uncertainty of about 10%. Even if some precise identification methods are introduced to obtain these constants, they still have errors. Therefore, the increase in joint torque prediction error will cause the load dynamic parameter identification error to accumulate. Summary of the invention

[0006] The purpose of the present invention is to solve the problem that the existing load dynamics parameter identification method has errors in obtaining the joint torque constant, which increases the joint torque prediction error and leads to the accumulation of load dynamics parameter identification errors. The present invention provides a robot load dynamics parameter identification method based on current dynamics.

[0007] A method for identifying robot load dynamics parameters based on current dynamics, the method comprising:

[0008] S1. Construct the joint combination friction model τ of the robot f , for τ f Optimize and solve to determine the robot's load-carrying state in the process of executing the excitation trajectory at each sampling time. f The corresponding combination coefficients α, V s and K v The value of α and V at the sampling time s and K v Construct the corresponding linear regression matrix H f ; where α is the joint velocity index vector, V s is the Stribeck velocity vector and K v is the inverse tangent coefficient vector;

[0009] When the robot executes the excitation trajectory under load, the linear regression matrix H at the sampling moment is constructed according to a set of joint coordinates corresponding to each sampling moment. l ; Each set of joint coordinates includes the angle position q, velocity and acceleration

[0010] S2, according to H at each sampling time l and H f , construct the matrix H e =[H l H f ];

[0011] Collect all joint current data corresponding to the two states at each sampling time when the robot executes the excitation trajectory in the loaded and unloaded states, and merge the data in the two states to form a set of joint current data corresponding to the current sampling time i e ;

[0012] S3, H at N sampling moments e Stack in chronological order to form a matrix i at N sampling moments e Stack in chronological order to form a matrix Where n is the total number of robot joints, c is H lThe number of elements in Re Nn×(c+3n) is a real number space with dimension Nn×(c+3n), Re Nn×1 is a real number space with dimension Nn×1;

[0013] S4. According to and Solve the constructed kinetic parameter identification model [χ l χ f ] T The OLS solution of l,0 χ f,0 ] T ∈Re (c+3n)×1 , according to [χ l,0 χ f,0 ] T Determine the initial value of the covariance matrix σ; where χ l is the load dynamics parameter at the current level, χ f is the friction variation parameter caused by load at the current level, χ l,0 is l The OLS solution of f,0 is f OLS solution of ;

[0014] S5, initialize the joint data weight vector P = ones (Nn, 1) and the stacked joint data weight vector P e =P·ones(1,c+3n);P∈Re Nn×1 , P e ∈Re Nn×(c+3n) , ones(Nn,1) is a column vector of N×n rows, and all elements of the column vector are 1; ones(1,c+3n) is a row vector of c+3n columns, and all elements of the row vector are 1;

[0015] S6, use the covariance matrix σ to and Weighted, the normalized motor current vector is obtained and the normalized linear regression matrix S7, according to P and Calculating the Matrix According to P e and Calculating the Matrix

[0016] S8. According to and Calculate [χ l χ f ] T ;

[0017] S9. According to and [χ l χ f ] T , calculate the normalized current residual vector R # ∈Re Nn×1 ;

[0018] S10, according to R # Update σ, according to the threshold δ and R # Update P, update P by the updated P e ;

[0019] S11, determine whether P converges, if yes, output the x in step S8 l , the result is no, return to step S6.

[0020] Preferably, in step S1, the joint combination friction model τ of the robot is constructed f The implementation is:

[0021] First, design the jth joint friction model τ f,j ,and

[0022]

[0023] Among them, F s,j 、F c,j 、F v,j and F b,j They represent the static friction coefficient, Coulomb friction coefficient, viscous friction coefficient and friction offset coefficient of the jth joint respectively. is the velocity of the jth joint, α j is the velocity index of the jth joint, K v,j is the inverse tangent coefficient of the jth joint, V s,j is the Stribeck velocity of the jth joint, and e is a natural constant;

[0024] Secondly, according to all joint friction models, the joint combination friction model τ is constructed f , where τ f =[τ f,1 ,τ f,2 ,…,τ f,n ].

[0025] Preferably, in step S1, τ is determined at each sampling time. f The corresponding combination coefficients α, V s and K v The implementation method of the value is:

[0026] S1-1, collect a set of joint current data i' corresponding to the current sampling time when the robot is executing the excitation trajectory in the loaded statee , and a set of joint coordinates corresponding to the current sampling time; where each set of joint current data i′ e Includes all joint current data under load at the current sampling moment

[0027] S1-2. According to a set of joint coordinates corresponding to the current sampling moment, construct the linear regression matrix H at the sampling moment l ;

[0028] S1-3, according to i′ e , H l and χ′, estimate the friction current i at the current sampling time f =i′ e -H l ·[χ′,0] T ;

[0029] χ′ is given by [χ l χ f ] T The current inertia parameter vector formed by the mass of each link of the robot, the first-order moment of each link, the Coriolis force of each link, the centrifugal force of each link, the gravity of each link, and the moment of inertia of the motor rotor of each joint;

[0030] S1-4, for τ f Constrain and optimize α, V s and K v , specifically:

[0031]

[0032] in, is the optimized variable set of the jth joint,

[0033] α=[α1,α2,……,α n ], α j is the velocity index of the jth joint;

[0034] K v =[K v,1 ,K v,2 ,……,K v,j ], K v,j is the inverse tangent coefficient of the jth joint;

[0035] V s =[V s,1 ,V s,2 ,……,V s,n ],V s,j is the Stribeck velocity of the jth joint;

[0036] K∈Re n×nis a diagonal matrix of constant coefficients of joint torques, and τ = Ki′ e , τ is the vector composed of all joint torques collected.

[0037] Preferably, in step S4, the expression of the kinetic parameter identification model constructed is:

[0038] Preferably, in step S4, solving [χ l χ f ] T The OLS solution of l,0 χ f,0 ] T The implementation is:

[0039] Preferably, in step S4, according to [χ l,0 χ f,0 ] T The implementation method for determining the initial value of the covariance matrix σ is:

[0040] According to [ l,0 χ f,0 ] T Determine the initial value of the current residual vector R0∈Re Nn×1 ,and Re Nn×1 is a real number space with dimension Nn×1;

[0041] According to R0, the initial value of the covariance matrix σ is obtained

[0042] Preferably, in step S6,

[0043] Preferably, in step S8,

[0044] Preferably, in step S9, R # ∈Re Nn×1 .

[0045] Preferably, in step S10, according to R # The implementation of updating σ is:

[0046]

[0047] According to the threshold δ and R # The implementation of updating P is:

[0048] P=P j ⊙Φ(R # ,δ);

[0049] Among them, Φ(R # ,δ) is a function of a vector of size Nn×1, and R # Elements in that exceed the threshold δ will be set to 0, otherwise they will be set to 1;

[0050] Update P by the updated P e The implementation is:

[0051] P e =P·ones(1,c+3n).

[0052] Advantages of the present invention:

[0053] The present invention proposes a method for identifying robot load dynamics parameters based on current dynamics. The present invention uses current instead of joint torque to identify load dynamics parameters. The method of the present invention is better than the comparative method in predicting effective load current dynamics parameters. It has a smaller root mean square error (RMSE) for both concentric loads and eccentric loads.

[0054] On the one hand, the method proposed in the present invention is derived from current dynamics, so there is no need for joint torque constants, which avoids the joint torque estimation error caused by inaccurate identification of joint torque constants and reduces the cumulative identification error. Therefore, unlike the method based on torque dynamics, the effective load identification method based on current level dynamics proposed in the present invention performs better in the estimation accuracy of effective load dynamics.

[0055] On the other hand, the proposed iterative weighted estimation helps to mitigate the negative impact of outliers in the measurement results, which also improves the identification accuracy of the effective load current dynamics. In addition, the use of a continuous nonlinear joint combined friction model at the velocity reversal point can effectively avoid the problem of friction force mutation when the joint velocity is reversed, which helps to improve the overall identification accuracy of the robot current dynamics. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 is a flow chart of the method for identifying robot load dynamics parameters based on current dynamics according to the present invention;

[0057] Figure 2 This is a schematic diagram of the end loading of the UR 10 robot; Figure 2 a is a schematic diagram of the structure when the robot end is not loaded. Figure 2 b is the structural diagram of the robot end when a concentric load is loaded. Figure 2 c is the structural diagram of the robot when an eccentric load is loaded at the end;

[0058] Figure 3 is the current dynamics reconstruction of the concentric load; where, Figure 3a to Figure 3 f are the comparison diagrams of the predicted values ​​and measured values ​​of the 1st to 6th joint currents respectively; the joint current is the current driving the joint motor, and the joint current includes the friction current (i.e. the current to overcome friction), the connecting rod inertia current (the current to overcome the connecting rod inertia), the joint inertia current (the current to overcome the joint inertia without loading and without considering friction), and the load current (the current to overcome the load resistance to the joint)

[0059] Figure 4 The schematic diagram of the distribution of the normalized current residuals of each joint obtained by the method of the present invention when δ=2.5;

[0060] Figure 5 The schematic diagram of the distribution of the normalized current residuals of each joint obtained by the method of the present invention when δ=3. DETAILED DESCRIPTION

[0061] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0062] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments may be combined with each other.

[0063] Specific implementation method 1. Combination Figure 1 The present embodiment is described. The method for identifying robot load dynamics parameters based on current dynamics described in the present embodiment includes:

[0064] S1. Construct the joint combination friction model τ of the robot f , for τ f Optimize and solve to determine the robot's load-carrying state in the process of executing the excitation trajectory at each sampling time. f The corresponding combination coefficients α, V s and K v The value of α and V at the sampling time s and K v Construct the corresponding linear regression matrix H f ; where α is the joint velocity index vector, V s is the Stribeck velocity vector and K v is the inverse tangent coefficient vector;

[0065] When the robot executes the excitation trajectory under load, the linear regression matrix H at the sampling moment is constructed according to a set of joint coordinates corresponding to each sampling moment.l ; Each set of joint coordinates includes the angle position q, velocity and acceleration

[0066] S2, according to H at each sampling time l and H f , construct the matrix H e =[H l H f ];

[0067] Collect all joint current data corresponding to the two states at each sampling time when the robot executes the excitation trajectory in the loaded and unloaded states, and merge the data in the two states to form a set of joint current data corresponding to the current sampling time i e ;

[0068] S3, H at N sampling moments e Stack in chronological order to form a matrix i at N sampling moments e Stack in chronological order to form a matrix Where n is the total number of robot joints, c is H l The number of elements in Re Nn×(c+3n) is a real number space with dimension Nn×(c+3n), Re Nn×1 is a real number space with dimension Nn×1;

[0069] S4. According to and Solve the constructed kinetic parameter identification model [χ l χ f ] T The OLS solution of l,0 χ f,0 ] T ∈Re (c+3n)×1 , according to [χ l,0 χ f,0 ] T Determine the initial value of the covariance matrix σ; where χ l is the load dynamics parameter at the current level, χ f is the friction variation parameter caused by load at the current level, χ l,0 is l The OLS solution of f,0 is f OLS solution of ;

[0070] S5, initialize the joint data weight vector P = ones (Nn, 1) and the stacked joint data weight vector P e =P·ones(1,c+3n);P∈ReNn×1 , P e ∈Re Nn×(c+3n) , ones(Nn,1) is a column vector of N×n rows, and all elements of the column vector are 1; ones(1,c+3n) is a row vector of c+3n columns, and all elements of the row vector are 1;

[0071] S6, use the covariance matrix σ to and Weighted, the normalized motor current vector is obtained and the normalized linear regression matrix

[0072] As an example,

[0073] S7, according to P and Calculating the Matrix According to P e and Calculating the Matrix

[0074] S8. According to and Calculate [χ l χ f ] T ; As an example,

[0075] S9. According to and [χ l χ f ] T , calculate the normalized current residual vector R # ∈Re Nn×1 ; As an example,

[0076] S10, according to R # Update σ, according to the threshold δ and R # Update P, update P by the updated P e ;

[0077] S11, determine whether P converges, if yes, output the x in step S8 l , the result is no, return to step S6.

[0078] In order to mitigate the impact of measurement noise or model uncertainty on outliers in the recognition results, a data weighting vector P and a matrix P based on P expansion are used. e , to calculate and in, In specific applications, if P and P e If an element in is 0, it means and If the corresponding element in is an outlier, the value will be set to 0, otherwise it is a non-outlier and will remain unchanged. After removing the outlier, and will be converted to and The operator ⊙ represents the multiplication of corresponding elements.

[0079] In specific applications, constructing a linear regression matrix can be achieved through existing technologies.

[0080] [ l χ f ] T is the robot current level dynamics parameter, including the current level load dynamics parameter χ l and the friction change parameter χ caused by the load at the current level f ; In specific applications, the load dynamics parameter χ l Specifically, it may include the mass of each connecting rod, the first-order moment of each connecting rod, the centrifugal force of each connecting rod, the Coriolis force of each connecting rod, the gravity of each connecting rod, the moment of inertia of each joint, the static friction parameter of each joint, the Coulomb friction parameter of each joint, the viscous friction parameter of each joint, the friction offset coefficient information of each joint, and the friction change parameter χ f Specifically, it may include F c,j 、F v,j and F b,j (j=1...n).

[0081] In this embodiment, the current is used instead of the joint torque to identify the load dynamic parameters. The method of the present invention is superior to the comparative method in predicting the effective load current dynamic parameters. It has a smaller root mean square error (RMSE) for both concentric loads and eccentric loads. Specifically, a new model for load dynamic parameter identification is derived based on current dynamics; a weighted iteration-based identification method is proposed to estimate the dynamic parameters at the current level, which helps to obtain a weighted least squares solution and eliminate outliers. In addition, in the process of dynamic parameter identification at the current level, a zero-speed continuous nonlinear joint combined friction model is introduced to improve the overall identification accuracy. The robot current level dynamic parameters [χ l χ f ] T Then, using this [χ l χ f ] T and the dynamic basis parameter vector χ of the robot current level under no-load conditions a , predict the estimated value of joint current, specifically through formula 9.

[0082] The present invention introduces a zero-speed continuous nonlinear joint combination friction model to improve the overall identification accuracy, and gives the joint combination friction model τ of the robot constructed in step S1. f The implementation is:

[0083] First, for the jth joint, the Stribeck friction model with velocity exponential term is used to characterize the viscous friction phenomenon, and the inverse tangent function is used instead of the sign function to improve the discontinuity of the friction model at the velocity reversal point. For the joint j, a joint friction model is proposed and used as follows:

[0084] Design the j-th joint friction model τ f,j ,and

[0085]

[0086] Among them, F s,j 、F c,j 、F v,j and F b,j They represent the static friction coefficient, Coulomb friction coefficient, viscous friction coefficient and friction offset coefficient of the jth joint respectively. is the velocity of the jth joint, α j is the velocity index of the jth joint, K v,j is the inverse tangent coefficient of the jth joint, V s,j is the Stribeck velocity of the jth joint, and e is a natural constant;

[0087] Secondly, according to all joint friction models, the joint combination friction model τ is constructed f , where τ f =[τ f,1 ,τ f,2 ,…,τ f,n ].

[0088] Specifically, in step S1, τ is determined at each sampling time. f The corresponding combination coefficients α, V s and K v The implementation method of the value is:

[0089] S1-1, collect a set of joint current data i' corresponding to the current sampling time when the robot is executing the excitation trajectory in the loaded state e , and a set of joint coordinates corresponding to the current sampling time; where each set of joint current data i′ e Includes all joint current data under load at the current sampling moment

[0090] S1-2. According to a set of joint coordinates corresponding to the current sampling moment, construct the linear regression matrix H at the sampling moment l;

[0091] S1-3, according to i′ e , H l and χ′, estimate the friction current i at the current sampling time f =i′ e -H l ·[χ′,0] T ;

[0092] χ′ is given by [χ l χ f ] T The current inertia parameter vector formed by the mass of each link of the robot, the first-order moment of each link, the Coriolis force of each link, the centrifugal force of each link, the gravity of each link, and the moment of inertia of the motor rotor of each joint, [χ l χ f ] T are the current level dynamic parameters including load dynamics and friction changes;

[0093] S1-4, for τ f Constrain and optimize α, V s and K v , specifically:

[0094]

[0095] in, is the optimized variable set of the jth joint,

[0096] α=[α1,α2,……,α n ], α j is the velocity index of the jth joint;

[0097] K v =[K v,1 ,K v,2 ,……,K v,j ], K v,j is the inverse tangent coefficient of the jth joint;

[0098] V s =[V s,1 ,V s,2 ,……,V s,n ],V s,j is the Stribeck velocity of the jth joint;

[0099] K∈Re n×n is a diagonal matrix of constant coefficients of joint torques, and τ = Ki′ e , τ is the vector composed of all joint torques collected.

[0100] In this preferred embodiment, the method for determining α, V is given. s and K v The implementation method of the value of can better handle the nonlinear parameter space and obtain a faster convergence speed by solving the optimization equation.

[0101] In step S4, the expression of the constructed kinetic parameter identification model is: The derivation process is:

[0102] For a serial robot with n rigid links, considering the dynamic characteristics of the links and the motor rotor, the dynamic differential equation can be described in the following concise form:

[0103]

[0104] In the formula, q, are the coordinates of the robot, respectively representing the robot joint angle position, velocity and acceleration, M(q) is the inertia matrix, I a is the equivalent inertia between the motor rotor and the transmission system, is the matrix composed of the centrifugal force and Coriolis force of the connecting rod, g(q) is the gravity vector of the connecting rod, and for the jth connecting rod, its mass, centrifugal force, Coriolis force, and gravity are M respectively. j , C j and g j . is the joint friction torque vector, and τ is the joint output torque vector.

[0105] For the jth joint, the Stribeck friction model with velocity exponential term is used to characterize the viscous friction phenomenon, and the inverse tangent function is used instead of the sign function to improve the discontinuity of the friction model at the velocity reversal point. For the joint j, a friction model is proposed and used as follows:

[0106]

[0107] In the formula, F s,j 、F c,j 、F v,j and F b,j Respectively represent the static friction coefficient, Coulomb friction coefficient, viscous friction coefficient and friction offset coefficient of the jth joint, all of which are linear friction coefficients. f,j is the friction torque of the jth joint, α j is the velocity index of the jth joint, V s,j is the Stribeck velocity of the jth joint, K v,j is the inverse tangent coefficient of the jth joint. j 、V s,j , K v,jCollectively referred to as the nonlinear friction coefficient of the j-th joint.

[0108] Formula (2) is a scalar equation. For a robot with n degrees of freedom, the friction scalar equation of each joint can be combined into a matrix form (j = 1, ..., n), that is, τ f =[τ f,1 ,τ f,2 ,…,τ f,n ], and then substitute it into formula (1), and linearize the formula obtained after substitution, that is, formula (1) can be written into the following linear form:

[0109]

[0110] Where ψ∈Re r×1 is the basis parameter set, L∈Re n×r is the robot dynamics torque level regression matrix, is q, α,K v ,V s Here, α, K v 、V s They are α j 、V s,j , K v,j The three parameters are collectively referred to as the combined friction coefficient. r represents the number of columns in the L matrix, that is, the number of basis parameters contained in ψ. For simplicity, each matrix will appear in the form of its first letter, for example, L will replace

[0111] Since the motor current i∈Re n×1 There is usually a linear relationship between and the joint torque τ:

[0112] τ=Ki (4);

[0113] Where K∈Re n×n is a diagonal matrix of constant coefficients of joint torques. For a certain model of robot, it can be obtained by measuring the joint torques τ and motor currents i during the operation of the robot and calculating them using formula (4). By substituting formula (4) into formula (1), and then dividing all the terms on the right side of the equal sign of formula (1) by K, an expression for the robot motor current i will be obtained. This expression contains the information of all joint currents, where the joint current expression for the jth joint is:

[0114]

[0115] Where M j ,I a,j , C j They are M and I respectively.a and the j-th row element of C, g j is the jth element of g, K j is the j-th diagonal element of K.

[0116] From formula (3), we can see that formula (5) can also perform similar linearization operations:

[0117] i j =L j ψ j (6);

[0118] In the formula, i j is the current of the motor of the jth joint, L j It can be obtained by extracting the j-th row element of L in formula (3),

[0119] For a robot with n joint degrees of freedom, these n scalar equations can be expressed as:

[0120] i1=L1ψ 1 ,i2=L2ψ 2 ,…,i n =L n ψ n ;

[0121] Then, by stacking the above n scalar equations, the matrix form of equation (6) can be obtained:

[0122]

[0123] If the linearly correlated columns in all matrices of formula (7) are removed, formula (7) can be further simplified:

[0124] i = Hχ (8);

[0125] Where χ∈Re r×1 is the basic parameter set at the current level, H∈Re n×r is the regression matrix at the current level. Similar to L, H is also a full-rank matrix. r represents the number of elements in. Similar to L, H is related to q, α, K v and V s There is a similar functional relationship between them.

[0126] If the robot executes the same excitation trajectory twice with and without a load at the end, the following set of dynamic equations can be obtained:

[0127]

[0128] Wherein, the subscripts a and b represent the two cases of carrying no load and carrying load respectively. In the two cases of executing the excitation trajectory, the measured robot joint currents are i and a ∈Re n×1 and i b ∈Re n×1 .H a ∈Re n×r and H b ∈Re n×r are the regression matrices calculated based on the joint coordinates of the excitation trajectory in the two cases. a ∈Re r×1 is the dynamic basis parameter vector of the robot current level under no-load conditions, including the robot inertia basis parameter vector χ a,i and the friction basis parameter vector χ when there is no load a,f , and χ a =[χ a,i , χ a,f ] T . χ l ∈Re c×1 is the load dynamics parameter at the current level, which includes the mass of each connecting rod, the first-order moment of each connecting rod, the centrifugal force of each connecting rod, the Coriolis force of each connecting rod, the gravity of each connecting rod, the moment of inertia of each joint, the static friction parameter of each joint, the Coulomb friction parameter of each joint, the viscous friction parameter of each joint, and the friction offset coefficient of each joint. l ∈Re n×c is its corresponding regression matrix. c is H l The number of elements in . f ∈Re 3n×1 is the friction change parameter caused by the load at the current level. Each joint has three dynamic friction parameters. For the jth joint, these three parameters are the F c,j 、F v,j and F b,j Divide by K j The result, H f ∈Re n×3n is its corresponding regression matrix.

[0129] If the robot executes the same excitation trajectory with and without a load, then theoretically H a ≈H b Combining it with formula (9) and sorting it out, we can get the following formula:

[0130] i b =i a +H l χ l +H f χ f (10);

[0131] The above formula can be expressed as the multiplication of two matrices:

[0132] i e =H e [ l χ f ] T (11);

[0133] In the formula, i e ∈Re n×1 Equal to i b -i a , H e ∈Re n×(c+3n) Equal to [H l H f ]. Among them, H e Still about q, α, K v and V s The function of H l It is about q, The function of H f It is about α, K v and V s function.

[0134] Formula (11) is a kinetic parameter identification model at one sampling time, and for practical applications, the present invention uses data at multiple sampling times to identify [χ l χ f ] T , so the formula (11) is replaced and transformed into:

[0135]

[0136] Formula (12) is the dynamic parameter identification model at the current level established by the present invention. For multiple moments e The stacking of H for multiple moments e of stacking.

[0137] When the robot executes the excitation trajectory with and without a load, after collecting N sets of sampling data (joint current, joint coordinates), the collected motor current vector and the calculated regression matrix It can be listed in the following ways:

[0138]

[0139] In the formula, i e,1 to i e,N Represent the motor current vectors measured from 1 to N, respectively, He,1 To H e,N Represent the regression matrices of the 1st to Nth measurements respectively.

[0140] In specific application, in step S4, solving [χ l χ f ] T The OLS solution of l,0 χ f,0 ] T The implementation is: By using The generalized inverse matrix of Solve [χ l,0 χ f,0 ] T , avoiding Irreversible problem, and can obtain [χ l,0 χ f,0 ] T The least squares solution of .

[0141] In step S4, according to [ l,0 χ f,0 ] T The implementation method for determining the initial value of the covariance matrix σ is:

[0142] According to [ l,0 χ f,0 ] T Determine the initial value of the current residual vector R0∈Re Nn×1 ,and Re Nn×1 is a real number space with dimension Nn×1;

[0143] According to R0, the initial value of the covariance matrix σ is obtained

[0144] In step S10, according to R # The implementation of updating σ is:

[0145]

[0146] According to the threshold δ and R # The implementation of updating P is:

[0147] P=P j ⊙Φ(R # ,δ);

[0148] Among them, Φ(R # ,δ) is a function of a vector of size Nn×1, and R # Elements exceeding the threshold δ will be set to 0, otherwise they will be set to 1. In specific applications, the threshold δ is set to 2.5;

[0149] Update P by the updated P e The implementation is:

[0150] P e =P·ones(1,c+3n).

[0151] Verification test:

[0152] The UR 10 robot is used to verify the method proposed in the present invention. The UR 10 robot has 6 joints. The end of the UR 10 robot is unloaded, loaded concentrically, and loaded eccentrically as shown in the following examples. Figure 2 (a) Figure 2 (b) Figure 2 (c) as shown.

[0153] The method proposed in this invention is used to identify the current plane load dynamics parameters χ when the robot end is loaded with concentric load and eccentric load. l and the friction change parameter χ caused by the load at the current level f The results are shown in Table 1 and Table 2 respectively.

[0154] Table 1 Load dynamic parameters and friction variation parameters when the robot end is loaded with concentric load

[0155]

[0156]

[0157]

[0158] Table 2 Load dynamic parameters and friction variation parameters when the robot end is loaded with eccentric load

[0159]

[0160]

[0161] From Table 1 and Table 2 above, l,1 To l,56 represents the load current dynamic parameters of the robot’s six joints, χ f,1 To l,18 It represents the friction variation parameters of the current level caused by the load of the six joints of the robot. From the data given in Table 1 and Table 2, it can be seen that the standard deviation of all parameter identification results is controlled within 0.01, so the identification results are relatively reliable.

[0162] In order to compare the dynamic parameter identification method proposed in the present invention with the methods used in the prior art (named as Method 1 (using joint torque to identify robot dynamic parameters), Method 2 (using weighted least squares method to simultaneously identify robot dynamic parameters and load parameters) and Method 3 (identifying robot dynamic parameters based on double weighting technology)), the RMSE of the joint current measured by the four methods was compared under the condition of concentric and eccentric loads at the end, as shown in Table 3 and Table 4, respectively.

[0163] Table 3 Current residual RMSE of concentric load identification using four different methods

[0164]

[0165] Table 4 Current residual RMSE of eccentric load identified by four different methods

[0166]

[0167] It can be found from Table 3 and Table 4 that, for each joint, the current residual identified by the method of the present invention is the smallest compared with other methods, so the method of the present invention has the highest accuracy compared with the other three methods.

[0168] Figure 3 The current in the middle vertical axis is the robot joint current i under load b .from Figure 3 It can be seen that, without considering the measurement noise, the predicted value and the measured value can fit well, so the method proposed in the present invention can more accurately predict the current of each joint when the robot carries a concentric load.

[0169] Figure 4 and Figure 5 The distribution of normalized current residuals of each joint when the threshold δ is given is based on Figure 4 and Figure 5 It can be seen that when the threshold δ is set to 2.5, the normal probability plot does not show "heavy tails", but rather presents a nearly straight line with a unit slope within the threshold range, indicating that the threshold setting is appropriate. On the contrary, when the value is set to 3.5, the normal probability plot shows "heavy tails", indicating that the normalized residual does not fully follow the normal distribution and the threshold setting is unreasonable.

[0170] Although the present invention is described herein with reference to specific embodiments, it should be understood that these embodiments are merely examples of the principles and applications of the present invention. It should therefore be understood that many modifications may be made to the exemplary embodiments and that other arrangements may be devised without departing from the spirit and scope of the present invention as defined by the appended claims. It should be understood that the various dependent claims and features described herein may be combined in a manner different from that described in the original claims. It should also be understood that features described in conjunction with individual embodiments may be used in other described embodiments.

Claims

1. A method for identifying robot load dynamics parameters based on current dynamics, characterized in that: The method includes: S1. Construct the joint combination friction model τ of the robot f , for τ f Optimize and solve to determine the robot's load-carrying state in the process of executing the excitation trajectory at each sampling time. f The corresponding combination coefficients α, V s and K v The value of α and V at the sampling time s and K v Construct the corresponding linear regression matrix H f ; where α is the joint velocity index vector, V s is the Stribeck velocity vector and K v is the inverse tangent coefficient vector; When the robot executes the excitation trajectory under load, the linear regression matrix H at the sampling moment is constructed according to a set of joint coordinates corresponding to each sampling moment. l ; Each set of joint coordinates includes the angle position q, velocity and acceleration S2, according to H at each sampling time l and H f , construct the matrix H e =[H l H f ]; Collect all joint current data corresponding to the two states at each sampling time when the robot executes the excitation trajectory in the loaded and unloaded states, and merge the data in the two states to form a set of joint current data corresponding to the current sampling time i e ; S3, H at N sampling moments e Stack in chronological order to form a matrix i at N sampling moments e Stack in chronological order to form a matrix Where n is the total number of robot joints, c is H l The number of elements in Re Nn×(c+3n) is a real number space with dimension Nn×(c+3n), Re Nn×1 is a real number space with dimension Nn×1; S4. According to and Solve the constructed kinetic parameter identification model [χ l χ f ] T The OLS solution of l,0 χ f,0 ] T ∈Re (c+3n)×1 , according to [χ l,0 χ f,0 ] T Determine the initial value of the covariance matrix σ; where χ l is the load dynamics parameter at the current level, χ f is the friction variation parameter caused by load at the current level, χ l,0 is l The OLS solution of f,0 is f OLS solution of ; S5, initialize the joint data weight vector P = ones (Nn, 1) and the stacked joint data weight vector P e =P·ones(1,c+3n);P∈Re Nn×1 , P e ∈Re Nn×(c+3n) , ones(Nn,1) is a column vector of N×n rows, and all elements of the column vector are 1; ones(1,c+3n) is a row vector of c+3n columns, and all elements of the row vector are 1; S6, use the covariance matrix σ to and Weighted, the normalized motor current vector is obtained and the normalized linear regression matrix S7, according to P and Calculating the Matrix According to P e and Calculating the Matrix S8. According to and Calculate [χ l χ f ] T ; S9. According to and [χ l χ f ] T , calculate the normalized current residual vector R # ∈Re Nn×1 ; S10, according to R # Update σ, according to the threshold δ and R # Update P, update P by the updated P e ; S11, determine whether P converges, if yes, output the x in step S8 l , the result is no, return to step S6.

2. The method for identifying robot load dynamics parameters based on current dynamics according to claim 1, characterized in that: In step S1, the joint combination friction model τ of the robot is constructed f The implementation is: First, design the jth joint friction model τ f,j ,and Among them, F s,j 、F c,j 、F v,j and F b,j They represent the static friction coefficient, Coulomb friction coefficient, viscous friction coefficient and friction offset coefficient of the jth joint respectively. is the velocity of the jth joint, α j is the velocity index of the jth joint, K v,j is the inverse tangent coefficient of the jth joint, V s,j is the Stribeck velocity of the jth joint, and e is a natural constant; Secondly, according to all joint friction models, the joint combined friction model τ is constructed f , where τ f =[τ f,1 ,τ f,2 ,…,τ f,n ].

3. The method for identifying robot load dynamics parameters based on current dynamics according to claim 2, characterized in that: In step S1, determine τ at each sampling time f The corresponding combination coefficients α, V s and K v The implementation method of the value is: S1-1, collect a set of joint current data i corresponding to the current sampling time when the robot executes the excitation trajectory in the loaded state e ′, and a set of joint coordinates corresponding to the current sampling time; where each set of joint current data i e 'Includes all joint current data under load at the current sampling time S1-2. According to a set of joint coordinates corresponding to the current sampling moment, construct the linear regression matrix H at the sampling moment l ; S1-3, according to i e ′、H l and χ′, estimate the friction current i at the current sampling time f =i e ′-H l ·[χ′,0] T ; χ′ is given by [χ l χ f ] T The current inertia parameter vector formed by the mass of each link of the robot, the first-order moment of each link, the Coriolis force of each link, the centrifugal force of each link, the gravity of each link, and the moment of inertia of the motor rotor of each joint; S1-4, for τ f Constrain and optimize α, V s and K v , specifically: Among them, θ j is the optimized variable set of the jth joint, θ j =[F c,j F v,j F s,j F b,j α j K v,j V s,j ]; α=[α1,α2,……,α n ], α j is the velocity index of the jth joint; K v =[K v,1 ,K v,2 ,……,K v,j ], K v,j is the inverse tangent coefficient of the jth joint; V s =[V s,1 ,V s,2 ,……,V s,n ],V s,j is the Stribeck velocity of the jth joint; K∈Re n×n is a diagonal matrix of constant coefficients of joint torques, and τ = Ki′ e , τ is the vector composed of all joint torques collected.

4. The method for identifying robot load dynamics parameters based on current dynamics according to claim 1, characterized in that: In step S4, the expression of the constructed kinetic parameter identification model is:

5. The method for identifying robot load dynamics parameters based on current dynamics according to claim 1, characterized in that: In step S4, solve [χ l χ f ] T The OLS solution of l,0 χ f,0 ] T The implementation is:

6. The method for identifying robot load dynamics parameters based on current dynamics according to claim 1, characterized in that: In step S4, according to [ l,0 χ f,0 ] T The implementation method for determining the initial value of the covariance matrix σ is: According to [ l,0 χ f,0 ] T Determine the initial value of the current residual vector R0∈Re Nn×1 ,and Re Nn×1 is a real number space with dimension Nn×1; According to R0, the initial value of the covariance matrix σ is obtained 7. The method for identifying robot load dynamics parameters based on current dynamics according to claim 1, characterized in that: In step S6, 8. The method for identifying robot load dynamics parameters based on current dynamics according to claim 1, characterized in that: In step S8, 9. The method for identifying robot load dynamics parameters based on current dynamics according to claim 1, characterized in that: In step S9, R # ∈Re Nn×1 .

10. The method for identifying robot load dynamics parameters based on current dynamics according to claim 1, characterized in that: In step S10, according to R # The implementation of updating σ is as follows: According to the threshold δ and R # The implementation of updating P is: P=P j ⊙Φ(R # ,d); Among them, Φ(R # ,δ) is a function of a vector of size Nn×1, and R # Elements in that exceed the threshold δ will be set to 0, otherwise they will be set to 1; Update P by the updated P e The implementation is: P e =P·ones(1,c+3n)。

Images