A method and system for online correction of trajectory accuracy of an electric discharge machining robot

Abstract

The application discloses an electric spark machining robot track precision online correction control method and relates to the technical field of robot control.The application collects multidimensional data such as joint motion, inter-electrode discharge state and electrode cumulative loss in real time to construct a unified state vector. The double-mechanism parameter identification method of the application updates inertia Coriolis force parameters online by combining a first mechanism with a forgetting factor, optimizes friction parameters by using a second mechanism through a particle swarm algorithm during a stable motion period to obtain a joint estimation dynamics vector. The track deviation type is determined based on the dynamics vector. If the track deviation type is discharge mutual interference or strong model drift, a track correction strategy that fuses space geometry and process state is executed. If the track deviation type is weak drift, torque compensation based on model uncertainty is implemented. The application realizes dynamic approximation of dynamics parameters, effectively balances track tracking precision and discharge gap stability and has the advantages of high precision, high stability and the like.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of robot control, and in particular to a method and system for online correction control of the trajectory precision of an electric spark machining robot. BACKGROUND

[0002] Electric spark machining is a special processing technology, which is widely used in processing high-hardness and complex curved surface parts. With the development of intelligent manufacturing, industrial robots are gradually introduced into the field of electric spark machining to replace traditional special machine tools to perform complex processing tasks due to their high degree of freedom, large working space, strong flexibility and relatively low cost. However, industrial robots essentially belong to a serial cantilever beam structure, and their rigidity is much lower than that of traditional numerical control machine tools. In addition, there are problems such as backlash, nonlinear friction and flexible deformation in the joint transmission system. When electric spark machining is performed, the robot needs to maintain a very high precision of end trajectory tracking capability to maintain a micron-level inter-electrode discharge gap. Once the actual trajectory deviates from the theoretical trajectory, it is easy to cause unstable discharge, short circuit or arc, which seriously affects the processing quality and efficiency. In addition, electric spark machining is a dynamic process accompanied by electrode wear. The mass distribution, inertia and joint friction characteristics of the tool electrode will change over time with the processing time and motion state. The traditional control method based on a fixed dynamic model is difficult to adapt to such dynamic changes, resulting in a decrease in trajectory precision over time.

[0003] The existing robot trajectory correction method mainly focuses on the field of contact processing. The patent document CN115213906B discloses a robot trajectory correction method and device. This technology establishes a robot dynamics model, calculates the theoretical force parameters of each joint motor in real time, and sets the joint floating parameters to make the robot joints passively float under the action of external force to adapt to the irregular changes of the workpiece surface, avoid overload and realize flexible contact. The above existing technology has obvious limitations when applied to high-precision electric spark machining scenes. Electric spark machining requires strict non-contact state control, i.e., maintaining a constant inter-electrode discharge gap, rather than non-contact floating. The existing technology relies on contact force feedback to achieve passive flexible adjustment, which will cause frequent short circuit backoff or excessive gap in electric spark machining, and cannot actively ensure the stability of the discharge gap. The dynamic model parameters of the existing technology are usually offline calibrated or fixed, and the unique time-varying factors in the electric spark machining process are not considered. The existing technology lacks precise identification of deviation types and targeted hierarchical correction strategies, making it difficult for the robot to simultaneously consider trajectory tracking precision and processing process stability in complex electric spark discharge environments. SUMMARY

[0004] The technical problem solved by the present application is that the prior art has obvious limitations when applied to high-precision electric spark machining scenes. Electric spark machining requires strict non-contact state control, i.e., maintaining a constant inter-electrode discharge gap, and the prior art relies on the feedback of contact force to achieve passive flexible adjustment, which can cause frequent short-circuit backoff or excessive gap in electric spark machining, and cannot actively ensure the stability of the discharge gap. The dynamic model parameters of the prior art are usually offline calibrated or fixed, and the unique time-varying factors in the electric spark machining process are not considered. The prior art lacks precise identification of deviation types and targeted hierarchical correction strategies, making it difficult for robots to simultaneously consider trajectory tracking accuracy and process stability in complex electric spark discharge environments.

[0005] To solve the above technical problems, the present application provides the following technical solutions: an electric spark machining robot trajectory accuracy online correction control method, comprising the following steps:

[0006] Step S1: running the machining trajectory and collecting multi-dimensional data;

[0007] Step S2: constructing a unified state vector according to the multi-dimensional data, obtaining an observation matrix, and updating the to-be-identified dynamic parameters using a double-mechanism parameter identification method combined with the observation matrix;

[0008] Step S3: comparing the theoretical running trajectory with the actual machining trajectory to obtain the trajectory deviation, determining the type of the deviation, and executing the trajectory online correction strategy according to the type of the deviation based on the updated to-be-identified dynamic parameters.

[0009] Preferably, step S1 comprises:

[0010] The multi-dimensional data is preprocessed joint motion data, inter-electrode discharge data, electrode wear data, and joint driving torque;

[0011] The joint motion data is the angle, angular velocity, and angular acceleration of the joint;

[0012] The inter-electrode discharge data includes inter-electrode voltage, machining current, pulse utilization rate, and inter-electrode gap;

[0013] The electrode wear data includes cumulative wear length calculated in real time based on a discharge energy loss model;

[0014] The joint driving torque is the product of the torque constant of the motor, the reduction ratio, and the armature current input to the motor;

[0015] The preprocessing includes denoising and outlier rejection;

[0016] Based on the preprocessed multi-dimensional data, a unified state vector is constructed;

[0017] The unified state vector comprises angles, angular velocities, angular accelerations, joint driving torques, inter-electrode gaps and cumulative loss lengths of the joints of the robot after pre-processing.

[0018] Preferably, the dual-mechanism parameter identification method comprises a first identification mechanism and a second identification mechanism, and the first identification mechanism comprises:

[0019] An observation matrix at the current time is obtained according to the angles, angular velocities and angular accelerations of the joints in the unified state vector;

[0020] A condition number of the observation matrix at the current time is calculated, and a dynamic update gain factor is calculated according to the condition number;

[0021] The mathematical expression of the dynamic update gain factor is:

[0022] ;

[0023] Wherein, is the dynamic update gain factor, is a preset adjustment constant, is the condition number of the observation matrix, is the observation matrix;

[0024] A forgetting factor is calculated based on the inter-electrode gap, and the mathematical expression of the forgetting factor is:

[0025] ;

[0026] Wherein, is the difference between the inter-electrode gap at the current time and the inter-electrode gap at the previous time, is the forgetting factor, is a reference gap reference value, is a sensitivity coefficient;

[0027] A gain matrix is calculated according to the observation matrix at the current time and the forgetting factor, and the mathematical expression of the gain matrix is:

[0028] ;

[0029] Wherein, is the forgetting factor, is the observation matrix, is a covariance matrix at the previous time, is the transpose matrix of the observation matrix.

[0030] The dual-mechanism parameter identification method comprises a first identification mechanism and a second identification mechanism, and the first identification mechanism further comprises:

[0031] an inertial Coriolis force parameter vector to be identified is updated on line by using a recursive least square method based on a gain matrix, a joint driving moment and an observation matrix at a previous time, the inertial Coriolis force parameter vector to be identified is composed of a mass of each connecting rod, a mass center position of each connecting rod and an inertial tensor of each connecting rod;

[0032] The mathematical expression for updating the inertial Coriolis force parameter vector to be identified is:

[0033] ;

[0034] wherein, is the inertial Coriolis force parameter vector to be identified, is a gain matrix, is a joint driving moment, is the inertial Coriolis force parameter vector to be identified at a previous time, is a dynamic updating gain factor, is a moment error vector, and the mathematical expression for the moment error vector is:

[0035] ;

[0036] wherein, is a joint driving moment vector measured at a current time, is an observation matrix at a current time, is the inertial Coriolis force parameter vector to be identified at a previous time;

[0037] The covariance matrix is initialized as a product of a preset constant and a unit matrix;

[0038] The mathematical expression for the covariance matrix is:

[0039] ;

[0040] wherein, is a forgetting factor, is an observation matrix, is a gain matrix, is a covariance matrix at a previous time.

[0041] Preferably, the step S2 further comprises a second identification mechanism, and the second identification mechanism comprises:

[0042] If the average condition number of the observation matrix continuously exceeds a preset condition number threshold within a preset time window, a parameter identification step of the friction parameter to be identified is started;

[0043] The parameter identification step adopts a particle swarm optimization algorithm, and specifically comprises:

[0044] initializing a particle swarm in a preset feasible region of friction parameters, wherein each particle represents a candidate friction parameter vector composed of a Coulomb friction coefficient and a viscous friction coefficient, and each particle is randomly assigned an initial velocity;

[0045] calculating a fitness function of each particle by a fitness function, and a mathematical expression of the fitness function is:

