Industrial robot dynamics parameter identification method

Abstract

The industrial robot dynamics parameter identification method comprises the following steps: establishing an industrial robot dynamics model through a Newton-Euler method; designing an excitation trajectory for guiding the movement of the industrial robot; optimizing the excitation trajectory based on a differential evolution algorithm; sampling input data and output data of the industrial robot, and optimizing the input data and output data set through a multi-cycle sampling and averaging method; and separating parameters to be identified of the industrial robot, and identifying the parameters through a convolutional neural network parameter. The method effectively improves the accuracy of the industrial robot dynamics parameter identification, and lays a good foundation for subsequent development of high-precision control of the industrial robot.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of industrial robots, and in particular to a method for identifying the dynamic parameters of an industrial robot. BACKGROUND

[0002] With the development of science and technology, industrial robots have been widely used in many industries, greatly improving productivity; at the same time, people's requirements for the performance of industrial robots are also getting higher and higher, requiring them to achieve more complex functions. Traditional PID control has been unable to meet the demand for high-precision control of industrial robots. Therefore, the method based on the mechanism model of industrial robots to achieve high-precision control will become the mainstream. The model-based industrial robot controller takes the robot dynamics model as the prior value, but the industrial robot is a multivariable, strongly coupled nonlinear system, and it is difficult to obtain the dynamics model through mathematical calculation. In view of this, a large number of researches on the identification of industrial robot dynamics parameters have been carried out at home and abroad, and a variety of effective industrial robot parameter identification methods have been proposed.

[0003] Zu Lianan et al. proposed a neural network-based mechanical arm network model modeling method in the patent application with the application number 202111654616.7 and the name "a neural network mechanical arm dynamics modeling method based on genetic algorithm optimization", which established a mechanical arm network model through a neural network based on genetic algorithm optimization. The established network model accurately reflects the relationship between the joint motion variables and the output torque of the mechanical arm. However, this method does not consider the internal structure of the identification object, but regards the entire system as a "black box", and cannot obtain the specific dynamics parameters inside the identified system. Therefore, it is impossible to obtain an accurate mechanism model of the industrial robot to achieve high-precision control of the industrial robot by relying on this method.

[0004] In the article titled "Industrial Robot Dynamics Parameter Identification and Software System Development", the author Qu Yidan established a dynamics model of the industrial robot through the Newton-Euler method, selected a finite Fourier series as the expression of the optimal excitation trajectory, obtained the parameters in the optimal excitation trajectory by minimizing the condition number of the degenerate matrix, combined the physical constraints of the joint position, velocity and acceleration with the excitation trajectory to make the obtained optimal excitation trajectory physically feasible, avoided the flexible characteristics of the excitation robot by specifying the frequency range in the excitation trajectory, and finally realized the dynamics parameter identification of the mechanical arm by using the least square method. This method can realize the identification of the minimum dynamics parameter set of the industrial robot under the condition that the observation data set is correct, but due to the characteristics of the least square method itself, when there are abnormal data in the observation data set, the identification parameter accuracy will be greatly reduced, resulting in an inaccurate dynamics model.

[0005] In the article entitled "Mechanical arm dynamics parameter identification and verification based on immune algorithm and neural network", the author Pan Qiang expands the traditional excitation trajectory design in the form of Fourier series, which will cause the mechanical arm to produce jitter in the start-stop process, resulting in reduced accuracy of the following trajectory. Through immune clone algorithm optimization, a better excitation trajectory with closed ends and strong noise resistance is obtained, and the neural network is used to identify the dynamics parameters. Compared with the traditional least square identification method, this method effectively improves the accuracy and reliability of the industrial robot dynamics parameter identification. However, the neural network used in this method is too simple, only one hidden layer is set, which leads to poor network performance, and may cause slow convergence or even non-convergence in the training process. SUMMARY

[0006] Therefore, the purpose of the present application is to provide an industrial robot dynamics parameter identification method to solve the problem of identifying key dynamics parameters of industrial robots, to help accurately establish the mechanism model of industrial robots, and to lay a good foundation for subsequent high-precision control of industrial robots. And overcome the technical problems of existing identification methods, such as unreasonable design of optimal excitation trajectory leading to mechanical arm jitter, low reliability of traditional parameter identification method, etc.

[0007] The industrial robot dynamics parameter identification method of the present application comprises the following steps:

[0008] Step one: establish the industrial robot dynamics model by Newton-Euler method as follows:

[0009]

[0010] Where: τ is the joint driving torque vector of the industrial robot; q is the joint angle vector of the industrial robot; is the joint angular velocity vector of the industrial robot; is the joint angular acceleration vector of the industrial robot; D(q) is the inertia matrix of the industrial robot, is the Coriolis force and centrifugal force matrix of the industrial robot, G(q) is the gravity matrix of the industrial robot, τ f is the friction torque;

