Robot Calibration Method Based on Adaptive Momentum LM Cascade B-Spline Interpolation Particle Swarm

Abstract

The present invention relates to the technical field of robot calibration, and discloses a robot calibration method based on adaptive momentum LM cascaded B-spline interpolation particle swarm, including the steps of: constructing a kinematic error model based on the LM model, constructing a loss function to optimize the D-H parameters of the kinematic error model; constructing a least squares objective function based on the kinematic error model, and obtaining a sub-optimal solution of the D-H parameters by using the AMLM algorithm; taking the sub-optimal solution as the central value of the initial population of BIPSO, generating a candidate solution set in the BIPSO algorithm, and obtaining the optimal solution of the robot D-H parameter error through the BIPSO algorithm. The present invention introduces a momentum term into the traditional LM algorithm to suppress the overshoot oscillation phenomenon, incorporates the B-spline interpolation technology into the classical particle swarm optimization algorithm PSO, cascades these two improved algorithms, and realizes a more accurate estimation of the geometric parameters of the robot.

Description

Technical Field

[0001] The present invention relates to the technical field of robot calibration, and particularly relates to a robot calibration method based on an adaptive momentum LM cascaded B-spline interpolation particle swarm. Background Art

[0002] With the wide popularization of robots in the manufacturing industry, their performance has become a key factor affecting production efficiency and product quality. Therefore, the industry's demand for improving the manufacturing quality and positioning accuracy of robots is increasing day by day. The positioning accuracy of a robot is mainly determined by two aspects: absolute positioning accuracy and repeat positioning accuracy. Absolute positioning accuracy refers to the accuracy with which the end effector of the robot reaches a specified target position in space; while repeat positioning accuracy measures the consistency of the robot in returning to the same position when performing the same task. The latter is mainly affected by the manufacturing process and mechanical structure of the robot. A large number of studies have shown that due to factors such as machining errors of components, installation deviations, and wear during long-term use, there are deviations between the actual mechanical structure parameters and the theoretical parameters of industrial robots. Such deviations are collectively referred to as geometric parameter errors, that is, D-H parameter errors. Research shows that geometric parameter errors account for about 90% of the total absolute positioning error of the robot. Therefore, accurately calibrating geometric parameters is the key to improving the positioning accuracy of the robot, making it an important research direction in the field of precision optimization.

[0003] To address this challenge, researchers are working on developing advanced error modeling techniques, optimization algorithms, and compensation strategies to further improve the applicability of robots in high-precision manufacturing tasks. The current robot calibration methods under several traditional algorithms mainly have the following defects: (1) Gradient-based algorithms, such as the least squares method LS, Levenberg-Marquardt (LM), have truncation errors and are prone to falling into local optima rather than global optima; (2) Filter-based algorithms, such as the extended Kalman filter EKF, particle filter PF, unscented Kalman filter UKF, highly depend on the selection of the covariance matrix; EKF and UKF assume that the noise is Gaussian distributed, which limits their effectiveness in complex environments; (3) Step size update algorithms, such as the beetle antennae search algorithm BAS, etc., have a convergence speed that varies with step size adjustment; without an effective update strategy, they may still fall into local optima; (4) Swarm intelligence algorithms, such as particle swarm optimization PSO, genetic algorithm GA, squirrel search algorithm SSA, etc., may lose diversity due to population aggregation, resulting in falling into local optima. Summary of the Invention

[0004] The present invention introduces a momentum term into the traditional Levenberg - Marquardt algorithm to suppress overshoot oscillation phenomena, incorporates B - spline interpolation technology into the classical Particle Swarm Optimization (PSO) algorithm, cascades these two improved algorithms, realizes a more accurate estimation of the geometric parameters of the robot, and provides a robot calibration method based on Adaptive Momentum LM (AMLM) cascaded with B - spline Interpolation Particle Swarm (BIPSO).

[0005] To achieve the above - mentioned invention purpose, the embodiments of the present invention provide the following technical solutions:

[0006] A robot calibration method based on Adaptive Momentum LM cascaded with B - spline Interpolation Particle Swarm includes the following steps:

[0007] Step 1: Based on the Denavit - Hartenberg model, construct a kinematic error model. According to the robot calibration principle, construct a loss function to optimize the D - H parameters of the kinematic error model.

[0008] Step 2: Based on the kinematic error model, construct a least - squares objective function, and use the AMLM algorithm to obtain a sub - optimal solution of the D - H parameters.

[0009] Step 3: Take the sub - optimal solution as the central value of the initial population of BIPSO, generate a candidate solution set in the BIPSO algorithm, and use the BIPSO algorithm to obtain the optimal solution of the D - H parameter error of the robot.

[0010] Compared with the prior art, the beneficial effects of the present invention are:

[0011] (1) To improve the robot positioning accuracy and enhance the manufacturing quality, this solution proposes a robot geometric error calibration method that combines Adaptive Momentum Levenberg - Marquardt and cascaded B - spline Interpolation Particle Swarm Optimization algorithms. First, by introducing a momentum term into the traditional Levenberg - Marquardt algorithm to suppress overshoot oscillation phenomena, a refined initial calibration of geometric errors is achieved. Second, by incorporating B - spline interpolation technology into the classical Particle Swarm Optimization (PSO) algorithm, the calibration process of robot geometric parameters is further optimized. By cascading these two improved algorithms, a more accurate estimation of geometric parameters is realized. Finally, calibration experiments are carried out on industrial robots to verify the effectiveness and accuracy of the proposed method. Experimental results show that this method significantly improves the robot positioning accuracy.

[0012] (2) In this solution, by introducing adaptive momentum into the traditional LM algorithm, AMLM takes into account historical update steps and dynamically adjusts the optimization direction. This improvement effectively suppresses over - correction, accelerates convergence, and improves the accuracy of robot parameter identification.

[0013] (3) By introducing B-spline interpolation into the standard PSO algorithm, this solution reconstructs the population distribution in BIPSO, ensures the continuity of the fitness landscape, and reduces the dependence on discrete particle positions. This modification significantly improves the calibration accuracy and convergence efficiency.

[0014] (4) The AMLM-BIPSO cascade in this solution further improves the calibration accuracy and convergence. By combining the fast convergence of AMLM with the enhanced global search ability of BIPSO, the proposed cascade method achieves excellent calibration performance. Compared with the state-of-the-art methods, AMLM-BIPSO improves the RMSE, Std, and Max error reduction by 10.75%, 11.11%, and 7.79% respectively. Description of the Drawings

[0015] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as limiting the scope. For those of ordinary skill in the art, other related drawings can be obtained based on these drawings without creative efforts.

[0016] Figure 1 Schematic diagram of the six-joint robot according to the embodiment of the present invention;

[0017] Figure 2 Calibration accuracy and total calculation time of different methods according to Embodiment 2 of the present invention, Figure 2 where (a) shows the root mean square error of different methods, Figure 2 where (b) shows the standard deviation of different methods, Figure 2 where (c) shows the maximum error of different methods, Figure 2 where (d) shows the total calculation time of different methods;

[0018] Figure 3 Training progress of different methods according to Embodiment 2 of the present invention;

[0019] Figure 4 Residual position error after using different methods according to Embodiment 2 of the present invention;

[0020] Figure 5 Flowchart of the method of the present invention. Detailed Embodiments

[0021] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Components of the embodiments of the present invention described and illustrated in the accompanying drawings here can be arranged and designed in various different configurations. Therefore, the detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed present invention, but merely represents selected embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative efforts fall within the scope of protection of the present invention.

[0022] It should be noted that similar reference numerals and letters denote similar items in the following drawings. Therefore, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings. At the same time, in the description of the present invention, terms such as "first", "second", etc. are only used for distinguishing descriptions and cannot be understood as indicating or implying relative importance or suggesting any such actual relationship or order between these entities or operations. In addition, terms such as "connected" and "coupled" can be directly connected between components or indirectly connected via other components.

[0023] Embodiment 1:

[0024] A robot system is as Figure 1 shown, and its main body is an ABB IRB 120 industrial robot. The operation of the robot is in the task planning stage. According to the requirements, the target path is designed, and then it is converted into code executable by the robot through post-processing and uploaded to the control system for execution. However, during long-term operation, factors such as mechanical wear, structural deformation, and thermal effects will cause offsets in the nominal D-H (Denavit-Hartenberg) parameters that define the kinematic structure of the robot. Such offsets cause deviations between the theoretical path and the actual execution path, directly affecting the machining accuracy and quality, specifically manifested as machining defects such as trajectory deviation, weld misalignment, uneven weld width, or poor fusion. To alleviate such problems, it is necessary to systematically determine the actual D-H parameters through robot calibration technology, accurately identify and compensate for kinematic parameter deviations. This process can significantly reduce the trajectory error, ensure a high degree of coincidence between the execution path and the planned path, and thus improve the machining accuracy and consistency.

[0025] Figure 1In the formula, {B} represents the robot base coordinate system (x0, y0, z0); {F} represents the coordinate system (x6, y6, z6) at the flange center of the robot end effector; {E} represents the coordinate system at the end of the wire encoder. When calibrating using the wire encoder, the position error of the robot end effector can be converted into the measurement error of the wire length. By accurately measuring the change in the wire length, the position deviation of the robot end effector can be further deduced and analyzed, so as to realize quantitative error evaluation and compensation. Therefore, the following relational formula can be established:

[0026] (1)

[0027] In formula (1), represents the deviation between the actual length and the theoretical length of the wire; represents the position of the flange center of the end effector in the robot base coordinate system, which is obtained through the robot kinematic model and contains the geometric parameter errors to be identified; represents the position of the end of the wire encoder in the robot base coordinate system, which is directly provided by the measurement software; represents the actual measured wire length obtained by the wire encoder; represents the two-norm.