[0046]

[0047] wherein, is a residual torque vector, is a model predicted torque calculated based on the friction parameter vector represented by the current particle, is a weight, is an observation matrix at the kth moment;

[0048] updating the velocity and position of all particles according to the individual optimal position of each particle in the iteration process and the global optimal position of all particles in the iteration process, and iterating until a maximum iteration number is satisfied to obtain an optimized friction parameter vector.

[0049] Preferably, the step S2 further comprises:

[0050] constructing a total dynamic parameter vector composed of inertial Coriolis force parameters and friction parameters;

[0051] the inertial Coriolis force parameter part is an initialized inertial Coriolis force parameter vector to be identified;

[0052] the friction parameter part is an initialized friction parameter to be identified;

[0053] in each control period, a first identification mechanism is executed to obtain an updated inertial Coriolis force parameter to be identified, and the inertial Coriolis force parameter to be identified is used to replace the corresponding part of the inertial Coriolis force parameter to be identified in the total dynamic parameter vector, and the friction parameter part remains unchanged;

[0054] when a parameter identification step of a second identification mechanism is triggered, the second identification mechanism is executed to obtain an optimized friction parameter vector, and the optimized friction parameter vector is used to replace the corresponding friction parameter part in the total dynamic parameter vector, and the inertial Coriolis force parameter part remains unchanged, thereby obtaining a joint estimated dynamic parameter vector.

[0055] Preferably, the step S3 comprises:

[0056] ​The difference between the theoretical end effector pose and the actual end effector pose is taken as a trajectory deviation, and it is determined whether the trajectory deviation exceeds a preset dynamic trajectory reference value, and if the trajectory deviation exceeds the preset dynamic trajectory reference interval threshold, a deviation type determination is performed;

[0057] When the trajectory deviation exceeds the preset dynamic trajectory reference threshold;

[0058] If the pulse utilization rate is less than 85%, it is determined that the discharge trajectory mutual interference deviation is determined, and a trajectory correction strategy is executed;

[0059] If the pulse utilization rate is greater than or equal to 85%, the relative change amount between the joint estimated kinetic parameter vector at the current time and the joint estimated kinetic parameter vector at the preset time is calculated;

[0060] If the numerical value of the relative change amount is greater than the preset change rate threshold, it is determined that the strong model drift deviation is determined, and a trajectory correction strategy is executed;

[0061] If the numerical value of the relative change amount is less than or equal to the preset change rate threshold, it is determined that the weak model drift deviation is determined, and a torque compensation strategy is executed.

[0062] Preferably, the trajectory correction strategy comprises:

[0063] According to the uncertainty of the current kinetic model, a comprehensive correction amount is calculated;

[0064] The mathematical expression of the comprehensive correction amount is:

[0065] ;

[0066] Wherein, is the deviation vector of the actual pose measured by the laser tracker and the theoretical pose; is the deviation value of the pulse utilization rate, is the unit tangent vector of the current trajectory, which ensures the rationality of the correction direction, and are the corresponding weight coefficients, respectively, is a preset reference adjustment step;

[0067] The product of the adaptive trajectory correction gain coefficient and the comprehensive correction vector is superimposed on the theoretical pose of the next control cycle to generate a corrected reference pose.

[0068] Preferably, the torque compensation strategy comprises:

[0069] Based on the joint estimated kinetic parameter vector at the current time, the theoretical torque of the robot is calculated;

[0070] calculating a friction compensation torque based on the friction parameter in the jointly estimated dynamics parameter vector and the current joint motion state;

[0071] calculating a model error compensation torque according to a deviation between the theoretical torque and the actual measured torque, and in combination with an uncertainty of the dynamics model at the current time;

[0072] A mathematical expression of the model error compensation torque is:

[0073] ;

[0074] wherein, the model error compensation torque is, the base integral gain coefficient is, the uncertainty weight coefficient is, the actual measured torque is, the theoretical torque is;

[0075] summing the theoretical torque, the friction compensation torque and the model error compensation torque to generate a final correction torque instruction sent to the motor driver.

[0076] An online trajectory precision correction control system for an electric spark machining robot, comprising a collection module, an identification module and a correction module;

[0077] The collection module is used to collect multi-dimensional data when a machining trajectory is running;

[0078] The identification module is used to construct a unified state vector and obtain an observation matrix according to the multi-dimensional data, and to update dynamics parameters to be identified by using a double-mechanism parameter identification method in combination with the observation matrix;

[0079] The correction module is used to compare a theoretical running trajectory with an actual machining trajectory to obtain a trajectory deviation, to determine a type of the deviation, and to execute a trajectory online correction strategy according to the type of the deviation based on the updated dynamics parameters to be identified.

[0080] The core innovation of the present application is to construct a dynamic evolution model based on double-mechanism parameter identification. Through the first identification mechanism, the forgetting factor adjustment recursive least square method based on the inter-electrode gap change is used to quickly update the inertia parameter when the discharge state fluctuates sharply (such as short circuit), and to suppress noise when the processing is stable, and to accurately track the changes of mass and inertia caused by electrode wear. Combined with the second identification mechanism, the particle swarm optimization algorithm is specially started to identify the nonlinear friction parameters independently in the time window of robot motion stability, effectively solving the problem of difficult accurate modeling of friction characteristics at low speed. This time-sharing and modular joint estimation strategy ensures that the dynamic model can accurately approximate the real physical state of the robot in all directions; for weak model drift, a torque compensation strategy based on model uncertainty is introduced to quickly eliminate small errors through feedforward control and improve dynamic response speed, and for strong model drift or discharge deterioration, a trajectory correction strategy containing a tangent unit vector is executed to ensure that the robot moves strictly along the tangent direction of the machining path during feed adjustment or back arc extinction, avoiding damage to the workpiece profile caused by blind adjustment. This closed-loop control logic combines spatial geometric information and discharge process state, effectively preventing short-circuit arc and ensuring the stability of the inter-electrode gap control. BRIEF DESCRIPTION OF DRAWINGS

[0081] Figure 1 A basic flowchart of an online correction control method for the trajectory accuracy of an electric discharge machining robot is provided for an embodiment of the present application. DETAILED DESCRIPTION

[0082] To make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the specific embodiments of the present application will be described in detail below with reference to the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present application, not all embodiments.

[0083] REFERENCE Figure 1 For an embodiment of the present application, an online correction control method for the trajectory accuracy of an electric discharge machining robot is provided, including the following steps:

[0084] Step S1: Run the machining trajectory and collect multi-dimensional data;

[0085] Step S2: According to the multi-dimensional data, construct a unified state vector and obtain an observation matrix, and update the to-be-identified dynamic parameters by using a double-mechanism parameter identification method combined with the observation matrix;

[0086] Step S3: Compare the theoretical running trajectory with the actual machining trajectory to obtain the trajectory deviation, and determine the type of deviation, and execute the trajectory online correction strategy according to the type of deviation based on the updated to-be-identified dynamic parameters.

[0087] This invention proposes a novel closed-loop control framework that solves the problem of traditional EDM robots struggling to balance trajectory accuracy and discharge stability. By collecting multi-dimensional data and constructing a unified state vector, the robot motion control and EDM process are treated as a coupled system. The dynamic parameters are updated online using a dual-mechanism parameter identification method, which can capture the physical characteristic drift of the robot caused by long-term operation or environmental changes in real time, avoiding the problem of offline calibration parameter failure. The deviation type is determined and a targeted correction strategy is matched.

[0088] Step S1 includes:

[0089] The multi-dimensional data includes preprocessed joint motion data, inter-electrode discharge data, electrode loss data, and joint driving torque.

[0090] Joint motion data includes joint angle, angular velocity, and angular acceleration;

[0091] Inter-electrode discharge data includes inter-electrode voltage, processing current, pulse utilization, and inter-electrode gap;

[0092] Electrode loss data includes the cumulative loss length calculated in real time based on the discharge energy loss model;

[0093] The joint drive torque is the product of the motor's torque constant, reduction ratio, and armature current input to the motor; preprocessing includes noise reduction and outlier removal;

[0094] A unified state vector is constructed based on the preprocessed multi-dimensional data.

[0095] The unified state vector includes the preprocessed robot joint angles, angular velocities, angular accelerations, joint driving torques, inter-electrode gaps, and cumulative electrode loss lengths.

[0096] In one specific embodiment of the present invention, the joint motion data includes the angle, angular velocity, and angular acceleration of each joint read from the robot control system, and all joints of the robot are numbered.

[0097] The inter-electrode gap is the equivalent distance obtained based on the inter-electrode voltage and electrical parameter detection circuit, and is obtained using a linear model, specifically including:

[0098] ;