[0011] Step two: design the excitation trajectory to guide the movement of the industrial robot:

[0012]

[0013] Where: represents the joint angle of a link in the industrial robot; represents the compensation amount of the excitation trajectory, which is a constant; N represents the order of Fourier series, N=5; a j represents the amplitude of the sine term in the Fourier series, b jrepresents the amplitude of cosine term in Fourier series; w f represents the frequency of Fourier series; k represents the number of period sampling; T s represents the sampling period;

[0014] The industrial robot dynamics model is linearized, which is shown as follows:

[0015]

[0016] Wherein: is the observation matrix, and P is a parameter set in the industrial robot system;

[0017] The optimization objective function of the excitation trajectory is designed as:

[0018] fitness=Cond(W)+Pub (4)

[0019] Wherein: fitness is the fitness, W represents a full rank observation coefficient matrix sampled k times in an excitation period, and the expression of W is as follows:

[0020]

[0021] Pub is a penalty function, and its expression is as follows:

[0022]

[0023] When the set excitation trajectory parameters meet the constraint condition, Pub is 0, when the set excitation trajectory parameters do not meet the constraint condition, Pub is far greater than the value of joint angle q i,max , the trajectory will be ignored;

[0024] The constraint conditions include:

[0025] Condition one:

[0026] Condition two:

[0027] Condition three: τ i,min ≤ τ i ≤ τ i,max

[0028] Condition four:

[0029] Wherein: q i,max represents the upper limit of the angle of joint i of the industrial robot connecting rod, q i,min represents the lower limit of the angle of joint i of the industrial robot connecting rod; represents the upper limit of the joint angle velocity of connecting rod i; represents the upper limit of the joint angular acceleration of connecting rod i; qi,init θi (0) represents the initial joint angle of the link i; T is the Fourier series period; τ i,max is the maximum value of the driving torque allowed for the link i, τ i,min is the minimum value of the driving torque allowed for the link i; {s(q(t))} represents the end motion position of the industrial robot at time t; S represents the motion space of the industrial robot;

[0030] The full-rank observation coefficient matrix W with the minimum condition number is obtained based on the differential evolution algorithm, and the condition number is specifically represented as:

[0031] κ(W) = ||W|| ||W -1 || (6)

[0032] Wherein:

[0033]

[0034]

[0035] ||W|| represents the norm of the full-rank observation coefficient matrix W, represents the stretching ability of the full-rank observation coefficient matrix W to the vector, and the larger the matrix norm is, the stronger the stretching ability is; ||W -1 || represents the compression ability of the full-rank observation coefficient matrix W to the vector, and the larger the matrix norm is, the stronger the compression ability is;

[0036] Step three: sampling the input data and the output data (τ) of the industrial robot to obtain the input data set and the output data set;

[0037] Step four: separating the parameters to be identified of the industrial robot to obtain the link i of the industrial robot The parameter vector θ * i :

[0038] θ * i = [m i ,I i ,v i ,c i ] (9)

[0039] Wherein: m i is the mass of the industrial robot mechanical arm link; I i is the moment of inertia at the centroid of the industrial robot mechanical arm link; v i is the viscous friction coefficient at the joint of the industrial robot link i; c i is the Coulomb friction coefficient at the joint of the industrial robot link i;

[0040] The new industrial robot dynamics model is obtained by combining equation (1) and equation (9):

[0041]

[0042] wherein: θ * is an industrial robot parameter vector matrix to be identified;

[0043] Based on formula (10), a neural network function for industrial robot dynamics parameter identification is constructed:

[0044]

[0045] wherein: is an industrial robot parameter vector matrix predicted by a neural network model;

[0046] θ * is identified by using a convolutional neural network, and the convolutional neural network comprises a convolutional layer Conv2d_D1 and a convolutional layer Conv2d_D2, a first BN layer is spliced after the Conv2d_D1, a first SiLU activation function is spliced after the first BN layer, the Conv2d_D2 is spliced after the first SiLU activation function, a second BN layer is spliced after the Conv2d_D2, a second SiLU activation function is spliced after the second BN layer, and the input end of the Conv2d_D2 and the output end of the second SiLU layer are connected by a Shortcut short connection.

[0047] θ is identified by using a convolutional neural network, and the convolutional neural network comprises a convolutional layer Conv2d_D1 and a convolutional layer Conv2d_D2, a first BN layer is spliced after the Conv2d_D1, a first SiLU activation function is spliced after the first BN layer, the Conv2d_D2 is spliced after the first SiLU activation function, a second BN layer is spliced after the Conv2d_D2, a second SiLU activation function is spliced after the second BN layer, and the input end of the Conv2d_D2 and the output end of the second SiLU layer are connected by a Shortcut short connection.