[0028] Figure 1 The ABB IRB 120 robot in [[ ]] is a six-joint robot, and the joint coordinate systems are respectively represented as: the first joint coordinate system (x1, y1, z1), the second joint coordinate system (x2, y2, z2), the third joint coordinate system (x3, y3, z3), the fourth joint coordinate system (x4, y4, z4), the fifth joint coordinate system (x5, y5, z5), the sixth joint coordinate system (x6, y6, z6). The sixth joint is also the end effector, specifically the coordinate system at the flange center of the end effector. The homogeneous transformation of the connecting rod between adjacent joints is represented as: the connecting rod between the previous joint r and the next joint r + 1 is the connecting rod r. Table 1 gives the initial kinematic parameters of this robot, and q1~q6 are the joint angles, which are directly read by the robot teach pendant.

[0029] Table 1 Initial kinematic parameter table of ABB IRB 120 robot

[0030]

[0031] Robot calibration includes four basic stages: modeling, measurement, parameter identification and error compensation. Among them, modeling plays a decisive role in improving the positioning accuracy of the robot. The Denavit-Hartenberg (D-H) model has long been used as the standard method for robot kinematic modeling. Therefore, this solution uses D-H parameters to construct the kinematic model.

[0032] The present invention is realized through the following technical solutions. As Figure 5 shown, a robot calibration method based on an adaptive momentum LM cascaded B-spline interpolation particle swarm includes the following steps:

[0033] Step 1: Based on the Denavit-Hartenberg model, construct a kinematic error model. According to the robot calibration principle, construct a loss function to optimize the D-H parameters of the kinematic error model.

[0034] Based on Figure 1 , using the four independent parameters in Table 1, the attitude of the (r + 1)-th joint coordinate system relative to the r-th joint coordinate system can be uniquely determined. Therefore, the transformation matrix between adjacent joint coordinate systems is expressed as:

[0035] (2)

[0036] In formula (2), represents the transformation matrix between the r-th joint coordinate system and the (r + 1)-th joint coordinate system; represents the length of link r; represents the offset distance of link r; represents the twist angle of link r; represents the rotation angle of link r.

[0037] For a six-joint robot, the total transformation matrix of the end effector relative to the base coordinate system is:

[0038] (3)

[0039] The relationship between the end effector position and the joint variables is expressed as:

[0040] (4)

[0041] In formula (4), is a rotation column vector; is a translation column vector, that is, the end effector position.

[0042] From this, the expression of the end effector position can be obtained:

[0043] (5)

[0044] The differential form of formula (5) is:

[0045] (6)

[0046] In formula (6), Denote the robot positioning error; Denote the Jacobian matrix; Denote the D-H parameter error of the robot.

[0047] (7)

[0048] In Equation (7), Denote the length error of link r; Denote the offset distance error of link r; Denote the twist angle error of link r; Denote the rotation angle error of link r; r = 1, 2, ..., 6.

[0049] According to the robot calibration principle, the following objective function can be constructed Optimize the D-H parameters:

[0050] (8)

[0051] In Equation (8), Denote the deviation function between the actual length and the theoretical length of the wire for the l-th sample data; L represents the number of sample data, l = 1, 2, …, L; D represents the nominal geometric parameters of the robot; X is the D-H parameter error of the robot.

[0052] Step 2, construct a least-squares objective function based on the kinematic error model, and use the AMLM algorithm to obtain a sub-optimal solution of the D-H parameters.

[0053] The traditional LM algorithm solves the non-linear least-squares problem by seeking a balance between the gradient descent method and the Newton method, and thus becomes a commonly used method for robot kinematic parameter identification. Its core mechanism lies in adjusting the damping factor , ensuring that the update step size has both local quadratic convergence and global stability. However, there are two prominent problems in the actual application of robot D-H parameter error calibration: First, when the curvature of the objective function changes violently or there are local flat regions, single-step updates are likely to lead to convergence to a sub-optimal local solution or oscillation phenomena, affecting the accuracy of kinematic parameter identification; Second, the LM algorithm usually relies on local quadratic approximation and may fail in highly non-linear kinematic calibration problems, especially when the condition number of the Jacobian matrix is high, which will lead to deterioration of numerical stability. Therefore, this solution introduces a momentum mechanism and uncertainty estimation into the LM algorithm to create the AMLM algorithm for robot kinematic parameter identification, realizing adaptive adjustment of the update direction, improving convergence stability, and ensuring higher accuracy in D-H parameter optimization.

[0054] The update step size formula of the traditional LM algorithm is:

[0055] (9)

[0056] In Equation (9), represents the increment vector of the D - H parameter error at the k - th step of update; represents the D - H parameter error vector at the k - th step; represents the Jacobian matrix; represents the identity matrix; is the objective function in Equation (8), represents the position error vector at the k - th step; represents the damping factor.