[0099] for The interpolar interval at any given moment The measured average voltage between electrodes. Set the short-circuit voltage threshold to 0V. The servo gain coefficient is an empirical value representing the number of microns of gap per volt of voltage, and is set to ;

[0100] The pulse utilization rate is the ratio of the number of effective discharge pulses to the total number of pulses;

[0101] The number of effective discharge pulses is the number of pulses that meet the conditions of a voltage drop ≥ 100 V and a current rise ≥ 5 A, and the total number of pulses is the total number of discharge pulses per unit time, which is collected in real time by the discharge state detection circuit;

[0102] In a single pulse cycle, the average voltage in the first 5 microseconds after the breakdown point, at which the voltage begins to drop, is taken as the no-load voltage, and the minimum voltage value in a 2-microsecond time window before and after the time point at which the voltage begins to drop is taken as the breakdown voltage. The difference between the no-load voltage and the breakdown voltage is the voltage drop value;

[0103] The no-load current is recorded as 0 A, and the maximum current value in a 2-microsecond time window before and after the time point at which the current begins to rise is taken as the breakdown current. The difference between the breakdown current and the no-load current is the current rise value;

[0104] The electrode wear data includes the cumulative wear length calculated by integrating the current value over time during the machining process through the electrode material coefficient;

[0105] The mathematical expression of the cumulative wear length is:

[0106] ;

[0107] wherein, is the cumulative wear length of the electrode, is the electrode material coefficient, which reflects the electrode volume wear rate under unit discharge energy. This coefficient is provided by the electrode manufacturer, is the average voltage value across the discharge gap between the tool electrode and the workpiece, is the machining current, which is collected in real time by a Hall current sensor connected in series in the output loop of the spark power supply. This value represents the actual current intensity flowing through the discharge channel between the electrode and the workpiece;

[0108] The tool electrode is the "tool" in the EDM system, i.e., the conductive tool mounted on the end effector of the industrial robot, and the workpiece is the target object to be machined, i.e., the metal part that the robot shapes;

[0109] The torque constant of the motor is a fixed parameter determined by the physical characteristics of the motor, with a unit of The motor is supplied by the motor manufacturer and has a reduction ratio of 100. This reduction ratio was chosen because electrical discharge machining robots typically require high rigidity and a large end-load capacity. This reduction ratio can increase the output torque and suppress cutting vibration.

[0110] Preprocessing includes:

[0111] The joint angle and current data were subjected to a 5th-order Butterworth low-pass filter with a cutoff frequency of 30Hz, and outliers were handled based on the 3σ criterion.

[0112] Based on the preprocessed multi-dimensional data, a unified state vector is constructed. The unified state vector is a 26-dimensional column vector, and its structure is defined as follows: ;

[0113] in, , , Here, represents the joint angle, angular velocity, and angular acceleration, respectively; τ is the joint driving torque calculated through driver current feedback; represents the inter-electrode gap; and L is the cumulative electrode loss length calculated in real time. This vector is updated every control cycle and serves as the input data for the subsequent dual-mechanism parameter identification algorithm. This is a transpose.

[0114] This invention defines in detail the composition and preprocessing of data, laying a data foundation for high-precision identification and control. By introducing electrode loss data based on the discharge energy loss model, it is possible to perceive the microscopic changes in the quality and shape of the tool electrode in real time. This is crucial for compensating for the dynamic model errors caused by electrode wear. A unified state vector including joint motion, driving torque, inter-electrode gap and loss length is constructed, enabling the subsequent algorithm to simultaneously process rigid body dynamics and discharge process characteristics, effectively eliminating noise and outliers, and ensuring the purity of the input data and the robustness of the algorithm.

[0115] The dual-mechanism parameter identification method includes a first identification mechanism and a second identification mechanism. The first identification mechanism includes:

[0116] The observation matrix at the current moment is obtained based on the angle, angular velocity, and angular acceleration of the joints in the unified state vector;

[0117] Calculate the condition number of the observation matrix at the current time, and dynamically update the gain factor based on the condition number;

[0118] The mathematical expression for dynamically updating the gain factor is:

[0119] ;

[0120] in, To dynamically update the gain factor, a preset adjustment constant, a condition number of the observation matrix, an observation matrix; a dynamic update gain factor is used to adjust the overall step size of the subsequent inertial Coriolis force parameter to be identified, so as to prevent divergence of the inertial Coriolis force parameter to be identified due to robot data illness;

[0121] a forgetting factor is calculated based on the inter-pole gap, and the mathematical expression of the forgetting factor is:

[0122] ;

[0123] wherein, a difference between the inter-pole gap at the current moment and the inter-pole gap at the previous moment, a forgetting factor, a reference gap reference value, a sensitivity coefficient; the forgetting factor is used to adjust the retention degree of the historical data, and solve the model lag problem when the processing state is suddenly changed, such as short circuit;

[0124] a gain matrix is calculated according to the observation matrix at the current moment and the forgetting factor, and the mathematical expression of the gain matrix is:

[0125] ;

[0126] wherein, a forgetting factor, an observation matrix, a covariance matrix at the previous moment, a transpose matrix of the observation matrix, determines the direction and basic weight of parameter update, is a normalization standardization factor, used to adjust the size of the gain, so as to prevent excessive update and ensure that the calculated gain will not cause system divergence due to too large covariance matrix;

[0127] the gain matrix is used to calculate the specific correction amount of the inertial Coriolis force parameter to be identified, so as to accurately map the torque residual error to each inertial Coriolis force parameter to be identified;

[0128] in one specific embodiment of the present application, for a 6-DOF robot, the step of obtaining the observation matrix according to the joint angle, angular velocity and angular acceleration comprises:

[0129] obtain the DH parameter table and the URDF file of the robot from the robot manufacturer, the DH parameter table and the URDF file define the length, torsion angle, offset and joint angle of each link of the robot, obtain the dynamic parameters of the robot, including the mass, center of mass position and inertia tensor of each link, the inertia tensor is a 3x3 matrix, used to describe the rotational inertia of the link;

[0130] A rigidBodyTree object is instantiated using robotics software library such as MATLAB's Robotics System Toolbox, each link of the robot is traversed to create corresponding rigidBody and rigidBodyJoint objects, the acquired DH parameter table or URDF file and dynamics parameters are assigned to the corresponding joint and link objects respectively, so as to build a complete, parameterized robot digital twin in the software;

[0131] Symbolic variables are defined for all joint states of the robot using MATLAB's symbolic math toolbox, including a 6-dimensional joint angle vector , an angular velocity vector , and an angular acceleration vector , and an inverse dynamics function such as inverseDynamics is called, taking the digital twin model and the defined symbolic variables as input, which automatically derives the required joint driving torque of the robot according to the Lagrange method, and outputs a 6x1 symbolic torque vector , each element of which is an analytical expression about symbolic variables , , and physical parameters of the robot, including inertia parameters and friction parameters, inertia parameters including the mass, center of mass and inertia tensor of each link, and friction parameters including viscous friction coefficient and Coulomb friction torque;

[0132] A symbolic vector containing all to-be-identified parameters is determined according to the identification requirements;

[0133] The to-be-identified parameters include to-be-identified inertia Coriolis force parameters and to-be-identified friction parameters;

[0134] Each analytical expression in the symbolic torque vector is algebraically reorganized with each element of the symbolic vector containing all to-be-identified physical parameters as a reference, and through symbolic operation, the coefficient term multiplied by in each equation is extracted, which is the element of the observation matrix , and all coefficients are extracted by systematically traversing all i and j, i is 1 to 6, j is 1 to n, and n is the dimension of β;

[0135] All extracted coefficients are combined into a 6xn symbolic observation matrix, each element is a function of , , ;

[0136] The symbolic observation matrix is converted into an independent function file by using the matlab Function of the code generation tool such as MATLAB, and the input of the function is specified as the known joint angle, angular velocity and angular acceleration, and the output is the numerical observation matrix ;

[0137] The mathematical expression of the calculation process of the condition number of the observation matrix is:

[0138] ;

[0139] wherein, is the condition number of the observation matrix, is the observation matrix, represents the matrix norm;

[0140] The mathematical expression of the dynamically updated gain factor is:

[0141] ;

[0142] wherein, is the dynamically updated gain factor, is a preset adjustment constant, is the condition number of the observation matrix, is the observation matrix;

[0143] The value of the preset adjustment constant is 100, and in the field of robot dynamics identification, 100 is usually regarded as the dividing line between good data and bad data, and the condition number reflects the quality of the data. If the robot is stationary, the condition number will be large, and the data is ill-conditioned. At this time, the formula makes the dynamically updated gain factor close to 0, and the algorithm pauses updating to prevent false data from damaging the model. When the robot is moving, the condition number is small, and the dynamically updated gain factor is close to 1, and the algorithm updates at full speed;