[0048] 1) The input data q obtained by sampling is divided into q by using a Slice operation, and θ * is initialized.

[0049] 2) q and the initialized θ * are input into the convolutional neural network, the output of the SiLU activation function is multiplied by the divided q , and finally θ is obtained by mapping through a linear layer.

[0050] 3) The divided q and θ are input into the convolutional neural network, and finally θ * is obtained by mapping through a linear layer.

[0051] 4) The divided q and the initialized θ * are input into the convolutional neural network, and finally θ is obtained by mapping through a linear layer.

[0052] 5) The divided q With initialized θ * Input the convolutional neural network, input the convolutional neural network, and obtain

[0053] 6) The obtained The neural network function is obtained by concatenating

[0054] 7) The gradient descent method is used again to identify the industrial robot dynamics parameters The loss function uses a cross-entropy loss function:

[0055]

[0056] The gradient of the loss function L with respect to is:

[0057]

[0058] The parameter matrix update rule is:

[0059]

[0060] Where: n represents the number of iterations; ε is the learning rate, is the momentum parameter, g n is the gradient of each iteration; and the finally identified industrial robot dynamics parameters are:

[0061]

[0062] Further, the step two) of obtaining the observation matrix W with the smallest condition number based on the differential evolution algorithm includes the following steps:

[0063] 1) Initialize the population:

[0064] Set the excitation trajectory parameters contained in the individuals in the population a j and b j , and set the initial population size;

[0065] 2) Calculate the fitness:

[0066] According to formula (4), the fitness of the individuals in the population is calculated, which provides a basis for subsequent comparison and selection operations;

[0067] 3) Compare the results:

[0068] Compare the fitness of the individuals in the population with the set optimal solution, if the requirement is met in advance, stop iteration, and derive the optimal solution;

[0069] 4) Mutation operation

[0070] The standard differential evolution algorithm is used to mutate the population, and the standard differential evolution algorithm is represented as:

[0071] V i,g+1 = X r1,g +F*(X r2,g -X r3,g ) (16)

[0072] Wherein: V i,g+1 represents the next generation of individuals produced after mutation; X r1,g , X r2,g and X r3,g represent individuals randomly taken from the old population; F represents the mutation factor, and the value of the mutation factor determines the convergence of the differential evolution algorithm;

[0073] 5) crossover operation:

[0074] The new individuals produced after mutation are combined with the old individuals in the population to produce new individuals to realize the crossover operation;

[0075] 6) selection operation:

[0076] After the mutation operation and the crossover operation, a new individual U i,g+1 is generated, the fitness of the new individual is calculated, and the new individual is compared with the old population individuals to determine whether the new individual can be retained to the next generation. The specific algorithm of the selection operation is:

[0077]

[0078] 7) iteration operation:

[0079] Continue to step 2) repeat the above steps until there is an individual in the population that meets the optimal result or the iteration number reaches the maximum iteration number g max .

[0080] Further, the step three) further comprises averaging the input data and output data (τ) of the industrial robot for M times of periodic sampling to eliminate measurement errors; the average driving torque average angle average angular velocity average angular acceleration of the connecting rod i at time k obtained by sampling are:

[0081]

[0082]

[0083]

[0084]

[0085] Wherein: M represents the number of sampling periods.

[0086] Advantages of the present application:

[0087] 1. In view of the problems that the traditional least square regression parameter identification method is greatly interfered by abnormal data and has poor stability, and the general neural network model identification cannot obtain specific dynamic parameters, the industrial robot dynamic parameter identification method provided by the present application is a convolutional neural network parameter identification method based on vector segmentation and recombination, which effectively improves the accuracy of industrial robot dynamic parameter identification and lays a good foundation for subsequent high-precision control of industrial robots.

[0088] 2. In view of the fact that the traditional excitation trajectory design does not consider the shaking of the industrial robot during the excitation process and lacks necessary constraint conditions, resulting in a serious shaking phenomenon during the excitation process, the industrial robot dynamic parameter identification method of the present application adds the joint speed and acceleration limits of the robot at the beginning and the end of the cycle to the constraint conditions, effectively reduces the shaking of the robot during the excitation process, and improves the reliability of parameter identification.

[0089] 3. The industrial robot dynamic parameter identification method of the present application designs the excitation trajectory as a five-order Fourier series and optimizes the trajectory through a differential evolution algorithm, which can provide suitable excitation trajectory parameters for the industrial robot. By using the periodic nature of the Fourier series and through the method of multi-cycle sampling and averaging, the influence of noise and other interference factors on the observation data set is effectively reduced, the accuracy of the observation data set is greatly improved, and the accuracy of parameter identification is further improved. BRIEF DESCRIPTION OF DRAWINGS

[0090] Figure 1 is a schematic diagram of an industrial robot connection;

[0091] Figure 2 is a flowchart of a differential evolution algorithm;

