Robot calibration method for adaptive momentum LM cascade B-spline interpolation particle swarm

Abstract

The invention relates to the technical field of robot calibration, and discloses an adaptive momentum LM cascade B-spline interpolation particle swarm robot calibration method, which comprises the following steps: constructing a kinematics error model based on an LM model, and constructing a loss function to optimize D-H parameters of the kinematics error model; constructing a least square objective function based on a kinematic error model, and obtaining a suboptimal solution of the D-H parameter by adopting an AMLM algorithm; and taking the suboptimal solution as a central value of a BIPSO initial population, generating a candidate solution set in a BIPSO algorithm, and obtaining an optimal solution of a robot D-H parameter error through the BIPSO algorithm. According to the method, a momentum item is introduced into a traditional LM algorithm to suppress an overshoot oscillation phenomenon, a B-spline interpolation technology is fused into a classic particle swarm optimization (PSO) algorithm, and the two improved algorithms are cascaded to realize more accurate estimation of 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. Repeat positioning accuracy, on the other hand, measures the consistency of the robot's ability to return 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 has shown that geometric parameter errors account for approximately 90% of the total absolute positioning error of robots. Therefore, accurately calibrating geometric parameters is the key to improving the positioning accuracy of robots, 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, and realizes a more accurate estimation of the geometric parameters of the robot, providing 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: Use 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 robot D - H parameter error.

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