[0144] The value of is 0.1, and this coefficient is used to calculate the forgetting factor, which maps the degree of change of the inter-electrode gap to the speed of the algorithm forgetting old data;

[0145] is 50μm, and the core logic of selecting this value is to ensure that the sensitivity of the forgetting factor algorithm is in the effective interval. In typical electric spark machining, the average voltage for maintaining stable discharge is usually between 20V~50V, and according to the formula , when the voltage is 25V, the corresponding physical gap is ;

[0146] therefore, It should be set to the nominal target clearance during processing, i.e., the clearance size under ideal conditions, which is 50μm;

[0147] The condition number of an observation matrix is ​​an indicator of the ill-conditioning of the matrix. The smaller the condition number, the lower the information redundancy of the observation data, and the higher the robustness and accuracy of parameter identification.

[0148] The role of the forgetting factor is that when Δg is large, i.e., a contact short circuit occurs, Making it smaller allows the algorithm to quickly forget old data, while maintaining 0.99 is used for averaging noise.

[0149] The first identification mechanism of this invention greatly improves the stability and adaptability of inertial parameter identification. By calculating the condition number of the observation matrix, the gain factor is dynamically adjusted and updated, which effectively prevents parameter divergence caused by ill-conditioned data when the robot moves smoothly. The forgetting factor is calculated based on the change of the inter-electrode gap, which enables the algorithm to quickly forget old data and adapt to the new state when the discharge state fluctuates violently (short circuit), while retaining more historical data to suppress noise when the processing is smooth.

[0150] The dual-mechanism parameter identification method includes a first identification mechanism and a second identification mechanism. The first identification mechanism further includes:

[0151] Based on the gain matrix, joint driving torque and observation matrix of the previous moment, the recursive least squares method is used to update the inertial Coriolis force parameter vector to be identified online. The inertial Coriolis force parameter vector to be identified consists of the mass of each link, the position of the center of mass of each link and the inertial tensor of each link.

[0152] The mathematical expression for updating the inertial Coriolis force parameter vector to be identified is:

[0153] ;

[0154] in, The vector of inertial Coriolis force parameters to be identified. Here is the gain matrix. For joint driving torque, This is the vector of inertial Coriolis force parameters to be identified at the previous moment. To dynamically update the gain factor, This is the torque error vector, representing the parameters used at the previous moment. The mathematical expression for the difference between the predicted torque and the actual measured torque τ(t), the torque error vector, is as follows:

[0155] ;

[0156] in, a joint driving torque vector measured at a current time, an observation matrix at a current time, an inertial Coriolis force parameter vector to be identified at a previous time;

[0157] The covariance matrix is initialized as a product of a preset constant and a unit matrix, and the dimension of the unit matrix is the same as the column number of the observation matrix;

[0158] The mathematical expression of the covariance matrix is:

[0159] ;

[0160] wherein, is a forgetting factor, is an observation matrix, is a gain matrix, is a covariance matrix at a previous time.

[0161] In specific embodiments of the application, the inertial Coriolis force parameter vector to be identified is composed of the mass of each connecting rod, the center of mass position of each connecting rod and the inertia tensor of each connecting rod;

[0162] The center of mass position of each connecting rod is a three-dimensional coordinate, and the coordinate system is the connecting rod coordinate system, which follows the Denavit-Hartenberg convention;

[0163] The inertia tensor of each connecting rod is originally a 3x3 symmetric matrix containing 6 independent elements, which is stacked into a 6x1 column vector in the process of constructing the inertial Coriolis force parameter vector to be identified, and the inertial Coriolis force parameter vector to be identified is composed of all the masses, center of mass positions and inertia tensors of all connecting rods;

[0164] The inertial Coriolis force parameter vector to be identified at time 0 and the covariance matrix are initialized, and the initialized covariance matrix is a diagonal matrix, and the elements on the diagonal are all 10 6 , which is an empirical value, and the diagonal elements of the covariance matrix represent variance, representing the initial guess of the error range;

[0165] The initialized inertial Coriolis force parameter vector to be identified is obtained from the DH parameter table or the URDF file of the robot of the robot manufacturer;

[0166] At the time when t is equal to 1, the observation matrix at the current time is obtained according to the angle, angular velocity and angular acceleration of the joint;

[0167] The forgetting factor is calculated based on the inter-polar gap, and the mathematical expression of the calculation process is:

[0168]

[0169] wherein, is the difference between the inter-electrode gap at the time t equal to 1 and the time t equal to 0, is the forgetting factor at the time t equal to 1, is equal to 0.1, is equal to 50 pm;

[0170] The gain matrix is calculated according to the observation matrix at the current time and the forgetting factor, and the mathematical expression for calculating the gain matrix is:

[0171]

[0172] wherein, is the forgetting factor, is the observation matrix, is the initialized covariance matrix;

[0173] The mathematical expression for updating the inertial Coriolis force parameter vector to be identified is:

[0174]

[0175] wherein, is the inertial Coriolis force parameter vector to be identified, is the gain matrix, is the joint driving torque, is the initialized inertial Coriolis force parameter vector to be identified;

[0176]

[0177] The calculation of is the covariance matrix at the time t equal to 1, which is used for is the update calculation of the inertial Coriolis force parameter vector to be identified at the time t equal to 2, The mathematical expression for the calculation process of the covariance matrix at the time t equal to 1 is:

[0178]

[0179] wherein, is the forgetting factor, is the observation matrix, is the gain matrix, is the initialized covariance matrix, is the unit matrix;

[0180] After the calculation at the time t equal to 1, the polarity of the inertial Coriolis force parameter to be identified is updated, and this process is repeated at the time t equal to 2, t equal to 3, and so on, so that the inertial Coriolis force parameter to be identified gradually converges to the true value.

[0181] ​​​The application realizes real-time online updating of inertial Coriolis force parameters (mass, center of mass, inertia) by using a recursive least square method, in electric spark machining, the mass and inertia of the tool electrode will continuously decrease with the loss, and the traditional fixed parameter model cannot describe this time-varying process, the method can accurately track the dynamic parameter changes of each connecting rod and electrode, and quantitatively estimate the uncertainty in real time through dynamic updating of the covariance matrix, which ensures that the robot dynamics model can always follow the changes of the actual physical system, and provides accurate mathematical model support for high-precision torque calculation.

[0182] The step S2 further comprises a second identification mechanism, which comprises:

[0183] If the average condition number of the observation matrix continuously exceeds the preset condition number threshold in the preset time window, a parameter identification step of the to-be-identified friction parameter is started;

[0184] The parameter identification step adopts a particle swarm optimization algorithm, and specifically comprises:

[0185] In a preset friction parameter feasible region, a particle swarm is randomly initialized, wherein each particle represents a candidate friction parameter vector, the friction parameter vector is composed of a Coulomb friction coefficient and a viscous friction coefficient, and an initial speed is randomly assigned to each particle;

[0186] The fitness function of each particle is calculated through a fitness function, and the mathematical expression of the fitness function is:

[0187] ;

[0188] Wherein, is a residual torque vector, is a model predicted torque calculated based on the friction parameter vector represented by the current particle, is a weight, is the observation matrix at the k moment; the smaller the fitness value is, the more accurate the friction parameter represented by the particle is, and the condition number penalty term is introduced in the formula to avoid over-optimizing the friction parameter and causing overfitting;

[0189] According to the individual optimal position of each particle in the iteration process and the global optimal position of all particles in the iteration process, the speed and position of all particles are iteratively updated, and the optimized friction parameter vector is obtained by iterating to meet the maximum iteration number.

[0190] In one specific embodiment of the application, if the average condition number of the observation matrix continuously exceeds the preset condition number threshold 500 in the preset time window 500 ms, the parameter identification step is started;

[0191] When the condition number of the observation matrix exceeds 500, the acceleration component in the data is weak, that is, the robot motion is smooth, at this time, the first mechanism inertia identification is based on the logic of The gain is reduced to below 0.16, the gain is basically invalid, and the second mechanism is switched to optimize the friction parameters by taking advantage of the opportunity of smooth robot motion, and the condition number threshold is 500 to realize the relay of the two algorithms;

[0192] The preset friction parameter feasible region is composed of the value range of the Coulomb friction coefficient and the viscous friction coefficient, the value range of the Coulomb friction coefficient is between 5% and 25% of the rated torque, and the viscous friction coefficient is between 1% and 10% of the Coulomb friction coefficient, the rated torque is provided by the manufacturer, if the rated torque of the joint motor is 100 Nm, the value range of the Coulomb friction coefficient is [5, 25] units Nm, and the value range of the viscous friction coefficient is [0.05, 2.5];