[0092] Figure 3 is a convolutional neural network structure diagram for industrial robot parameter identification;

[0093] Figure 4 is a flowchart of industrial robot dynamic parameter identification. DETAILED DESCRIPTION

[0094] The present application will be further described below in combination with the drawings and examples.

[0095] The industrial robot dynamic parameter identification method in the present embodiment includes the following steps:

[0096] Step one: an industrial robot dynamic model is established by Newton-Euler method as follows:

[0097]

[0098] Where: τ is the joint driving torque vector of the industrial robot; q is the joint angle vector of the industrial robot; This refers to the angular velocity vector of the industrial robot joints. Let be the joint angular acceleration vector of the industrial robot; D(q) is the inertia matrix of the industrial robot. Let G(q) be the Coriolis force and centrifugal force matrix of the industrial robot, and let τ be the gravity matrix of the industrial robot. f This is the frictional torque.

[0099] The following is based on Figure 1 Taking the articulated robot shown as an example, the steps for establishing an industrial robot dynamics model using the Newton-Euler method are explained in detail:

[0100] Figure 1 The diagram shows a link connection of an industrial robot. A coordinate system is established at each link node, and the transformation matrix between adjacent coordinate systems is obtained using the DH modeling method. This enables the kinematic correlation between links in an industrial robot system.

[0101] According to Newton's third law, the net force f acting on the center of mass of link i in an industrial robot is... ic for:

[0102] f ic =m i a ic ①

[0103] Where: m is the mass of the connecting rod; a ic Let be the linear acceleration at the center of mass of link i.

[0104] According to Euler's equations, the net torque on the center of mass of link i in the industrial robot is... c n ic for:

[0105]

[0106] in: c I i Let ω be the inertial tensor in the coordinate system {c} established with the center of mass as the origin for robot link i; i Let be the angular velocity of robot link i relative to the base coordinate system; Let be the angular acceleration of robot link i relative to the polar coordinate system.

[0107] By using the Newton-Euler method in a forward recursion, the angular velocity of the i-th link of the industrial robot arm relative to coordinate system i can be obtained. i ω iangular acceleration linear acceleration i a i The iterative relationship between i-1 and i:

[0108] angular velocity i ω i The iterative relationship:

[0109]

[0110] angular acceleration The iterative relationship:

[0111]

[0112] linear acceleration i a i The iterative relationship:

[0113]

[0114] wherein: represents the rotation matrix between the i-1 coordinate system and the i coordinate system, which describes the transformation relationship between the two coordinate systems; represents the joint angular velocity of the robot link i; represents the projection direction of the vector in the coordinate system i to the Z-axis direction in the coordinate system i; represents the joint angular acceleration of the robot link i; i-1 P i represents the vector from the origin of the coordinate system i-1 to the origin of the coordinate system i.

[0115] By Newton-Euler method, the joint moment of the robot link i is obtained by inverse recursion i n i and the iterative relationship between the joint moment of the link i+1 i+1 n i+1

[0116]

[0117] wherein: i r ci represents the distance from the robot link i to the joint.

[0118] The obtained joint moment i n i is projected to the Z-axis direction in the coordinate system i, and the driving moment τ of the industrial robot link i is obtained i :

[0119]

[0120] ​Simultaneous equations ①-⑦, taking a three-joint RRR industrial robot as an example, the dynamic equation thereof is established:

[0121]

[0122] Wherein: τ = [τ1, τ2, τ3] T is the joint driving torque vector of the industrial robot; q = [θ1, θ2, θ3] is the joint angle vector of the industrial robot; is the joint angular velocity vector of the industrial robot; is the joint angular acceleration vector of the industrial robot.

[0123] D(q) is the inertia matrix of the industrial robot, which is specifically represented as a 3-order symmetric matrix as follows:

[0124]

[0125] Wherein:

[0126] D 11 = I2+I3+m2r2 2 +2m3r3l2c2c 2 +m3r3 23 c 2 23 3 +m2l2c2 2

[0127] D 22 = I2+I3+m2r2 2 +2m3r3l2c3+m3r3 2 +m3l2 2

[0128] D 33 = I3+m3r3 2

[0129] D 23 = D 32 = I3+m3r3 2 +m3r3l2c3

[0130] D 12 = D 13 = D 21 = D 31 = 0

[0131] is the Coriolis force and centrifugal force matrix of the industrial robot system, which is specifically represented as a 3-dimensional column vector:

[0132]

[0133] Wherein:​

[0134]

[0135]

[0136]

[0137] G(q) is the gravity matrix of the industrial robot system, which is specifically a 3-dimensional column vector:

[0138]

[0139] Wherein:

[0140] G1 = [(m2r2+m3l2)c2+m3r3c 23 ]s1g

[0141] G2 = [(m2r2+m3l2)s2+m3r3s 23 ]c1 2 g