[0057] After introducing momentum into the traditional LM algorithm, a new AMLM algorithm is created, and its update step formula is:

[0058] (10)

[0059] In Equation (10), represents the D - H parameter error vector updated after adding momentum at the k - th step; represents the momentum coefficient at the k - th step, , and its value is adaptively adjusted according to the current local curvature of the objective function and the descent performance of the objective function.

[0060] Momentum coefficient adjustment strategy: When the current step k has good descent performance (such as the actual descent amount is significantly higher than the expected descent amount), can be appropriately increased to utilize the advantage of the previous direction to accelerate convergence; conversely, when the descent performance is poor, should be decreased to make the update more dependent on the gradient information calculated in the current iteration, and avoid the algorithm from drifting in the wrong direction. Update using the following strategy:

[0061] (11)

[0062] In Equation (11), represents the preset maximum value, such as ; exp(.) represents the natural exponential function; represents the momentum scaling factor; represents the error reduction rate at the k - th step, and its value is determined by the following formula:

[0063] (12)

[0064] In Equation (12), represents the position error vector at the k - th step; represents the position error vector at the (k - 1) - th step. If , it indicates that the error has decreased significantly and a large momentum is maintained; if , indicating that the error decreases slowly and the momentum is reduced to prevent oscillations.

[0065] Therefore, the position update formula of the AMLM algorithm is as follows:

[0066] (13)

[0067] In formula (13), It represents the suboptimal solution obtained by the AMLM algorithm, that is, the k+1th step DH parameter error vector obtained by the AMLM algorithm.

[0068] Step 3: Use the suboptimal solution as the center value of the BIPSO initial population, generate a candidate solution set in the BIPSO algorithm, and use the BIPSO algorithm to obtain the optimal solution of the robot DH parameter error.

[0069] The suboptimal solution of the k+1th step obtained by the AMLM algorithm As the center value of the initial population of the BIPSO algorithm, the population iteration of the BIPSO algorithm is performed to obtain the optimal solution G of the robot DH parameter error in the k+1 step best The population consists of N individuals, each with dim dimensions. The initial position of each individual is determined by the following formula:

[0070] (14)

[0071] In formula (14), i=1,2,…,N; j=1,2,…,dim; represents the position of the i-th individual in the j-th dimension; represents the central value of the i-th individual in the j-th dimension; represents the upper limit of the i-th individual in the j-th dimension; represents the lower limit of the i-th individual in the j-th dimension; rand represents a random number generated in the interval [0,1].

[0072] Calculate the fitness function of each individual:

[0073] (15)

[0074] In formula (15), represents the fitness function of the i-th individual; Represents the position of the i-th individual in all dimensions.

[0075] Update the optimal position of each individual to:

[0076] (16)

[0077] In formula (16), Denote the optimal position of individual \(i\) in all dimensions at the \(t\)-th iteration; Denote the position of individual \(i\) in all dimensions at the \(t\)-th iteration; Denote the position error vector of individual \(i\) in all dimensions at the \(t\)-th iteration.

[0078] Update the global optimal solution as:

[0079] (17)

[0080] In formula (17), \(G\) t Denote the global optimal solution at the \(t\)-th iteration; The function represents finding The minimum value of the objective function, where \(i = 1, 2, \ldots, N\), and assign the corresponding To \(G\) t .

[0081] Update the velocity and position of the individual using the PSO equation:

[0082] (18)

[0083] In formula (18), Denote the velocity of individual \(i\) in all dimensions at the \((t + 1)\)-th iteration; Denote the velocity of individual \(i\) in all dimensions at the \(t\)-th iteration; Denote the inertia weight; \(c1\) denotes the individual learning factor; \(c2\) denotes the social learning factor; \(r1\) and \(r2\) are random variables obeying the uniform distribution \(U(0, 1)\); Denote the position of individual \(i\) in all dimensions at the \(t\)-th iteration; Denote the position of individual \(i\) in all dimensions at the \((t + 1)\)-th iteration.

[0084] Take the positions of the updated \(N\) individuals As \(N\) value points respectively, that is, the fitted B-spline interpolation curve must pass through these value points:

[0085] (19)

[0086] In formula (19), \(P_t\) denotes the \(N\) value points constructed at the \((t + 1)\)-th iteration.

[0087] Using the \(N\) Values in \(P_t\) as value points, interpolate \(M\) interpolation points, and \(M\) is the total number of points generated by resampling or smoothing. The goal is to given \(N\) value points , \(i = 1, 2, \ldots, N\), use the B-spline interpolation function to generate a smooth and continuous interpolation curve and sample \(M\) interpolation points on it.

[0088] Using the chord length parameter method, obtain the parameter vector :

[0089] Initial chord length:

[0090] (20)

[0091] Accumulated chord length:

[0092] (21)

[0093] Normalize the parameter vector t i to the interval [0, 1]:

[0094] (22)

[0095] In Equation (22), represents the total accumulated chord length.

[0096] Since the B-spline interpolation curve is required to pass through all the given N shape points Pt and interpolate M interpolation points, it is necessary to determine N control points C b (b = 0, 1, …, N - 1), that is, the number of control points is equal to the number of shape points. So that:

[0097] (23)

[0098] In Equation (23), represents the B-spline interpolation curve, and there is:

[0099] (24)

[0100] In Equation (24), represents the b-th B-spline basis function of order o with respect to the parameter t i and is recursively defined as:

[0101] (25)

[0102] In Equation (25), t b , t b+1 , t i , t b+o , t b+o+1 represent the b-th, (b + 1)-th, i-th, (b + o)-th, and (b + o + 1)-th parameter vectors respectively.

[0103] The key equation to be solved is:

[0104] (26)

[0105] In Equation (26), A represents ti Evaluated B-spline basis function matrix; the \(i\)-th row in \(A\) is expressed as:

[0106]

[0107] \(C\) represents the vector of control points:

[0108]

[0109] \(P\) t represents the vector of given value points:

[0110]

[0111] By solving Equation (26), the vector \(C\) of control points can be obtained.

[0112] After solving the vector \(C\) of control points, \(M\) interpolation points are generated, and the \(m\)-th new interpolation point is calculated as:

[0113] (27)

[0114] In Equation (27), represents the position of the \(m\)-th interpolation point; \(t\) m represents the interpolation parameter of uniform segmentation; \(M\) represents the number of interpolation points.

[0115] By Equation (27), \(M\) interpolation points (i.e., the candidate solution set) are generated, enhancing the exploration ability of the Particle Swarm Optimization (PSO) algorithm.

[0116] Calculate the fitness function of each interpolation point in the candidate solution set :

[0117] (28)

[0118] The \(M\) fitness values of the interpolation points (fitness function) are sorted in ascending order:

[0119]

[0120] Allocate the first \(N\) values to :

[0121] (29)

[0122] (30)

[0123] In Equation (29), represents sorted in ascending order. represents the 1st to \(N\)th values sorted in ascending order.

[0124] Assign the minimum fitness value to G best :

[0125] (31)

[0126] In Equation (31), represents the first value sorted in ascending order; G best represents the final optimization result of the BIPSO algorithm, which is the error of the robot D-H parameters finally identified.

[0127] In summary, this solution integrates the advantages of the AMLM algorithm and the BIPSO algorithm through the AMLM-BIPSO algorithm, enhancing the accuracy of robot D-H parameter error identification. The AMLM algorithm improves the local search efficiency, reduces overcompensation and oscillation by adding momentum for dynamic adjustment of updates, thereby improving the stability and accuracy of convergence. At the same time, the BIPSO algorithm reconstructs the population distribution through B-spline interpolation, ensuring the continuity of the search space and enhancing the global exploration ability. The cascade of AMLM and BIPSO effectively balances local optimization and global optimization, achieving high-precision D-H parameter error identification, while accelerating convergence and improving the overall robustness.

[0128] Example 2:

[0129] In this example, an experimental verification is carried out on the technical solution of Example 1. The experimental device includes an ABB IRB120 industrial robot (arm) with a repeatability of 0.01 mm. In addition, a wire-pulling encoder with a resolution of 0.004 mm and a host computer are equipped.

[0130] To ensure full coverage of the robot's workspace, sampling points are selected. A total of 1042 points are manually taught, covering the entire workspace of the robot. Subsequently, these sampling points are measured using the wire-pulling encoder. The wire-pulling encoder is used to display the cable length, and at the same time record the position coordinates, joint angles and corresponding cable lengths of the robot end effector. The data is saved in the host computer. After completing the data acquisition, a calibration system is applied to process the collected data set to determine the optimal kinematic parameters.

[0131] The main concern of this solution is the accuracy of robot calibration. The root mean square error (RMSE), standard deviation (STD) and maximum error (MAX) are used as key evaluation criteria:

[0132] (26)

[0133] In Equation (26), L i represents the actual observed value; represents the predicted value; n represents the number of samples.

[0134] To demonstrate the advantages of this solution, this embodiment compares and analyzes this solution with multiple existing algorithms:

[0135] Model M1: Extended Kalman filter (EKF);

[0136] Model M2: Levenberg-Marquardt (LM);

[0137] Model M3: Particle Swarm Optimization (PSO);

[0138] Model M4: Unscented Kalman filter and variable step-size Levenberg-Marquardt (UKF-VSLM);

[0139] Model M5: Extended Kalman filter and Improved Covariance Matrix Adaptive Evolution Strategy (EKF-ICMA-ES);

[0140] Model M6: Levenberg-Marquardt and Beetle Antennae Search Algorithm (LM-BASA);

[0141] Model M7: Levenberg-Marquardt with momentum correction (AMLM);