[0193] The velocity of the particle determines the change amount of the Coulomb friction coefficient and the viscous friction coefficient of the particle in the next iteration, so the velocity of the particle includes the change velocity of the Coulomb friction coefficient and the change velocity of the viscous friction coefficient, and the range of the initial velocity is set to ±10% corresponding to the width of the feasible region space;

[0194] For example, the value range of the initial velocity of the change velocity of the Coulomb friction coefficient is [-2, 2];

[0195] The value range of the initial velocity of the change velocity of the viscous friction coefficient is [-0.245, 0.245];

[0196] ;

[0197] is the difference between the measured torque at t and ;

[0198] The mathematical expression of is:

[0199] ;

[0200] Wherein, is the friction model predicted torque, is the current joint angular velocity, is the Coulomb friction coefficient vector, is the viscous friction coefficient vector, from the current position of the particle, is the sign function, 1 when the velocity is positive, and -1 when the velocity is negative;

[0201] According to the individual optimal position and the global optimal position of the particle, the particle velocity and position are updated;

[0202] The individual optimal position refers to the position with the lowest fitness value, i.e., the minimum error, experienced by each particle from the beginning to the current iteration, and the position represents the friction parameter vector as a coordinate in a parameter space, which is a two-dimensional coordinate system, with the x-axis being the Coulomb friction coefficient and the y-axis being the viscous friction coefficient;

[0203] The global optimal position is the position with the lowest fitness value found by all particles in the entire particle swarm in all historical iterations. In each iteration, after the individual optimal position of each particle is determined, the individual optimal positions of all particles are compared to find the position corresponding to the particle with the lowest fitness value as the global optimal position, which is the target pursued by the entire swarm;

[0204] The mathematical expression for updating the velocity is:

[0205]

[0206] wherein w is the inertia weight, c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. w is the inertia weight, c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position.

[0207] c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position.

[0208] c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. c1 and c2 are learning factors, rand1 and rand2 are random numbers in the range of [0, 1], and x is the individual optimal position, g is the global optimal position. Is a random variable uniformly distributed between 0 and 1, for the uniform distribution on the interval [0,1], the center of the probability density, that is, the average value, is located in the middle of the interval, that is, 0.5, the expected value of c1*r1 and c2*r2 is 1, which means that the particle will approach the individual optimal position with a step size of one step, and at the same time, it will approach the global optimal position with a step size of one step, which makes the particle keep balance between individual experience and group learning, and the particle will And The value set to be equal means that the algorithm equally values the experience of the particle itself and the experience of the whole group;

[0209] The mathematical expression of position updating is as follows:

[0210]

[0211] Wherein, Pi(t) is the position vector of particle i at the current iteration t, that is, the friction parameter vector, Pi(t-1) is the position vector of particle i at the last iteration t-1, Vi(t) is the velocity vector of particle i at the current iteration t.

[0212] The second identification mechanism of the application uses a time window in which the robot motion is stable (the condition number of the observation matrix is high), and specially starts the friction parameter identification based on the particle swarm optimization algorithm, avoiding the case that the inertia force covers the friction force characteristics when moving at high speed, the PSO algorithm can more accurately fit the nonlinear coulomb friction and viscous friction coefficient through global optimization, and the condition number is introduced as a penalty term of the fitness function, which effectively prevents overfitting and significantly improves the accuracy of joint friction compensation when processing at low speed.

[0213] Step S2 further includes:

[0214] A total dynamics parameter vector is constructed, which is composed of inertia coriolis force parameters and friction parameters;

[0215] The inertia coriolis force parameter part is an initialized inertia coriolis force parameter vector to be identified;

[0216] The friction parameter part is an initialized friction parameter to be identified;

[0217] In each control period, the first identification mechanism is executed to obtain the updated inertia coriolis force parameter to be identified, and the inertia coriolis force parameter to be identified in the total dynamics parameter vector is replaced by the inertia coriolis force parameter to be identified, and the friction parameter part remains unchanged;

[0218] ​When the parameter identification step of the second identification mechanism is triggered, the second identification mechanism is executed to obtain an optimized friction parameter vector, and the optimized friction parameter vector is used to overwrite the corresponding friction parameter part of the total dynamic parameter vector, while the inertial Coriolis force parameter part remains unchanged, thereby obtaining a joint estimated dynamic parameter vector.

[0219] In one specific embodiment of the present application, one control cycle refers to the time unit required for the robot control system to perform a complete process from data acquisition to parameter identification to execution of the correction instruction;

[0220] For a 6-DOF electric spark machining robot, the initial inertial Coriolis force parameter vector and the friction parameters including the Coulomb friction coefficient and the viscous friction coefficient are obtained according to the default dynamic parameter model provided by the simulation or control software of the robot manufacturer,

[0221] The initial inertial Coriolis force parameter vector and the initial friction parameter vector are spliced to obtain a total dynamic parameter vector;

[0222] After the robot starts to perform the machining task, the first identification mechanism is executed to update the inertial Coriolis force parameter vector, and the new inertial Coriolis force parameter vector is used to overwrite the elements of the inertial Coriolis force parameter in the total dynamic parameter vector, while the friction parameter part remains unchanged;

[0223] When it is monitored that the cumulative tool electrode loss length has exceeded the preset threshold value of 0.5 mm, the second identification mechanism is executed to initialize a particle swarm in a preset feasible region, and to use the multi-dimensional data collected in the past 5 seconds to perform 50 iterations of optimization to obtain an optimized friction parameter vector;

[0224] The set of optimized friction parameter vectors is used to overwrite the last two elements in the total dynamic parameter vector, and at this time, the inertial Coriolis force parameter part remains the value at the moment before triggering.

[0225] Step S2 also strategically realizes complementary updating of the inertial parameters and the friction parameters, by respectively updating different parts in the total dynamic vector under different control cycles and working conditions, parameter drift or mutual interference caused by a single algorithm attempting to solve all problems at the same time is avoided, and this time-sharing and modular updating mechanism ensures that the mass inertia part and the friction part can both be maintained in the latest and most accurate state, thereby providing a high-credibility full-parameter model for subsequent trajectory correction and torque compensation.

[0226] Step S3 includes:

[0227] The difference between the theoretical end effector pose and the actual end effector pose is taken as the trajectory deviation, and it is determined whether the trajectory deviation exceeds the preset dynamic trajectory reference value, and if the trajectory deviation exceeds the preset dynamic trajectory reference interval threshold value, deviation type determination is performed.

[0228] when the trajectory deviation exceeds a preset dynamic trajectory reference threshold value;

[0229] if the pulse utilization rate is less than 85%, it is determined that the discharge trajectory mutual interference deviation is large, and a trajectory correction strategy is executed;

[0230] if the pulse utilization rate is greater than or equal to 85%, the relative change amount between the joint estimation kinetic parameter vector at the current time and the joint estimation kinetic parameter vector at the preset time is calculated;

[0231] if the numerical value of the relative change amount is greater than a preset change rate threshold value, it is determined that the strong model drift deviation is large, and a trajectory correction strategy is executed;

[0232] if the numerical value of the relative change amount is less than or equal to the preset change rate threshold value, it is determined that the weak model drift deviation is small, and a torque compensation strategy is executed.

[0233] In one specific embodiment of the present application, a laser tracker is used to obtain the actual end effector pose, including position coordinates and a quaternion, the quaternion is converted into a rotation matrix, which is used for subsequent calculation of trajectory deviation;

[0234] The trajectory deviation includes a position deviation and an attitude deviation, the position deviation is the Euclidean distance between the three-dimensional coordinates of the theoretical end effector position and the three-dimensional coordinates of the actual position;

[0235] The attitude deviation is obtained in the following manner:

[0236] The product of the actual rotation matrix and the transpose matrix of the theoretical attitude matrix is taken as an error rotation matrix, and the structure of the error rotation matrix is as follows:

[0237]

[0238] wherein, represents the element value of the first row, the first column of the matrix.

[0239] The error rotation matrix is converted into a rotation angle and a rotation axis, and the product of the rotation angle and the rotation axis is taken as the attitude deviation;

[0240] The mathematical expression for calculating the rotation angle is:

[0241]

[0242] wherein, is the rotation angle, is the error rotation matrix, is a trace;

[0243] The components of the rotation axis are calculated and normalized to obtain the rotation axis, specifically including:​​​​

[0244] Rotation axis is a unit vector , which indicates the direction of rotation, based on the inverse transformation of the Rodriguez formula, the anti-symmetric part of the matrix is used to extract the axis component.

[0245] The complete calculation formula is as follows:

[0246] ;

[0247] Wherein, the meaning is to extract the component of rotation around the X axis by using the anti-symmetric property of the matrix;

[0248] The meaning of is to extract the component of rotation around the Y axis;

[0249] The meaning of is to extract the component of rotation around the Z axis