[0142] G3 = m3r3c1 2 s 23 g

[0143] m2, m3 are the masses of the industrial robot arm links 2, 3; I2, I3 are the moments of inertia at the centers of mass of the industrial robot arm links 2, 3; r2, r3 are the distances from the centers of mass of the industrial robot arm links 2, 3 to the joints.

[0144] The dynamic model of formula 7 does not consider the effect of friction force, in order to make the mechanism model more realistic, the friction torque item τ f is added, which is specifically:

[0145]

[0146] Wherein: f C is the Coulomb friction coefficient matrix; f v is the viscous friction coefficient matrix.

[0147] Thus, a complete industrial robot dynamic model is established:

[0148]

[0149] Step two: only the appropriate excitation trajectory is designed, the correct input output (τ) data can be obtained, so as to improve the parameter identification accuracy. This embodiment adopts five-order Fourier series as the parameter identification excitation trajectory of the industrial robot.

[0150] Designing excitation trajectory to guide industrial robot motion:

[0151]

[0152] wherein: q represents a link joint angle in industrial robot; represents a compensation amount of excitation trajectory, which is a constant; N represents Fourier series order, N=5; a j represents a sine term amplitude in Fourier series, b j represents a cosine term amplitude in Fourier series; w f represents Fourier series frequency; k represents periodic sampling number; T s represents sampling period; w f j in jkT and a j and b j have the same meaning, which is Fourier series order.

[0153] Analyzing the industrial robot dynamics model established in step 1), it can be known that four components of driving torque τ are all related to industrial robot link angle q, angular velocity angular acceleration , and all contain all inertia parameters of industrial robot system, so the industrial robot dynamics model can be linearized, and the specific form is as follows:

[0154]

[0155] wherein: is an observation matrix, and P is a parameter set in industrial robot system.

[0156] The optimization objective function of excitation trajectory design is:

[0157] fitness=Cond(W)+Pub (4)

[0158] wherein: fitness is fitness, W represents a full rank observation coefficient matrix sampled k times in one excitation period, and the expression of W is as follows:

[0159]

[0160] Pub is a penalty function, and its expression is as follows:

[0161]

[0162] When the set excitation trajectory parameters meet the constraint condition, Pub is 0, when the set excitation trajectory parameters do not meet the constraint condition, Pub is far greater than a value of joint angle q i,max , and the trajectory will be ignored. The constraint conditions include:

[0163] Condition one:

[0164] Condition two:

[0165] Condition three: τ i,min ≤ τ i ≤ τ i,max

[0166] Condition four:

[0167] Wherein: q i,max represents the upper limit of the angle of the joint of the industrial robot link i, q i,min represents the lower limit of the angle of the joint of the industrial robot link i; represents the upper limit of the joint angle velocity of the link i; represents the upper limit of the joint angular acceleration of the link i; q i,init represents the initial joint angle of the link i; T is the Fourier series period; τ i,max is the maximum value of the driving torque allowed for the link i, τ i,min is the minimum value of the driving torque allowed for the link i; {s(q(t))} represents the end motion position of the industrial robot at time t; S represents the motion space of the industrial robot;

[0168] Based on the differential evolution algorithm, a full-rank observation coefficient matrix W with the smallest condition number is obtained, and the condition number is specifically represented as:

[0169] κ(W) = ||W|| ||W -1 || (6)

[0170] Wherein:

[0171]

[0172]

[0173] ||W|| represents the norm of the full-rank observation coefficient matrix W, represents the stretching ability of the full-rank observation coefficient matrix W to the vector, and the larger the matrix norm is, the stronger its stretching ability is; ||W -1 || represents the compression ability of the full-rank observation coefficient matrix W to the vector, and the larger the matrix norm is, the stronger its compression ability is.

[0174] The identification principle of the industrial robot system is to design an excitation trajectory, sample k times within one excitation period kT s , and based on the established industrial robot dynamics model and the collected input and output data, use the identification algorithm to deduce the internal parameters of the system. After sampling k times, we can get:

[0175] Y = W·P (I)

[0176] wherein: Y and W represent joint torque matrix and full rank observation coefficient matrix respectively, which are sampled k times in one excitation period:

[0177]

[0178]

[0179] Considering that there are noise and other interference factors in actual observation, which cause the error of the identified parameter set P to be large, in order to improve the stability of the linear equation Y=W·P, the observation matrix condition number needs to be minimized by setting reasonable excitation trajectory parameters, and the following is inferred:

[0180] The observation matrix W is a full rank matrix, let W -1 R, then the equation Y=W·P can be rewritten as:

[0181] RY=P (II)

[0182] Adding noise δY to the right side of equation , the following equation can be obtained:

[0183] R(Y+δY)=P+δP (III)

[0184] RδY=δP (IV)