[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 algorithm. 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. The 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 of 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. BRIEF 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 a six-joint robot according to an 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 DESCRIPTION OF THE EMBODIMENTS

[0021] The following will clearly and completely describe the technical solutions in the embodiments of the present invention 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. The 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 following 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 implying 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 through other components.

[0023] Embodiment 1:

[0024] The 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, robot calibration technology is required to systematically determine the actual D-H parameters, 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 rope encoder. When calibrating using the wire rope encoder, the position error of the robot end effector can be converted into the measurement error of the wire rope length. By accurately measuring the change in the wire rope 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 expression can be established: (1)

[0026] In formula (1), represents the deviation between the actual length and the theoretical length of the wire rope; 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 rope encoder in the robot base coordinate system, which is directly provided by the measurement software; represents the actual measured wire rope length obtained by the wire rope encoder; represents the two-norm.

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

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

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

[0030] The present invention is achieved 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, comprising the following steps:

[0031] Step 1, construct a kinematic error model based on the Denavit-Hartenberg model, and according to the robot calibration principle, construct a loss function to optimize the D-H parameters of the kinematic error model.

[0032] 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: (2)

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

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

[0035] The relationship between the end effector position and the joint variables is expressed as: (4)

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

[0037] From this, the expression of the end effector position can be obtained: (5)

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

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

[0040] (7)

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

[0042] According to the robot calibration principle, the following objective function can be constructed to optimize the D-H parameters: (8)

[0043] In Equation (8), represents the deviation function between the actual length and the theoretical length of the wire rope of 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.

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

[0045] 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 at the same time. 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 areas, 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, it will lead to deterioration of numerical stability. Therefore, in this solution, a momentum mechanism and uncertainty estimation are introduced 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.

[0046] The update step size formula of the traditional LM algorithm is: (9)

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

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

[0049] In Equation (10), denotes the updated D-H parameter error vector after adding momentum at the k-th step; denotes 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.

[0050] 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 to avoid the algorithm drifting in the wrong direction. Update using the following strategy: (11)

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

[0052] In Equation (12), denotes the position error vector at the k-th step; denotes 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 , it indicates that the error has decreased slowly and the momentum is reduced to prevent oscillation.

[0053] Therefore, the position update formula of the AMLM algorithm is as follows: (13)

[0054] In Equation (13), Denote the sub-optimal 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.

[0055] Step 3: Use 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 robot D-H parameter error.

[0056] The sub-optimal solution at the (k + 1)-th step obtained by the AMLM algorithm is used as the central value for generating the initial population of the BIPSO algorithm, and population iteration of the BIPSO algorithm is performed to obtain the optimal solution G of the robot D-H parameter error at the (k + 1)-th step best . The population includes N individuals, and each individual has dim dimensions. The initial position of each individual is determined by the following formula: (14)

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

[0058] Calculate the fitness function of each individual: (15)

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

[0060] Update the optimal position of each individual to: (16)

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

[0062] Update the global optimal solution to: (17)

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

[0064] Update the velocity and position of the individual using the PSO equation: (18)

[0065] In Equation (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; \(c_1\) represents the individual learning factor; \(c_2\) represents the social learning factor; \(r_1\) and \(r_2\) are random variables following the uniform distribution \(U(0, 1)\); represents the position of individual \(i\) in all dimensions at the \(t\)-th iteration; represents the position of individual \(i\) in all dimensions at the \((t + 1)\)-th iteration.

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

[0067] In Equation (19), \(P_t\) represents the \(N\) value points constructed at the \((t + 1)\)-th iteration.

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

[0069] Use the chord length parameter method to obtain the parameter vector :

[0070] Initial chord length: (20)

[0071] Accumulated chord length: (21)

[0072] Normalize the parameter vector \(t\) i to the interval \([0, 1]\): (22)

[0073] In Equation (22), represents the total cumulative chord length.

[0074] Since the B-spline interpolation curve is required to pass through all the given N shape points Pt and interpolate M interpolation points, N control points C b (b = 0, 1, …, N−1) need to be determined, that is, the number of control points is equal to the number of shape points, such that: (23)

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

[0076] 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: (25)

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

[0078] The key equation to be solved is: (26)

[0079] In Equation (26), A represents the matrix of B-spline basis functions evaluated at t i ; the i-th row in A is expressed as:

[0080] C represents the vector of control points:

[0081] P t represents the vector of given shape points:

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

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

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

[0085] Equation (27) will generate M interpolation points (i.e., the candidate solution set), enhancing the exploration ability of the Particle Swarm Optimization (PSO) algorithm.

[0086] Calculate the fitness function of each interpolation point in the candidate solution set : (28)

[0087] Interpolate the fitness values of the M interpolation points Sort them in ascending order:

[0088] Assign the first N values to : (29) (30)

[0089] In Equation (29), means to sort in ascending order. represents the 1st to the N-th values sorted in ascending order.

[0090] Assign the minimum fitness value to G best : (31)

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

[0092] 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 by adding momentum for dynamic adjustment of updates, reducing overcompensation and oscillation, thereby enhancing 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 overall robustness.

[0093] Embodiment 2:

[0094] In this embodiment, an experimental verification is conducted on the technical solution of Embodiment 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 also equipped.

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

[0096] The main focus 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: (26)

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

[0098] To demonstrate the advantages of this solution, this embodiment conducts a comparative analysis of this solution with multiple existing algorithms:

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

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

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

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

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

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

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

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

[0107] Model M9: This scheme (AMLM - BIPSO).

[0108] Table 2 shows the kinematic parameters calibrated by this scheme; Table 3 shows the calibration results and calculation times of Models M1 - M9; Figure 2 The calibration accuracy and total calculation time of different methods are shown, 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 The training progress of different methods is shown; Figure 4 The residual position error after using different methods is shown. Figure 2 and Figure 4 In it, BC represents Before Calibration, representing the method before calibration.

[0109] Table 2 Table of kinematic parameters calibrated by this scheme

[0110] Table 3 Table of calibration results and calculation times of Models M1 - M9

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

[0112] (a) The AMLM-BIPSO method (model M9) showed significant improvement in calibration accuracy. After calibration, the RMSE, Std, and Max of M9 were reduced to 0.332, 0.272, and 0.840 respectively, compared with 2.177, 2.006, and 3.99 before calibration, representing 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 accuracy optimization, providing a solid technical foundation for the geometric parameter calibration of robotic systems.

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

[0114] (c) As shown in Table 3 and Figure 2 as indicated, the proposed AMLM-BIPSO cascade method M9 of this scheme was 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 single methods have a lighter computational burden, while the significantly higher costs of M3 and M8 indicate that some single methods can bring a greater computational burden. Among the cascade methods, the efficiency of M9 was 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 cascade structure, it is still more efficient than some methods while maintaining the advantages of cascade methods.

[0115] (d) M9 was significantly superior to other calibration methods in terms of convergence efficiency. As Figure 3As 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 rate and a larger final error, indicating limited calibration effectiveness. Among the cascaded methods (M4, M5, and M6), M9 still achieves faster convergence and a lower final RMSE, demonstrating its excellent optimization ability. These results highlight the effectiveness of M9 in improving calibration accuracy while maintaining high convergence efficiency.

[0116] (e) To ensure accurate robot calibration, M9 was used for kinematic parameter identification, and the results are shown in Table 2. Additionally, 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 conducted, as Figure 4 shown. The results indicate that the robot positioning error decreased significantly after calibration. Among all methods (M1–M8), M9 achieved the highest calibration accuracy, demonstrating the effectiveness of the proposed method.

[0117] (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 calibration accuracy and convergence efficiency. By introducing this mechanism, M7 effectively accelerates convergence and reduces the calibration error compared to M2. The results show that M7 achieves a lower RMSE and 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 enhance its performance, making M7 a more effective calibration method.

[0118] (g) As shown in Table 3 and Figure 2 , Figure 3 , Figure 4 shown, M8 (BIPSO) is an improved version of M3 (PSO). It combines B-spline interpolation, inspired by the KAN neural network, to enhance population distribution reconstruction. 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, introducing interpolation requires additional particle calculations, resulting in an increased 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 for robot calibration.

[0119] (h) M9 is a cascade of M7 (AMLM) and M8 (BIPSO), and it performs better than single methods in terms of 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, further improving the accuracy and convergence efficiency. However, the increased computational complexity results in a higher time cost. Nevertheless, the improvement in calibration accuracy and efficiency of M9 demonstrates its effectiveness in high-precision robot calibration.

[0120] As described above, only the specific embodiments of the present invention are provided, 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 cascaded B-spline interpolation particle swarm, characterized in that: It includes the following steps: Step 1: Construct a kinematic error model based on the Denavit-Hartenberg model. According to the robot calibration principle, construct a loss function to optimize the D-H parameters of the kinematic error model; 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; 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 robot D-H parameter error.

2. The robot calibration method based on an adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 1, characterized in that: In the said Step 1, the steps of constructing a kinematic error model based on the Denavit-Hartenberg model include: According to the position deviation of the robot end effector, establish the following relationship: （1） In Equation (1), represents the deviation between the actual length and the theoretical length of the wire rope; represents the position at the flange center of the end effector in the robot base coordinate system; represents the position at the end of the wire rope encoder in the robot base coordinate system; represents the actual measured wire rope length obtained by the wire rope encoder; represents the two - norm; Transformation matrix between adjacent joint coordinate systems It is expressed as: （2） 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; 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 position of the end effector 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; From this, the expression of the position of the end effector 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 Equation (7), represents the length error of link r; represents the offset distance error of link r; represents the twist angle error of link r; represents the rotation angle error of link r; r = 1, 2,..., 6.

3. The robot calibration method based on an adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 2, characterized in that: In the said Step 1, the steps of constructing a loss function to optimize the D-H parameters of the kinematic error model according to the robot calibration principle include: According to the robot calibration principle, the following objective function is constructed Optimize the D-H parameters: （8） In formula (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, \(l = 1, 2, \ldots, L\); \(D\) represents the nominal geometric parameters of the robot.

4. The robot calibration method based on an adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 3, characterized in that: In the said Step 2, the steps of constructing a least squares objective function based on the kinematic error model and using the AMLM algorithm to obtain a sub-optimal solution of the D-H parameters include: 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, a new AMLM algorithm is created, 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 formula (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.

5. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 4, 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 Equation (12), represents the position error vector at the k-th step; represents the position error vector at the (k - 1)-th step.

6. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 4, characterized in that: In the said Step 3, the steps of taking the sub-optimal solution as the central value of the initial population of BIPSO and generating a candidate solution set in the BIPSO algorithm include: Take the sub-optimal solution as the central value of the initial population of BIPSO, 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; Take the position of the updated individual as the value point 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.

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 steps of taking the sub-optimal solution as the central value of the initial population of BIPSO, 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 include: The sub-optimal solution at the (k + 1)-th step obtained by the AMLM algorithm is used as the central value for generating the initial population of the BIPSO algorithm, and population iteration of the BIPSO algorithm is performed. 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 Equation (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 subject to the uniform distribution U(0, 1); represents the position of individual i in all dimensions at the (t + 1)-th iteration.

8. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 7, characterized in that: The said steps of taking the position of the updated individual as the value point to generate a B-spline interpolation curve include: 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: （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, use the B-spline interpolation function to generate a smooth and continuous interpolation curve, and sample M interpolation points on it; Using the chord length parameter method, obtain the parameter vector : Initial chord length: （20） 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), denotes the B-spline interpolation curve, and we have: （24） In Equation (24), represents the $b$-th B-spline basis function of order $o$ with respect to the parameter $t$, i which 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.

9. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 8, wherein: The said steps of solving the control point vector based on the B-spline interpolation curve and thereby generating a candidate solution set include: 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 control points can be obtained; After solving the vector C of control points, M interpolation points are generated, 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) to form a candidate solution set.

10. The robot calibration method of the adaptive momentum LM cascaded B-spline interpolation particle swarm according to claim 9, wherein: 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; means the 1st to Nth values arranged 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 and is the error of the robot D-H parameters finally identified.

Images