[0250] The normalization factor ensures that the calculated vector length is 1;

[0251] The attitude deviation is a rotation vector, the size is the angle, and the direction is the axis;

[0252] ;

[0253] When the trajectory deviation exceeds the preset dynamic trajectory reference threshold, including when the position deviation is greater than the preset reference threshold 0.2mm or the modulus of the attitude deviation is greater than the preset reference threshold 1°, the deviation type determination is performed, and the dynamic trajectory reference threshold is determined according to the requirement of machining accuracy;

[0254] When the trajectory deviation exceeds the preset dynamic trajectory reference threshold;

[0255] If the pulse utilization rate is less than 85%, it is determined that the discharge trajectory mutual interference deviation is determined, and the trajectory correction strategy is executed; 85% is an empirical value;

[0256] If the pulse utilization rate is greater than or equal to 85%, the relative change quantity between the joint estimation kinetic parameter vector at the current time and the joint estimation kinetic parameter vector at the preset time is calculated;

[0257] If the numerical value of the relative change quantity is greater than the preset change rate threshold 10%, it is determined that the strong model drift deviation is determined, and the trajectory correction strategy is executed;

[0258] If the numerical value of the relative change quantity is less than or equal to the preset change rate threshold 10%, it is determined that the weak model drift deviation is determined, and the torque compensation strategy is executed.

[0259] ​The application sets a dynamic trajectory reference threshold and a pulse utilization rate threshold, so that the system can accurately identify whether the deviation is caused by discharge state deterioration (discharge trajectory mutual interference), model parameter drastic change (strong model drift) or slight disturbance (weak model drift), the refined classification determination logic is the premise of realizing hierarchical control, and the machining efficiency is maximized under the premise of ensuring machining safety.

[0260] The trajectory correction strategy includes:

[0261] According to the uncertainty of the current moment dynamics model, the comprehensive correction amount is calculated;

[0262] The mathematical expression of the comprehensive correction amount is:

[0263] ;

[0264] Wherein, is the deviation vector of the actual pose measured by the laser tracker and the theoretical pose; is the deviation value of the pulse utilization rate, is the unit tangent vector of the current trajectory, which ensures the rationality of the correction direction, and are the corresponding weight coefficients, is a preset reference adjustment step 1mm;

[0265] The product of the adaptive trajectory correction gain coefficient and the comprehensive correction vector is added to the theoretical pose of the next control period to generate the corrected reference pose.

[0266] In one embodiment of the application, when the recursive least squares method is used to update the dynamic parameters in step S2, the covariance matrix is also updated synchronously. The trace of the covariance matrix, i.e. the sum of the diagonal elements, is used as an uncertainty index to quantify the confidence of the system in the current dynamic parameter estimation. The larger the value of the uncertainty index, the higher the uncertainty of the dynamic model.

[0267] The adaptive trajectory correction gain coefficient is calculated according to the dynamic model uncertainty index, and the mathematical expression is:

[0268] ;

[0269] Wherein, is a preset basic trajectory correction gain coefficient 1.0, representing the reference correction strength when the model is completely determined, is a preset uncertainty weight coefficient 50, used to adjust the influence degree of model uncertainty on the gain, is the trace of the covariance matrix as the uncertainty index;

[0270] The comprehensive correction vector is calculated according to the following formula:

[0271] ;

[0272] Wherein, is the deviation vector of the actual pose measured by the laser tracker and the theoretical pose; is the normalized deviation value of the pulse utilization rate, is the unit tangent vector of the current trajectory, which ensures the rationality of the correction direction, and are the corresponding weight coefficients, respectively;

[0273] The deviation value of the pulse utilization rate is the difference between the pulse utilization rate and the pulse utilization rate target value, and the pulse utilization rate target value is a preset value of 85%, because if the pulse utilization rate is greater than or equal to 85%, the correction is performed;

[0274] is the position correction weight, which is 0.6, meaning that 60% of the current error is corrected in each cycle, and the remaining error is left for correction in the next cycle, which ensures the tracking speed and smooths the measurement noise;

[0275] is the process correction weight, which is 0.5, which gives a clear physical meaning, i.e. for every 100% of utilization deviation, there is a 0.5mm feed adjustment;

[0276] is a vector composed of the position deviation vector and the attitude deviation vector, which is a 6-dimensional vector;

[0277] The position deviation vector is the vector difference between the column vector constructed by the three-dimensional coordinates of the theoretical end effector position and the column vector constructed by the three-dimensional coordinates of the actual position ( );

[0278] The attitude deviation vector is ( );

[0279] is the tangent direction of the current motion path of the robot tool end, which is to ensure that when the robot needs to back off or slow down, this action must be along the just-cut trajectory, and cannot run perpendicular to the trajectory, otherwise the workpiece shape will be cut badly. The calculation method is:

[0280] ;

[0281] That is, the current instantaneous velocity vector divided by the size of the velocity;

[0282] is the linear velocity vector of the current end effector , the linear velocity vector is calculated by mapping the joint space motion to the cartesian space through the robot Jacobian matrix, the real-time calculation of the current rigidBodyTree Jacobian matrix is realized by using the built digital twin, the Jacobian matrix describes the linear relationship between joint velocity and end cartesian velocity,

[0283] ;

[0284] is a matrix (for 6-axis robots).

[0285] is the joint angular velocity vector of .

[0286] The calculation result is the linear velocity vector of the end effector in space ;

[0287] is the current feed speed, which is the length (norm) of the linear velocity vector;

[0288] In the calculation process of the comprehensive correction vector, the linear velocity vector of , so contains three components ( ), which is inconsistent with the dimension of , and cannot be calculated, so the dimension of needs to be expanded to for subsequent calculation;

[0289] The corrected reference pose of the next control cycle is generated, and its calculation formula is:

[0290] ;

[0291] The corrected pose will replace the original theoretical pose as the tracking target of the controller in the next cycle, realizing dynamic and forward-looking correction of the future trajectory;

[0292] The original theoretical pose comes from the machining file (G code), which contains position and attitude, the position is , describing the coordinates of the tool tip in space, and the attitude is , describing the angle of the tool electrode, usually represented by Euler angles or rotation vectors, the original theoretical pose is .

[0293] The present application ensures the feed adjustment based on the pulse utilization rate by introducing the tangential unit vector, strictly along the machining path, effectively preventing the damage to the workpiece shape caused by blind tool withdrawal, at the same time, combining the model uncertainty to calculate the adaptive gain coefficient, so that the correction amount contains both the geometric error correction of the laser tracker feedback and the process displacement required to maintain the discharge gap, realizing the forward-looking dynamic compensation of the reference pose at the next moment, and significantly improving the geometric precision and discharge stability of the machining.

[0294] The torque compensation strategy includes:

[0295] Based on the joint estimation of the dynamic parameter vector at the current moment, the theoretical torque of the robot is calculated;

[0296] Based on the friction parameter in the joint estimation of the dynamic parameter vector and the current joint motion state, the friction compensation torque is calculated;

[0297] According to the deviation of the theoretical torque and the actual measured torque, and combined with the uncertainty of the dynamic model at the current moment, the model error compensation torque is calculated;

[0298] The mathematical expression of the model error compensation torque is:

[0299] ;

[0300] Wherein, is the model error compensation torque, is the basic integral gain coefficient, is the uncertainty weight coefficient, is the actual measured torque, is the theoretical torque; is the basic integral gain coefficient, and the physical unit is inverse second , the numerical value is The coefficient is to restore the torque-time accumulation (impulse moment) generated by integration to an instantaneous torque compensation value;

[0301] The theoretical torque, friction compensation torque and model error compensation torque are summed up to generate the final correction torque instruction sent to the motor driver.

[0302] In one embodiment of the present application, the updated joint estimation dynamic parameter vector of step S2 is called, and the theoretical torque vector required by the robot motion is derived based on the expected joint motion state through the inverse dynamics equation, the expected joint motion state is the state that the robot should be in at the next moment planned by the robot controller according to the machining G code, specifically including three vectors, which are the expected angle, the expected angular velocity and the expected angular acceleration;

[0303] The newly recognized Coulomb friction and viscous friction parameters are extracted from the jointly estimated dynamics parameter vector, and the current joint angular velocity is combined to accurately calculate the friction compensation torque vector. The calculation of the friction compensation torque vector is used for independent calculation and compensation of the friction force in the joint bearing and the internal friction of the speed reducer. The calculation formula is:

[0304] ;

[0305] The Coulomb friction coefficient vector is obtained by the second mechanism recognition of step S2, is a sign function, and the output is +1 when the speed is positive and -1 when the speed is negative, The viscous friction coefficient vector is obtained by the second mechanism recognition of step S2, is the current joint angular velocity;