[0142] Model M8: PSO algorithm combined with B-spline interpolation;

[0143] Model M9: This solution (AMLM-BIPSO).

[0144] Table 2 shows the kinematic parameters calibrated by this solution; Table 3 shows the calibration results and calculation times of Models M1 - M9; Figure 2 It shows the calibration accuracy and total calculation time of different methods, Figure 2 (a) in it shows the root mean square error of different methods, Figure 2 (b) in it shows the standard deviation of different methods, Figure 2 (c) in it shows the maximum error of different methods, Figure 2 (d) in it shows the total calculation time of different methods; Figure 3 It shows the training progress of different methods; Figure 4 It shows the residual position error after using different methods.Figure 2 and Figure 4 In it, BC represents Before Calibration, which stands for the method before calibration.

[0145] Table 2 Kinematic parameter table for calibration of this scheme

[0146]

[0147] Table 3 Calibration results and calculation time table of models M1 - M9

[0148]

[0149] The main findings of the experimental results are summarized as follows:

[0150] (a) The AMLM - BIPSO method (model M9) shows significant improvement in calibration accuracy. After calibration, the RMSE, Std, and Max of M9 are reduced to 0.332, 0.272, and 0.840 respectively, compared with 2.177, 2.006, and 3.99 before calibration, with improvements of 84.75%, 86.44%, and 78.95% respectively. These significant improvements indicate that the AMLM - BIPSO method is highly effective in error suppression and precision optimization, providing a solid technical foundation for the geometric parameter calibration of robot systems.

[0151] (b) Among the evaluated models M1 - M8, M9 shows the highest accuracy in robot calibration. As shown in Table 3 and Figure 2 , the RMSE of M9 is 0.332, the Std is 0.272, and the Max is 0.840; compared with the most accurate alternative method M7, whose RMSE is 0.372, Std is 0.306, and Max is 0.911, M9 has improvements of 10.75%, 11.11%, and 7.79% respectively. These results highlight the effectiveness of the proposed method in improving the accuracy of robot calibration.

[0152] (c) As shown in Table 3 and Figure 2As shown, the AMLM-BIPSO cascading method M9 proposed in this solution is superior to M3 (140.216), M5 (120.643), and M8 (160.258) in terms of computational efficiency, but slightly inferior to M1 (89.635), M2 (87.512), M4 (92.147), M6 (106.215), and M7 (92.156). The lower time costs of M1 and M2 indicate that the computational burden of a single method is relatively light, while the significantly higher costs of M3 and M8 indicate that some single methods can bring a relatively large computational burden. Among the cascading methods, the efficiency of M9 is comparable to that of M4 and M6, but superior to M5. These results show that although M9 incurs a moderate computational cost due to its cascading structure, it is still more efficient than some methods while maintaining the advantages of the cascading method.

[0153] (d) M9 is significantly superior to other calibration methods in terms of convergence efficiency. As Figure 3 shown, the RMSE position error of M9 drops rapidly in the first few iterations and stabilizes at a lower error level than other methods. In contrast, single methods such as M1, M2, M3, M7, and M8 have a slower convergence speed and a larger final error, indicating that their calibration effects are limited. Among the cascading methods (M4, M5, and M6), M9 still achieves faster convergence and a lower final RMSE, demonstrating its excellent optimization ability. These results highlight the effect that M9 improves calibration accuracy while maintaining high convergence efficiency.

[0154] (e) To ensure accurate robot calibration, M9 was used for kinematic parameter identification, and the results are shown in Table 2. In addition, 120 measurement points were collected from an ABB IRB120 industrial robot. After completing the calibration process, a comparative analysis of the calibration results of different methods was carried out, as Figure 4 shown. The results show that the robot positioning error decreases significantly after calibration. Among all methods (M1–M8), M9 achieves the highest calibration accuracy, demonstrating the effectiveness of the proposed method.

[0155] (f) As shown in Table 3 and Figure 2 、 Figure 3 、 Figure 4 shown, M7 (AMLM) is an improved version of M2 (LM). By introducing adaptive momentum, it significantly improves the calibration accuracy and convergence efficiency. By introducing this mechanism, M7 effectively accelerates convergence and reduces the calibration error compared with M2. The results show that M7 achieves a lower RMSE and a faster error reduction, especially in the early iterations, highlighting its improved optimization ability. These findings confirm that integrating adaptive momentum into the LM algorithm can improve its performance, making M7 a more effective calibration method.

[0156] (g) As shown in Table 3 and Figure 2 、Figure 3 , Figure 4 As shown in Figure 4 , M8 (BIPSO) is an improved version of M3 (PSO), which combines B-spline interpolation. Inspired by the KAN neural network, it enhances the reconstruction of population distribution. This modification improves the search ability by ensuring the continuity of the fitness landscape, thus reducing the dependence on discrete particle positions. Therefore, M8 performs better than M3 in terms of calibration accuracy and convergence speed. However, the introduction of interpolation requires additional particle calculations, resulting in an increase in computational time cost. Despite this trade-off, the improvement in the optimization performance of M8 demonstrates the effectiveness of integrating B-spline interpolation into the PSO framework to improve robot calibration.