[0185] Taking the norm of both sides of equation (II):

[0186] ||RδY||=||δP||(V)

[0187] Substitute equation (7) and equation (8):

[0188] ||R||||δY||≥||δP|| (VI)

[0189] ||W||||P||≥||Y|| (VII)

[0190] Multiply equation (VI) and equation (VII):

[0191]

[0192] Equation (VIII) shows the lower bound of the influence of observation data Y on the identified parameter P, and the larger the condition number k(W) is, the larger the lower bound deviation range of P is. According to the above derivation method, the following equation can be obtained:

[0193]

[0194] Equation (IX) shows the upper bound of the influence of observation data Y on the identified parameter P, and the larger the condition number k(W) is, the larger the upper bound deviation range of P is.

[0195] Simultaneous equations (VIII) and (IX):

[0196]

[0197] According to formula (X), the identification parameter deviation depends on the condition number of the full rank observation coefficient matrix W. In order to ensure the accuracy of the identification parameters, the Fourier series parameters should be reasonably set when designing the excitation trajectory a j , b j Make k(W) as small as possible under the constraint condition.

[0198] This embodiment is based on the differential evolution algorithm to obtain the full rank observation coefficient matrix W with the minimum condition number. The differential evolution algorithm is a heuristic random search algorithm based on the difference between groups. The algorithm is similar to the genetic algorithm process. New individuals are generated from the existing population through selection, compilation and crossover operation. The fitness (fitness) is used to judge the individual. The lower the fitness value, the stronger the adaptability of the individual, and the better the individual. In each iteration, the individuals with poor adaptability are eliminated to obtain the optimal solution. The differential evolution algorithm process is shown in Figure 2 Based on the differential evolution algorithm, the observation matrix W with the minimum condition number includes the following steps:

[0199] 1) Initialize the population:

[0200] Set the excitation trajectory parameters contained in the individual in the population a j and b j , and set the initial population size.

[0201] 2) Calculate the fitness:

[0202] According to formula (4), the fitness of the individual in the population is calculated, which provides a basis for subsequent comparison and selection operation;

[0203] 3) Compare the results:

[0204] Compare the fitness of the individual in the population with the set optimal solution. If the requirement is met in advance, stop iteration and derive the optimal solution;

[0205] 4) Mutation operation

[0206] The standard differential evolution algorithm is used for mutation of the population, and the standard differential evolution algorithm is represented as:

[0207] V i,g+1 =X r1,g +F*(X r2,g -X r3,g ) (9)

[0208] Where: Vi,g+1 represents the next generation of individuals produced after mutation; X r1,g , X r2,g and X r3,g represents the individual randomly taken out in the old population; F represents the mutation factor, whose value determines the convergence of the differential evolution algorithm, and if F is too small, the optimization algorithm may diverge, and if F is too large, the convergence may be too slow, and in the embodiment, F is taken as 0.7;

[0209] 5) crossover operation:

[0210] The new individual produced after mutation is combined with the old individual in the population to produce a new individual to realize the crossover operation;

[0211] 6) selection operation:

[0212] After the mutation operation and the crossover operation, a completely new individual U i,g+1 is produced, the fitness of the new individual is calculated, and the new individual is compared with the old population individual to determine whether the new individual can be retained to the next generation, and the selection operation is specifically as follows:

[0213]

[0214] 7) iteration operation:

[0215] Continue to step 2) repeat the above steps until there is an individual in the population that meets the optimal result or the iteration number reaches the maximum iteration number g max .

[0216] Through the above operation, the observation matrix with the minimum condition number can be obtained, so that the key parameters in the five-order Fourier series are obtained, and the optimization of the excitation trajectory is realized.

[0217] Step three: sampling the input data and output data (τ) of the industrial robot to obtain the input data and output data set.

[0218] Step four: separating the parameters to be identified of the industrial robot to obtain the parameter vector θ * i to be identified of the link i of the industrial robot

[0219] θ * i = [m i , I i , v i , c i ] (11)

[0220] wherein: m i is the mass of the industrial robot arm link; I iis the rotational inertia of the industrial robot mechanical arm link at the center of mass; v i is the viscous friction coefficient of the industrial robot link i at the joint; c i is the Coulomb friction coefficient of the industrial robot link i at the joint.

[0221] The new industrial robot dynamics model is obtained by combining equations (1) and (11) as follows:

[0222]

[0223] where: θ * is the vector matrix of the industrial robot parameters to be identified;

[0224] Based on equation (12), a neural network function for identifying the dynamics parameters of the industrial robot is constructed as follows:

[0225]

[0226] where: is the vector matrix of the industrial robot parameters predicted by the neural network model.