[0306] The mathematical expression of the model error compensation torque is:

[0307] ;

[0308] wherein, is the model error compensation force;

[0309] is the basic integral gain coefficient, which is 1.0, and 1.0 is a safe starting point, meaning that 100% of the cumulative error is eliminated per second;

[0310] is the uncertainty weight coefficient, which is a preset value, and the value range is 10-100. The selection of the uncertainty weight coefficient depends on the order of magnitude of the trace of the covariance matrix;

[0311] is the actual measured torque, which represents the actual output value of the motor driver through the current loop feedback, is the theoretical torque calculated based on the current dynamics parameters. The current joint state of the joint angle, joint angular velocity and joint angular acceleration is read from the robot controller. The mass of each link, the position of the mass center of each link, the inertia tensor of each link and the friction coefficient are extracted from the updated parameters in step S2, and substituted into the standard equation of inverse dynamics. The mathematical expression is automatically executed by an algorithm library such as Newton-Euler recursive algorithm, which is:

[0312] ;

[0313] wherein, is the inertial force term, is the Coriolis force and centrifugal force term, is the gravity term, is the friction force term, using the Coulomb friction coefficient vector recognized by the second mechanism of S2 Substitute calculation is carried out;

[0314] The theoretical torque vector, the friction compensation torque vector and the error compensation torque vector are summed to generate a final correction torque instruction sent to the motor driver;

[0315] The instruction combines the dynamic model prediction, friction characteristic compensation and adaptive error suppression strategy, and can significantly improve the trajectory tracking accuracy of the robot under complex working conditions.

[0316] The correction torque instruction is in the form of a digital integer, which is sent to the driver through a specific object dictionary such as Torque Offset 0x60B2 of the EtherCAT real-time bus, and the driver converts it into a current instruction, and finally controls the power electronic elements to generate physical torque through PWM.

[0317] The correction strategy of the application not only compensates the theoretical torque and friction torque, but also introduces a model error compensation term based on the trace of the covariance matrix (uncertainty), which means that when the model is less certain, the compensation weight of the controller will be adjusted accordingly, thereby suppressing errors while ensuring the stability of the control system. This compensation method directly acting on the motor current loop reacts faster than simple trajectory planning and can effectively eliminate small trajectory errors caused by friction or slight parameter fluctuations.

[0318] An electric spark machining robot trajectory precision online correction control system, comprising a collection module, an identification module and a correction module;

[0319] The collection module is used to collect multi-dimensional data when running the machining trajectory;

[0320] The identification module is used to construct a unified state vector according to the multi-dimensional data, and obtain an observation matrix, and a double mechanism parameter identification method is used to update the to-be-identified dynamic parameters combined with the observation matrix;

[0321] The correction module is used to compare the theoretical running trajectory with the actual machining trajectory to obtain a trajectory deviation, and determine the type of deviation, and based on the updated to-be-identified dynamic parameters, execute a trajectory online correction strategy according to the type of deviation.

[0322] The core innovation of the application is to construct a dynamic evolution model based on double mechanism parameter identification. Through the first identification mechanism, the forgetting factor adjustment recursive least square method based on the inter-electrode gap change is used, so that the algorithm can quickly update the inertia parameter when the discharge state fluctuates sharply (such as short circuit), and can suppress noise when the processing is stable, and accurately track the changes of mass and inertia caused by electrode wear. Combined with the second identification mechanism, the particle swarm optimization algorithm is specially started to identify the nonlinear friction parameters independently in the time window of robot motion stability, effectively solving the problem of difficult accurate modeling of friction characteristics at low speed processing. This time-sharing and modular joint estimation strategy ensures that the dynamic model can comprehensively and accurately approximate the real physical state of the robot in real time. For weak model drift, a torque compensation strategy based on model uncertainty is introduced to quickly eliminate small errors through feedforward control and improve dynamic response speed. For strong model drift or discharge deterioration, a trajectory correction strategy containing a tangent unit vector is executed to ensure that the robot moves strictly along the tangent direction of the processing path during feed adjustment or back arc elimination, avoiding the damage of workpiece contour caused by blind adjustment. This closed-loop control logic combining spatial geometric information and discharge process state effectively prevents short-circuit arc and ensures the stability of the inter-electrode gap control.

[0323] Those skilled in the art will appreciate that embodiments of the application can be provided as methods, systems or computer program products. Accordingly, the application can be embodied in the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the application can be embodied in the form of a computer program product embodied on one or more computer-usable storage media having computer-usable program code embodied thereon. The storage media can be any type of volatile or non-volatile storage device, or a combination thereof, such as static random access memory (SRAM), electrically erasable programmable read-only memory (EEPROM), erasable programmable read-only memory (EPROM), programmable read-only memory (PROM), read-only memory (ROM), magnetic storage, flash memory, magnetic disk or optical disk. These computer program instructions can also be stored in a computer readable storage medium that can guide a computer or other programmable data processing device to work in a specific way, so that the instructions stored in the computer readable storage medium produce an instruction device including the instructions, which implement the functions described in the flowcharts Figure 1one or more processes and / or blocks Figure 1 the function specified in the one or more blocks.

[0324] It should be noted that the above-mentioned embodiments are only used to illustrate the technical solutions of the present application, not to limit the present application. Although the present application has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present application can be modified or equivalent replaced without departing from the spirit and scope of the technical solutions of the present application, which should be covered in the scope of the claims of the present application.

Claims

1. An electric discharge machining robot trajectory accuracy online correction control method, characterized by, The method comprises the following steps: Step S1: running a machining trajectory and collecting multi-dimensional data; Step S2: constructing a unified state vector according to the multi-dimensional data, obtaining an observation matrix, and updating to-be-identified dynamic parameters by using a double-mechanism parameter identification method combined with the observation matrix; Step S3: comparing a theoretical running trajectory with an actual machining trajectory to obtain a trajectory deviation, determining a type of the deviation, and executing a trajectory online correction strategy according to the type of the deviation based on the updated to-be-identified dynamic parameters; The double-mechanism parameter identification method comprises a first identification mechanism and a second identification mechanism, and the first identification mechanism comprises: An observation matrix at a current time is obtained according to angles, angular velocities and angular accelerations of joints in the unified state vector; A condition number of the observation matrix at the current time is calculated, and a dynamic update gain factor is calculated according to the condition number; A mathematical expression of the dynamic update gain factor is as follows: ； wherein, is a dynamic update gain factor, is a preset adjustment constant, is a condition number of the observation matrix, is an observation matrix; A forgetting factor is calculated based on an inter-electrode gap, and a mathematical expression of the forgetting factor is as follows: ； wherein, is a difference between the current time and the last time of the interpolar gap, is a forgetting factor, is a reference gap reference value, is a sensitivity coefficient; A gain matrix is calculated according to the observation matrix at the current time and the forgetting factor, and a mathematical expression of the gain matrix is as follows: ； wherein is a forgetting factor, is an observation matrix, is a covariance matrix of the previous time instant, is a transpose matrix of the observation matrix; The second identification mechanism comprises: If an average condition number of the observation matrix continuously exceeds a preset condition number threshold within a preset time window, a parameter identification step of to-be-identified friction parameters is started; The parameter identification step adopts a particle swarm optimization algorithm, and specifically comprises: In a preset friction parameter feasible region, a particle swarm is randomly initialized, wherein each particle represents a candidate friction parameter vector composed of a Coulomb friction coefficient and a viscous friction coefficient, and each particle is randomly assigned an initial speed; An adaptability function of each particle is calculated by an adaptability function, and a mathematical expression of the adaptability function is as follows: ； wherein, is a residual torque vector, is a model predicted torque calculated based on a friction parameter vector represented by the current particle, is a weight, is an observation matrix at time k; According to an individual optimal position of each particle in an iteration process and a global optimal position of all particles in the iteration process, the speed and position of all particles are iteratively updated, and the iteration is performed until a maximum iteration number is satisfied, and an optimized friction parameter vector is obtained; The step S2 further comprises: A total dynamic parameter vector is constructed, and the vector is composed of inertial Coriolis force parameters and friction parameters; The inertial Coriolis force parameter part is an initialized to-be-identified inertial Coriolis force parameter vector; The friction parameter part is an initialized to-be-identified friction parameter; In each control cycle, the first identification mechanism is executed to obtain updated to-be-identified inertial Coriolis force parameters, and the to-be-identified inertial Coriolis force parameter vector is used to replace the corresponding to-be-identified inertial Coriolis force parameter part of the total dynamic parameter vector, and the friction parameter part remains unchanged; When the parameter identification step of the second identification mechanism is triggered, the second identification mechanism is executed to obtain an optimized friction parameter vector, and the optimized friction parameter vector is used to replace the corresponding friction parameter part of the total dynamic parameter vector, and the inertial Coriolis force parameter part remains unchanged, thereby obtaining a joint estimated dynamic parameter vector.