[0157] (h) M9 is a cascade of M7 (AMLM) and M8 (BIPSO), which performs better than single methods in both calibration accuracy and convergence speed. As shown in Table 6 and Figure 2 , Figure 3 , Figure 4 As shown in Figure 4 , M9 utilizes the adaptive momentum of M7 and the B-spline interpolation of M8 to further improve the accuracy and convergence efficiency. However, the increased computational complexity leads to a higher time cost. Nevertheless, the improvement in the calibration accuracy and efficiency of M9 proves its effectiveness in high-precision robot calibration.

[0158] As described above, it is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention can easily think of changes or substitutions, which should be covered within the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.

Claims

1. A robot calibration method based on an adaptive momentum LM-level cascade B-spline interpolation particle swarm, characterized in that: It includes the following steps: Step 1: Based on the Denavit-Hartenberg model, construct a kinematic error model. According to the robot calibration principle, construct a loss function to optimize the D-H parameters of the kinematic error model; In the said Step 1, the step of constructing a loss function to optimize the D-H parameters of the kinematic error model according to the robot calibration principle includes: According to the robot calibration principle, the following objective function is constructed Optimize the D-H parameters: （8） In Equation (8), represents the deviation function between the actual length and the theoretical length of the wire for the \(l\)-th sample data; \(L\) represents the number of sample data, where \(l = 1, 2, \ldots, L\); \(D\) represents the nominal geometric parameters of the robot; Step 2: Based on the kinematic error model, construct a least-squares objective function, and use the AMLM algorithm to obtain a sub-optimal solution of the D-H parameters; The said Step 2 specifically includes: The update step formula of the LM algorithm is: （9） In Equation (9), represents the D-H parameter error increment vector for the k-th step update; represents the D-H parameter error vector for the k-th step; represents the Jacobian matrix; represents the identity matrix; is the objective function in Equation (8), represents the position error vector for the k-th step; represents the damping factor; After introducing momentum into the traditional LM algorithm, create a new AMLM algorithm, and its update step formula is: （10） In formula (10), represents the D-H parameter error vector updated after adding momentum at the k-th step; represents the momentum coefficient at the k-th step, ; The position update formula of the AMLM algorithm is as follows: （13） In Equation (13), represents the suboptimal solution obtained by the AMLM algorithm, that is, the D-H parameter error vector at the (k + 1)-th step obtained by the AMLM algorithm; Step 3: Use the sub-optimal solution as the central value of the BIPSO initial population, generate a candidate solution set in the BIPSO algorithm, and use the BIPSO algorithm to obtain the optimal solution of the robot D-H parameter error.

2. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 1, wherein: In the said Step 1, the step of constructing a kinematic error model based on the Denavit-Hartenberg model includes: According to the position deviation of the robot end effector, establish the following relationship: （1） In formula (1), represents the deviation between the actual length and the theoretical length of the cable; represents the position at the center of the flange of the end effector in the base coordinate system of the robot; represents the position at the end of the cable encoder in the base coordinate system of the robot; represents the actual measured cable length obtained by the cable encoder; represents the two-norm; Transformation matrix between adjacent joint coordinate systems It is expressed as: （2） In Equation (2), represents the transformation matrix between the r-th joint coordinate system and the (r + 1)-th joint coordinate system; represents the length of link r; represents the offset distance of link r; represents the twist angle of link r; represents the rotation angle of link r; For a six-axis robot, the total transformation matrix of the end effector's position relative to the base coordinate system is as follows: （3） The relationship between the end effector position and the joint variables is expressed as: （4） In formula (4), is a rotation column vector; is a translation column vector, i.e., the end effector position; Thus, the expression of the end effector position can be obtained: （5） The differential form of Equation (5) is: （6） In formula (6), represents the robot positioning error; represents the Jacobian matrix; represents the D-H parameter error of the robot; （7） In formula (7), represents the length error of link r; represents the offset distance error of link r; represents the torsional angle error of link r; represents the angular displacement error of link r; r = 1, 2,..., 6.

3. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 1, characterized in that: Momentum coefficient The adjustment strategy is as follows: （11） In formula (11), represents a preset maximum value; exp(.) represents the natural exponential function; represents the momentum scaling factor; represents the error reduction rate at the k-th step, and its value is determined by the following formula: （12） In formula (12), represents the position error vector at the k-th step; represents the position error vector at the (k - 1)-th step.

4. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 1, characterized in that: In the said Step 3, the step of using the sub-optimal solution as the central value of the BIPSO initial population and generating a candidate solution set in the BIPSO algorithm includes: Use the sub-optimal solution as the central value of the BIPSO initial population, determine the initial position of each individual in the population, and update the optimal position of the individual according to the calculated individual fitness function; update the speed and position of the individual at the iteration step based on the optimal position; Use the updated position of the individual as the value points to generate a B-spline interpolation curve; Solve the control point vector based on the B-spline interpolation curve, thereby generating a candidate solution set.

5. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 4, characterized in that: The said step of using the sub-optimal solution as the central value of the BIPSO initial population, determining the initial position of each individual in the population, and updating the optimal position of the individual according to the calculated individual fitness function; updating the speed and position of the individual at the iteration step based on the optimal position includes: Take the sub - optimal solution at the (k + 1)-th step obtained by the AMLM algorithm as the central value for generating the initial population of the BIPSO algorithm, and perform population iteration of the BIPSO algorithm. This population consists of N individuals, and each individual has dim dimensions; The initial position of each individual is determined by the following formula: （14） In formula (14), i = 1, 2, …, N; j = 1, 2, …, dim; represents the position of the i-th individual in the j-th dimension; represents the central value of the i-th individual in the j-th dimension; represents the upper limit of the i-th individual in the j-th dimension; represents the lower limit of the i-th individual in the j-th dimension; rand represents a random number generated within the interval [0, 1]; Calculate the fitness function of each individual: （15） In formula (15), represents the fitness function of the i-th individual; represents the position of the i-th individual in all dimensions; Update the optimal position of each individual to: （16） In formula (16), represents the optimal position of individual i in all dimensions at the t-th iteration; represents the position of individual i in all dimensions at the t-th iteration; represents the position error vector of individual i in all dimensions at the t-th iteration; Update the global optimal solution to: （17） In Equation (17), G t represents the global optimal solution at the t-th iteration; The function represents finding the minimum value of the objective function, where i = 1, 2, …, N, and assign the corresponding to the minimum value of the objective function to G t ; Use the PSO equation to update the speed and position of the individual: （18） In formula (18), represents the velocity of individual i in all dimensions at the (t + 1)-th iteration; represents the velocity of individual i in all dimensions at the t-th iteration; represents the inertia weight; c1 represents the individual learning factor; c2 represents the social learning factor; r1 and r2 are random variables following the uniform distribution U(0, 1); represents the position of individual i in all dimensions at the (t + 1)-th iteration.

6. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 5, wherein: The said step of using the updated position of the individual as the value points to generate a B-spline interpolation curve includes: The positions of the updated N individuals are used as N profile points respectively, that is, the fitted B-spline interpolation curve must pass through these profile points: （19） In Equation (19), Pt represents the N value points constructed in the (t + 1)-th iteration; Based on the given N shape value points , where i = 1, 2, …, N, a smooth and continuous interpolation curve is generated using the B-spline interpolation function, and M interpolation points are sampled on it; Using the chord length parameter method, obtain the parameter vector : Initial chord length: （20） [[ID=Z7]]Accumulated chord length: （21） Normalize the parameter vector t i to the interval [0, 1]: （22） In formula (22), represents the total cumulative chord length; The B-spline interpolation curve passes through all the given N shape points Pt and interpolates M interpolation points. It is necessary to determine N control points C b , b = 0, 1, …, N - 1, that is, the number of control points is equal to the number of shape points; such that: （23） In formula (23), representing the B-spline interpolation curve, there is: （24） In formula (24), represents the b-th B-spline basis function of order o with respect to the parameter t i and is recursively defined as: （25） In formula (25), t b , t b+1 , t i , t b+o , t b+o+1 respectively represent the b-th, (b + 1)-th, i-th, (b + o)-th, and (b + o + 1)-th parameter vectors.

7. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 6, characterized in that: The said step of solving the control point vector based on the B-spline interpolation curve and thereby generating a candidate solution set includes: Solve the equation: （26） In formula (26), A represents t i The evaluated B-spline basis function matrix; the i-th row in A is expressed as: C represents the vector of control points: P t Vector representing the given value points: By solving Equation (26), the vector C of the control points can be obtained; After solving the vector C of the control points, generate M interpolation points, and the m-th new interpolation point is calculated as: （27） In formula (27), represents the position of the m-th interpolation point; t m represents the interpolation parameter of uniform segmentation. M interpolation points will be generated through Equation (27), constituting a candidate solution set.

8. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 7, characterized in that: In step 3, the steps of obtaining the optimal solution of the D-H parameter error of the robot by using the BIPSO algorithm include: Calculate the fitness function of each interpolation point in the candidate solution set : （28） Interpolate the fitness function at M fitness value interpolation points Sort in ascending order: Assign the first N values to : （29） （30） In formula (29), means arranging in ascending order; represents the 1st to Nth values sorted in ascending order; Assign the minimum fitness value to G best : （31） In formula (31), represents the first value sorted in ascending order; G best represents the final optimization result of the BIPSO algorithm, which is the error of the robot D-H parameters finally identified.

Images