[0227] The convolutional neural network is used to identify θ * , which includes a convolutional layer Conv2d_D1 and a Conv2d_D2, a first BN layer spliced after Conv2d_D1, a first SiLU activation function spliced after the first BN layer, a Conv2d_D2 spliced after the first SiLU activation function, a second BN layer spliced after Conv2d_D2, a second SiLU activation function spliced after the second BN layer, and a Shortcut short-circuit between the input end of Conv2d_D2 and the output end of the second SiLU layer.

[0228] First, the convolutional neural network is used to identify the parameters of , including the following steps:

[0229] 1) The input data obtained by sampling is divided into q, Initialize θ * ;

[0230] 2) Input q and the initialized θ * into the convolutional neural network, multiply the output of the SiLU activation function by the divided , and finally map through a linear layer.

[0231] 3) Input the divided q and , and the initialized θ * into the convolutional neural network, and finally map

[0232] 4) The segmented q is combined with the initialized θ * The input is mapped to the output of the convolutional neural network through a linear layer

[0233] 5) The segmented is combined with the initialized θ * The input is mapped to the output of the convolutional neural network through a linear layer

[0234] 6) The obtained is concatenated to obtain the neural network function

[0235] 7) Gradient descent is used to identify the dynamics parameters of the industrial robot The loss function uses the cross-entropy loss function:

[0236]

[0237] The gradient of the loss function L with respect to is:

[0238]

[0239] The parameter matrix update rule is:

[0240]

[0241] Where: n represents the number of iterations; ε is the learning rate, is the momentum parameter, g n is the gradient of each iteration; the identified dynamics parameters of the industrial robot are:

[0242]

[0243] The designed excitation trajectory obtained from step two is a five-order Fourier series, so the angular velocity and angular acceleration of link i joint are expressed as:

[0244]

[0245]

[0246] From equations (18) and (19), it is easy to know that the angular velocity and angular acceleration of each joint of the industrial robot are still periodic functions. In theory, the angular velocity and angular acceleration The same time in each cycle is the same, but due to the presence of noise and other interference factors, resulting in measurement error. For this improvement to the above embodiment, the step three) further comprises M cycle sampling input data of industrial robot And the output data (τ) are averaged to eliminate the measurement error; the average driving torque of the connecting rod i at time k obtained by sampling Average angle Average angular velocity Average angular acceleration Is:

[0247]

[0248]

[0249]

[0250]

[0251] Wherein: M represents the number of sampling cycles.

[0252] Finally, it should be pointed out that the above embodiments are only used to illustrate the technical solutions of the present application and are not limiting. Although the present application has been described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that the technical solutions of the present application can be modified or replaced by equivalents without departing from the spirit and scope of the present application, and they should be covered in the scope of the claims of the present application.

Claims