2. The electric discharge machining robot trajectory accuracy on-line correction control method according to claim 1, characterized by, The step S1 comprises: The multi-dimensional data is preprocessed joint motion data, inter-electrode discharge data, electrode loss data and joint driving torque; The joint motion data is angles, angular velocities and angular accelerations of joints. The inter-electrode discharge data includes inter-electrode voltage, machining current, pulse utilization rate and inter-electrode gap; The electrode loss data includes cumulative loss length calculated in real time based on a discharge energy loss model; The joint driving torque is the product of torque constant of the motor, speed reduction ratio and armature current input to the motor; The preprocessing includes denoising and outlier rejection, and a unified state vector is constructed based on the preprocessed multi-dimensional data, elements in the unified state vector including angles, angular velocities, angular accelerations of joints of the robot after preprocessing, joint driving torque, inter-electrode gap and cumulative loss length of the electrode.

3. The control method of the electric discharge machining robot trajectory precision on-line correction according to claim 1, wherein The double-mechanism parameter identification method includes a first identification mechanism and a second identification mechanism, and the first identification mechanism further includes: Based on the gain matrix, joint driving torque and observation matrix of the previous moment, the inertial Coriolis force parameter vector to be identified is updated online by using the recursive least square method, and the inertial Coriolis force parameter vector to be identified is composed of the mass of each connecting rod, the mass center position of each connecting rod and the inertia tensor of each connecting rod; The mathematical expression for updating the inertial Coriolis force parameter vector to be identified is: ； wherein is a vector of inertial Coriolis force parameters to be identified, is a gain matrix, is a joint driving torque, is a vector of inertial Coriolis force parameters to be identified at the previous time, is a dynamic update gain factor, is a torque error vector, the mathematical expression of which is: ； wherein, is a joint driving torque vector measured at the current time instant, is an observation matrix at the current time instant, is an inertial Coriolis force parameter vector to be identified at the previous time instant; The covariance matrix is initialized as the product of a preset constant and a unit matrix; The mathematical expression of the covariance matrix is: ； wherein is a forgetting factor, is an observation matrix, is a gain matrix, is a covariance matrix of the previous time instant.

4. The electric discharge machining robot trajectory accuracy on-line correction control method according to claim 3, characterized by, The step S3 includes: The difference between the theoretical end effector pose and the actual end effector pose is taken as the trajectory deviation, and it is judged whether the trajectory deviation exceeds the preset dynamic trajectory reference value, if the trajectory deviation exceeds the preset dynamic trajectory reference interval threshold, the deviation type is determined; When the trajectory deviation exceeds the preset dynamic trajectory reference threshold; If the pulse utilization rate is less than 85%, it is determined that the discharge trajectory mutual interference deviation occurs, and the trajectory correction strategy is executed; If the pulse utilization rate is greater than or equal to 85%, the relative change quantity between the joint estimation dynamics parameter vector at the current moment and the joint estimation dynamics parameter vector at the preset moment is calculated; If the numerical value of the relative change quantity is greater than the preset change rate threshold, it is determined that the strong model drift deviation occurs, and the trajectory correction strategy is executed; If the numerical value of the relative change quantity is less than or equal to the preset change rate threshold, it is determined that the weak model drift deviation occurs, and the torque compensation strategy is executed.

5. The electrical discharge machining robot trajectory accuracy on-line correction control method according to claim 4, characterized by, The trajectory correction strategy includes: According to the uncertainty of the dynamics model at the current moment, an adaptive trajectory correction gain coefficient and a comprehensive correction quantity are calculated; The mathematical expression of the adaptive trajectory correction gain coefficient is: ； wherein, is a preset base trajectory correction gain coefficient, is a preset uncertainty weight coefficient for adjusting the influence degree of model uncertainty on the gain, is an uncertainty index; The mathematical expression of the comprehensive correction quantity is: ； wherein, is a deviation vector of the actual pose measured by the laser tracker and the theoretical pose; is a deviation value of the pulse utilization rate, is a unit tangent vector of the current trajectory, ensuring the rationality of the correction direction, and are corresponding weight coefficients, respectively, is a preset reference adjustment step. The product of the adaptive trajectory correction gain coefficient and the comprehensive correction quantity is superimposed on the theoretical pose of the next control cycle to generate the corrected reference pose.

6. The electrical discharge machining robot trajectory accuracy on-line correction control method according to claim 5, wherein The torque compensation strategy includes: Based on the joint estimation dynamics parameter vector at the current moment, the theoretical torque of the robot is calculated; Based on the friction parameter in the joint estimation dynamics parameter vector and the current joint motion state, the friction compensation torque is calculated; According to the deviation between the theoretical torque and the actual measured torque, and combined with the uncertainty of the dynamics model at the current moment, the model error compensation torque is calculated; The mathematical expression of the model error compensation torque is: ； wherein, is a model error compensation torque, is a base integral gain coefficient, is an uncertainty weight coefficient, is an actual measured torque, is a theoretical torque; Summing the theoretical torque, the friction compensation torque and the model error compensation torque generates a final correction torque instruction sent to the motor driver.

7. An electric discharge machining robot trajectory accuracy online correction control system, characterized by, The method comprises a collection module, an identification module and a correction module; The collection module is configured to collect multi-dimensional data when the machining trajectory is running. The identification module is configured to construct a unified state vector according to the multi-dimensional data, obtain an observation matrix, and update the to-be-identified dynamic parameters by using a double-mechanism parameter identification method combined with the observation matrix. The correction module is configured to compare the theoretical trajectory with the actual machining trajectory to obtain a trajectory deviation, determine the type of the deviation, and execute a trajectory online correction strategy according to the type of the deviation based on the updated to-be-identified dynamic parameters. The double-mechanism parameter identification method comprises a first identification mechanism and a second identification mechanism. The first identification mechanism comprises the following steps: An observation matrix at the current time is obtained according to the angle, angular velocity and angular acceleration of the joint in the unified state vector. A condition number of the observation matrix at the current time is calculated, and a dynamic update gain factor is calculated according to the condition number. ； wherein, is a dynamic update gain factor, is a preset adjustment constant, is a condition number of the observation matrix, is an observation matrix; The mathematical expression of the dynamic update gain factor is as follows: ； wherein, is the difference between the current time and the last time of the inter-pole gap, is a forgetting factor, is a reference gap reference value, is a sensitivity coefficient; A forgetting factor is calculated based on the inter-pole clearance, and the mathematical expression of the forgetting factor is as follows: ； wherein, is a forgetting factor, is an observation matrix, is a covariance matrix of the previous time, is a transpose matrix of the observation matrix; A gain matrix is calculated according to the observation matrix at the current time and the forgetting factor, and the mathematical expression of the gain matrix is as follows: The second identification mechanism comprises the following steps: If the average condition number of the observation matrix continuously exceeds a preset condition number threshold within a preset time window, a parameter identification step for the to-be-identified friction parameters is started. The parameter identification step adopts a particle swarm optimization algorithm, and specifically comprises the following steps: A particle swarm is randomly initialized in a preset friction parameter feasible region, wherein each particle represents a candidate friction parameter vector composed of a Coulomb friction coefficient and a viscous friction coefficient, and each particle is randomly assigned an initial speed. ； wherein, is a residual torque vector, is a model predicted torque calculated based on a friction parameter vector represented by the current particle, is a weight, is an observation matrix at time k; The fitness function of each particle is calculated by using a fitness function, and the mathematical expression of the fitness function is as follows: The speed and position of each particle are iteratively updated according to the individual optimal position of each particle in the iteration process and the global optimal position of all particles in the iteration process, and the iteration is performed until a maximum iteration number is satisfied, and an optimized friction parameter vector is obtained. The identification module is further configured to: A total dynamic parameter vector is constructed, which is composed of inertial Coriolis force parameters and friction parameters. The inertial Coriolis force parameter part is an initialized to-be-identified inertial Coriolis force parameter vector. The friction parameter part is an initialized to-be-identified friction parameter. In each control cycle, the first identification mechanism is executed to obtain updated to-be-identified inertial Coriolis force parameters, and the to-be-identified inertial Coriolis force parameter part of the total dynamic parameter vector is replaced by the updated to-be-identified inertial Coriolis force parameter vector, and the friction parameter part remains unchanged. When the parameter identification step of the second identification mechanism is triggered, the second identification mechanism is executed to obtain an optimized friction parameter vector, and the corresponding friction parameter part of the total dynamic parameter vector is replaced by the optimized friction parameter vector, and the inertial Coriolis force parameter part remains unchanged, thereby obtaining a joint estimated dynamic parameter vector.

Images