1. A method for identifying dynamics parameters of an industrial robot, characterized in that: The method comprises the following steps: Step one: establish the industrial robot dynamics model by Newton-Euler method as follows: wherein: τ is a joint driving torque vector of the industrial robot; q is a joint angle vector of the industrial robot; is a joint angular velocity vector of the industrial robot; is a joint angular acceleration vector of the industrial robot; D(q) is an inertia matrix of the industrial robot, is a Coriolis and centrifugal force matrix of the industrial robot, G(q) is a gravity matrix of the industrial robot, τ f is a friction torque; Step two: design the excitation trajectory to guide the industrial robot movement: wherein: θ represents a link joint angle in an industrial robot; represents a compensation amount of the excitation trajectory, which is a constant; N represents a Fourier series order, N = 5; a j represents a sine term amplitude in the Fourier series; b j represents a cosine term amplitude in the Fourier series; w f represents a Fourier series frequency; k represents a periodic sampling number; T s represents a sampling period; Linearize the industrial robot dynamics model, and the specific form is as follows: wherein: is an observation matrix, P is a parameter set in the industrial robot system; The optimization objective function of the excitation trajectory is: fitness=Cond(W)+Pub (4) Wherein: fitness is the fitness, W represents a full rank observation coefficient matrix sampled k times in an excitation period, and the expression of W is as follows: Pub is a penalty function, and its expression is as follows: When the set excitation trajectory parameter meets the constraint condition, Pub is 0, and when the set excitation trajectory parameter does not meet the constraint condition, Pub is far greater than a value of the joint angle q i,max , and the trajectory will be ignored. The constraint conditions include: Condition one: Condition two: Condition three: τ i,min ≤ τ i ≤ τ i,max Condition Four: wherein: q i,max represents the upper limit of the angle of the joint of the industrial robot link i, q i,min represents the lower limit of the angle of the joint of the industrial robot link i; represents the upper limit of the joint angular velocity of the link i; represents the upper limit of the joint angular acceleration of the link i; q i,init represents the initial joint angle of the link i; T is the period of the Fourier series; τ i,max is the maximum value of the driving torque allowed for the link i, τ i,min is the minimum value of the driving torque allowed for the link i; {s(q(t))} represents the end motion position of the industrial robot at time t; S represents the motion space of the industrial robot; Based on the differential evolution algorithm, the full rank observation coefficient matrix W with the minimum condition number is obtained, and the condition number is specifically represented as: K(W) = ||W|| ||W -1 || (6) Wherein: ||W|| represents the norm of the full-rank observation coefficient matrix W, represents the stretching ability of the full-rank observation coefficient matrix W to the vector, and the larger the matrix norm, the stronger the stretching ability; ||W -1 || represents the compression ability of the full-rank observation coefficient matrix W to the vector, and the larger the matrix norm, the stronger the compression ability; Step three: sampling input data of industrial robot and output data (τ) to obtain input data and output data sets; Step four: separate the parameters to be identified by the industrial robot, and obtain the parameter vector θ to be identified of the connecting rod i of the industrial robot * i : θ * i = [m i , I i , v i , c i ] (9) where: m i is the industrial robot manipulator link mass; I i is the industrial robot manipulator link mass center moment of inertia; v i is the industrial robot link i joint viscous friction coefficient; c i is the industrial robot link i joint coulomb friction coefficient; The new industrial robot dynamics model obtained by combining equation (1) and equation (9) is: wherein: θ * is a vector matrix of parameters to be identified for the industrial robot; Based on equation (10), the neural network function for industrial robot dynamics parameter identification is constructed: wherein: is a matrix of parameter vectors of the industrial robot predicted by the neural network model; A convolutional neural network is used to recognize θ * , which includes a convolution layer Conv2d_D1 and a convolution layer Conv2d_D2. The Conv2d_D1 is spliced with a first BN layer, the first BN layer is spliced with a first SiLU activation function, the first SiLU activation function is spliced with the Conv2d_D2, the Conv2d_D2 is spliced with a second BN layer, the second BN layer is spliced with a second SiLU activation function, and the input end of the Conv2d_D2 and the output end of the second SiLU layer are connected through a Shortcut. First, a convolutional neural network is used to perform parameter identification, including the following steps: 1) Input data sampled by Slice operation Slice into q, Initialize θ * ; 2) set q and initialized theta * The input convolutional neural network multiplies the output of the SiLU activation function with the split ​ 3) q and and initialized θ * The input convolutional neural network is finally mapped by a linear layer 4) q is divided by the initialized theta * The input convolutional neural network is finally mapped by a linear layer 5) cut the with the initialized θ * inputting the convolutional neural network, inputting the convolutional neural network, obtaining 6) The resulting Concatenate the resulting neural network functions 7) Again, the gradient descent method is used to identify the dynamics parameters of industrial robots The loss function adopts a cross-entropy loss function: The gradient of the loss function L with respect to is: The parameter matrix update rule is: where: n represents the number of iterations; ε is the learning rate, is a momentum parameter, g n is the gradient for each iteration; The finally identified industrial robot dynamics parameters are:

2. The industrial robot dynamics parameter identification method according to claim 1, characterized in that: The step two based on the differential evolution algorithm to obtain the observation matrix W with the minimum condition number comprises the following steps: 1) Initialize the population: setting the parameters of the excitation trajectories contained by the individuals in the population a j and b j and setting the number of the initial population; 2) Calculate the fitness: Calculate the fitness of the individuals in the population according to formula (4), which provides the basis for subsequent comparison and selection operation; 3) Compare the results: Compare the fitness of the individuals in the population with the fitness corresponding to the set optimal solution, if the requirement is reached in advance, stop iteration, and export the optimal solution; 4) Mutation operation The population is mutated by using the standard differential evolution algorithm, and the standard differential evolution algorithm is represented as: V i,g+1 = X r1,g + F*(X r2,g - X r3,g ) (16) where: V i,g+1 represents the next generation of individuals produced after variation; X r1,g , X r2,g and X r3,g represent individuals randomly taken from the old population; F represents a variation factor, the value of which determines the convergence of the differential evolution algorithm; 5) Cross operation: Combine the new individuals generated after mutation with the old individuals in the population to generate new individuals to realize cross operation; 6) Selection operation: After the mutation operation and the crossover operation, a new individual U is generated i,g+1 The fitness of the new individual is calculated, compared with the old population individuals, and judged whether the new individual can be reserved to the next generation. The selection operation specific algorithm is: 7) Iteration operation: Continuing to step 2) repeat the above steps until an individual in the population of individuals meets the optimal result or the number of iterations reaches a maximum number of iterations g max .

3. The industrial robot dynamics parameter identification method according to claim 1, characterized in that: The step three further comprises averaging the input data sampled M times for the industrial robot and the output data (τ) to eliminate measurement errors; the average driving moment of the link i at time k resulting from the sampling average angle average angular velocity average angular acceleration is: Wherein: M represents the number of sampling periods.

Images