Kinematic Calibration Method for Large-Sized 7-DoF Hybrid Spray-Painting Robots

: Large-sized seven-degrees-of-freedom (7-DoF) hybrid spray-painting robots combine ample working space and high ﬂexibility, making them lucrative for the spray painting of aircraft and rocket surfaces. However, their kinematic calibration is hindered by gravitational deformation, which problem is addressed in this study by introducing a rigid-ﬂexible coupling error modeling method. The latter combines the ﬁnite element method (FEM) and stiffness matrix method to assess the spatial gravitational deformation of a hybrid robot, which is then introduced into a geometric error model to establish the rigid-ﬂexible coupling error identiﬁcation model. Given many redundant parameters in the identiﬁcation model for 7-DoF robots, these parameters are classiﬁed and simpliﬁed using the nonlinear least-square regularization method for parameter identiﬁcation. Combining the inverse solution of 7-DoF spray-painting robots with dynamic characteristics considered, an error compensation method for 7-DoF robots is proposed. The kinematic calibration test results strongly indicate that position errors are signiﬁcantly reduced with gravity compensation taken into consideration, and error convergence speed increases, demonstrating that the kinematic calibration method is feasible and can effectively improve the accuracy of spray-painting robots. The mean errors in the X-and Y-directions are reduced by 20 and 17%, respectively, compared to the conventional method. The proposed method is instrumental in the accurate kinematic calibration of large-sized 7-DoF hybrid robots.


Introduction
With the rapid development of aviation and spacecraft industries, large-sized aeroand spacecraft production is gradually increasing.Their surface is usually protected from erosion and heat by special coatings [1], with thickness and uniformity satisfying stringent requirements [2].Compared to manual operation, robot spray painting is more stable, highefficient, and uniform, and thus has a higher coating quality, making automated painting a vital development direction [3,4].With the advantages of good dynamic characteristics, increased flexibility, high acceleration, ample working space, and the ability to keep the spray gun perpendicular to the painted surface and avoid obstacles, a hybrid robot with seven degrees of freedom (7-DoF) [5,6] is a lucrative choice for the spray painting of large-sized airplanes and rockets.
However, in practice, the accuracy of spray painting by 7-DoF hybrid robots is significantly deteriorated by the following factors: (1) The robot has many joints, which may have relatively large errors in manufacturing and assembly, and the significant geometric errors can result in a relatively low absolute position accuracy [7,8].
Machines 2023, 11, 20 2 of 23 (2) Due to the large size of the robot, which can extend up to 4 m, gravitational deformation significantly impacts the accuracy of this mechanism.
Thus, improving the kinematic accuracy of this mechanism becomes a problem to be solved.The accuracy problem limits the operational capability of the spray-painting robot and makes it challenging to perform spray-painting tasks with a demanding requirement for accuracy.Therefore, improving the absolute kinematic accuracy of the spray-painting robots is of great significance.
Generally, the above problem can be mitigated by kinematic calibration [9,10].However, for the proposed 7-DoF hybrid spray-painting robot, conventional kinematic calibration, which considers only geometric errors in error modeling, is hardly sufficient because gravitational deformation also significantly impacts the terminal position error of the proposed large-sized robot with a cantilever-like terminal telescopic mechanism.The above gravitational deformation's effect on the terminal error leads to inaccurate terminal posture error modeling.
This paper attempts to solve this problem via the proposed calibration method for the 7-DoF spray-painting robot using a rigid-flexible coupling error model combining a rigid error model and an analytical stiffness model of the robot.Such a combination allows one to consider both geometric and gravity deformation errors in the calibration process.The proposed method's verification proves its higher accuracy compared with the conventional kinematic calibration method.

Related Works
Kinematic calibration of robots has been widely used to improve their kinematic accuracy [9,10].Generally, kinematic calibration involves the following four steps: (i) error modeling, (ii) posture measurement, (iii) parameter identification, and (iv) error compensation.Error modeling establishes the relationship between the kinematic error and terminal actuator posture error.Herein, geometric error [11][12][13] is mainly caused by manufacturing and assembly errors, encoder zero-point errors, and other errors, and is the leading cause of terminal posture error.Kinematic modeling is the foundation of geometric error modeling.In kinematic modeling, the most widely used Denavit-Hartenberg (D-H) kinematic model [14] does not apply to robots with two adjacent parallel joints, as the small angular error in the two parallel joints can lead to a significant error in the terminal position of the robot.To address this issue, Hayati [15] proposed the Modified Denavit-Hartenberg (MD-H) model, where an additional angle β was added around the Y-axis at the parallel joints to avoid the singularity problem caused by the slight deflection of the two parallel axes.Additionally, both the S model [16] and the Complete and Parametrically Continuous(CPC) model [17] were used in kinematic modeling.In establishing geometric errors, the more DoF the robot has, the more redundancy errors are introduced.
Since the redundancy errors will reduce the identifiability of errors, the redundancy errors need to be eliminated.Khalil [18] screened and eliminated redundancy errors by analyzing the correlation of errors between adjacent joints.The measurement method dramatically affects the efficiency and accuracy of the calibration.A variety of devices [19][20][21][22] have been used for the measurement of terminal errors.Laser trackers with an extensive measurement range, high accuracy, and good flexibility are more suitable for measuring large-sized robots with high posture flexibility.Parameter identification is the process of using identification algorithms to solve the error model parameters based on the measured data, and the least square method [23,24], maximum likelihood estimation [25], and genetic algorithms [26,27] are used in the parameter identification process.Patel [28] obtained the error sensitivity of Hexapod by approximate linear error mapping and eliminated redundant parameters based on the obtained error sensitivity.Additionally, for error models with complex covariance between parameters, the least squares method faces the problem of unstable identification.The regularization method [29,30] is an ideal identification algorithm with high robustness.Error compensation is the last step in the calibration of geometric parameters.In this step, the identified errors are compensated into the robot control model, improving the robot terminal accuracy.
The terminal errors of the proposed 7-DoF hybrid spray-painting robot is strongly biased by gravitational deformation, whose effects needs to be accurately estimated and compensated for in the calibration process.For the estimation of gravitational deformation, a static stiffness model is required.The finite element method (FEM) [31,32], virtual joint method (VJM) [33,34], and matrix structure analysis (MSA) [35,36] are the three most commonly used methods to establish the static stiffness model.The gravitational deformation results obtained by using FEM are relatively accurate.Still, the modeling process is complicated and computationally inefficient.The entire FEM model must be re-modeled when the mechanism position or posture changes, significantly affecting the error model's efficiency.Proposed by Gosselin [37,38], the VJM takes the connecting rod as a rigid body and the virtual spring at the joint as the flexibility of the connecting rod and joint.However, when multidimensional equivalent springs are used, the virtual joint's stiffness parameters are difficult to solve.Combining the idea of FEM, MSA represents complex structural parts by using simple approximate units and then translates the analytical stiffness model of robots according to the deformation coordination condition and the superposition principle.However, this method requires profound simplification of structural parts, which affects the modeling accuracy.To solve the problem of modeling accuracy of MSA, Klimchik et al. [39,40] proposed the method of using FEM to extract the stiffness matrix of a complex structure and then solve it by MSA, which promoted the application of the structural matrix stiffness modeling method.Overall, the existing techniques for static stiffness modeling of robots and relevant theories considering the impact of gravity need to be revised, and few studies consider both the geometric error model and the gravitational deformation error.Therefore, based on geometric error kinematic calibration, how to efficiently and accurately estimate gravitational deformation, and consider the gravitational deformation in the establishment of the error model to improve model accuracy has become an urgent problem, especially for large-sized robots.
Error compensation is the last step of calibration.Two methods have been widely used for error compensation [41].The first method involves the modification of structural parameters in the kinematic equations and is usually used in calibration, considering only structural errors [42].The second method requires modification of commands and postures.It is generally used in calibration considering non-geometric errors [43] or when the errors of structural parameters are not included in the kinematic equations [44].Given that the proposed 7-DoF hybrid spray-painting robot has redundant DoFs and numerous compensation methods for error compensation, which are difficult to solve, a comprehensive compensation strategy was proposed by combining the inverse kinematic optimization method.
This study established a rigid-flexible coupling error model for the kinematic calibration of 7-DoF hybrid spray-painting robots.Based on the MD-H stiffness error model, the commercial ANSYS software package was used to extract the stiffness matrix of each part.MSA was used to obtain the analytical stiffness model of hybrid robots, thereby separating the position errors induced by geometric errors and those caused by gravity.Considering gravity's impact, the terminal's position was corrected, improving the kinematic model's accuracy and error compensation performance.To address the problem of multiple solutions for error compensation, combining the inverse kinematic solution-solving method that considers dynamic characteristics, an error compensation method with the same objective as the inverse kinematic optimization was proposed.The rest of the paper is organized as follows: Section 1 introduces the structure of the 7-DoF spray-painting robot and establishes a forward kinematic model for robots; Section 2 establishes a kinematic model with geometric deviation taken into consideration, analyzes the redundancy of geometric errors, and eliminates the redundancy errors in the error model by combining FEM and the structural matrix method; Section 3 establishes a static stiffness model for robots and analyzes the impact of gravity on robots, followed by correction of the error model; In Section 4, given the redundancy feature of the 7-DoF spray painting robot, a 7-DoF inverse kinematic solution-solving strategy involving error compensation is proposed, and then, calibration experiments are conducted to verify the correctness of the proposed method and the experimental results are analyzed.In the last section, conclusions are presented.

Introduction to Robot Configuration
The 3D model of the large-sized 7-DoF redundancy spray-painting robot is shown in Figure 1.The robot contains seven kinematic joints, including two motion joints and five rotation joints, and is a 7-DoF kinematic redundancy robot.
Machines 2023, 11, x FOR PEER REVIEW 4 of 24 the structural matrix method; Section 3 establishes a static stiffness model for robots and analyzes the impact of gravity on robots, followed by correction of the error model; In Section 4, given the redundancy feature of the 7-DoF spray painting robot, a 7-DoF inverse kinematic solution-solving strategy involving error compensation is proposed, and then, calibration experiments are conducted to verify the correctness of the proposed method and the experimental results are analyzed.In the last section, conclusions are presented.

Introduction to Robot Configuration
The 3D model of the large-sized 7-DoF redundancy spray-painting robot is shown in Figure 1.The robot contains seven kinematic joints, including two motion joints and five rotation joints, and is a 7-DoF kinematic redundancy robot.The base of the robot is connected to the pedestal through joint #1, which makes the robot capable of moving as a whole in the direction of the axis of the sprayed workpiece to realize the spraying of the whole surface of the sprayed workpiece.The base is connected to the waist through joint #2, which controls the robot's swing in the lateral direction.The waist is connected to the main arm through joint #3, and the main arm is connected to the forearm through joint #4.Joints #3 and #4 are two rotation joints with their axes parallel and which control the angle of pitch of the main arm and forearm respectively.Joint #5 on the forearm controls the rotational movement of the forearm around its own axis.Joint #6 controls the telescopic movement of the forearm, which can increase its working space and flexibility without increasing the size of the robot.The spray gun is installed at the terminal of the forearm, and the direction of the spray gun is the same as the axis  7 , so joint #7 can control the angle between the pointing direction of the spray gun and the axis of the forearm.Joints #5 and #7 are two rotation joints with axes perpendicular to each other, which can ensure that the spray gun points in any direction in the space.
As can be observed from the 3D model of the robot shown in Figure 1, the drive motor of joint #4 utilizes a parallelogram four-rod mechanism on the main arm to transfer the motion to the forearm.By doing so, the motor on the forearm can be transferred to the head of the main arm, which can reduce the weight on the terminal as well as the inertia in the moving process, thus improving the flexibility and operational stability of the robot.The base of the robot is connected to the pedestal through joint #1, which makes the robot capable of moving as a whole in the direction of the axis of the sprayed workpiece to realize the spraying of the whole surface of the sprayed workpiece.The base is connected to the waist through joint #2, which controls the robot's swing in the lateral direction.The waist is connected to the main arm through joint #3, and the main arm is connected to the forearm through joint #4.Joints #3 and #4 are two rotation joints with their axes parallel and which control the angle of pitch of the main arm and forearm respectively.Joint #5 on the forearm controls the rotational movement of the forearm around its own axis.Joint #6 controls the telescopic movement of the forearm, which can increase its working space and flexibility without increasing the size of the robot.The spray gun is installed at the terminal of the forearm, and the direction of the spray gun is the same as the axis x 7 , so joint #7 can control the angle between the pointing direction of the spray gun and the axis of the forearm.Joints #5 and #7 are two rotation joints with axes perpendicular to each other, which can ensure that the spray gun points in any direction in the space.
As can be observed from the 3D model of the robot shown in Figure 1, the drive motor of joint #4 utilizes a parallelogram four-rod mechanism on the main arm to transfer the motion to the forearm.By doing so, the motor on the forearm can be transferred to the head of the main arm, which can reduce the weight on the terminal as well as the inertia in the moving process, thus improving the flexibility and operational stability of the robot.The parallelogram four-rod mechanism is simplified in the subsequent section to reduce complexity.Since the driving connecting rod of the forearm near the head of the main arm and the forearm itself are opposite sides of the parallelogram, the motion of the two is Machines 2023, 11, 20 5 of 23 similar.Therefore, in kinematic modeling, the four-rod mechanism can be ignored, and only the movement of the forearm relative to the main arm is considered.

Forward Kinematic Modeling
The kinematic homogeneous transformation matrix of the robot is established using the D-H method, as shown in Figure 2.
Machines 2023, 11, x FOR PEER REVIEW 5 of 24 The parallelogram four-rod mechanism is simplified in the subsequent section to reduce complexity.Since the driving connecting rod of the forearm near the head of the main arm and the forearm itself are opposite sides of the parallelogram, the motion of the two is similar.Therefore, in kinematic modeling, the four-rod mechanism can be ignored, and only the movement of the forearm relative to the main arm is considered.

Forward Kinematic Modeling
The kinematic homogeneous transformation matrix of the robot is established using the D-H method, as shown in Figure 2. According to the D-H method, the coordinate system of each rod is shown in Figure 2. Therefore, the relative position relationship between two adjacent rods can be described by four parameters: distance ai−1 from axis Zi−1 to axis Zi along the direction of axis Xi−1, rotating angle αi−1 from axis Zi−1 to axis Zi around axis Xi−1, distance di * from axis Xi−1 to axis Xi along the direction of axis Zi, and rotating angle θi from axis Xi−1 to axis Xi about axis Zi.As a result, the homogeneous transformation matrix from the coordinate system Oi−1-Xi−1Yi−1Zi−1 to Oi−1-XiYiZi can be expressed as sin cos cos cos sin sin sin sin cos si 0 0

T T T T T T
(2)

The MD-H Error Model
Since the third and fourth axes of the 7-DoF robot are parallel, the shortcomings of the D-H model shall be considered in error modeling.Therefore, the D-H model is improved by using 5-parameter MD-H modeling, which adds a rotation parameter βi about the Y-axis to transform axis Zi−1 to Zi, thus avoiding the defects of the D-H model.The other parameter definitions of the improved model are the same as those of the D-H model.The joint deviation di * is set to zero if the axes of two adjacent connecting rods are parallel.The rotation angle βi is set to zero if the axes of two adjacent shores are not According to the D-H method, the coordinate system of each rod is shown in Figure 2. Therefore, the relative position relationship between two adjacent rods can be described by four parameters: distance a i−1 from axis Z i−1 to axis Z i along the direction of axis X i−1 , rotating angle α i−1 from axis Z i−1 to axis Z i around axis X i−1 , distance d i * from axis X i−1 to axis X i along the direction of axis Z i , and rotating angle θ i from axis X i−1 to axis X i about axis Z i .As a result, the homogeneous transformation matrix from the coordinate system

The MD-H Error Model
Since the third and fourth axes of the 7-DoF robot are parallel, the shortcomings of the D-H model shall be considered in error modeling.Therefore, the D-H model is improved by using 5-parameter MD-H modeling, which adds a rotation parameter β i about the Y-axis to transform axis Z i−1 to Z i , thus avoiding the defects of the D-H model.The other parameter definitions of the improved model are the same as those of the D-H model.The joint deviation d i * is set to zero if the axes of two adjacent connecting rods are parallel.The rotation angle β i is set to zero if the axes of two adjacent shores are not parallel.Then the MDH parameters of the robot are shown in Table 1.The coordinate transformation matrix of the improved model can be expressed as Table 1.MD-H parameters of the robot.
Joint No.

Rigid-Flexible Coupling Error Model 4.1. Kinematic Error Model Kinematic Error Model
By taking the differential of Equation ( 3), we can obtain Then, the total differential of Equation ( 4) can be calculated as To simplify Equation ( 5), d i−1 T i needs to be expressed as a function of i−1 T i .Suppose that where D α , D a , D θ , D d* , and D β are coefficient matrices.
Combining Equations ( 5) and (6), we get The differentiation by coordinates yields where ∆ i is the differential transformation matrix of the ith joint in the (i − 1)th joint coordinate system and can be expressed as where dx i, dy i, and dz i correspond to differential translation, and δx i, δy i, and δz i to differential rotation.
Combining Equations ( 7) and (8), we get where The terminal position error of the robot is the sum of the errors of different joints.Since the terminal position of the robot is measured in practical measurement, it is necessary to transform the errors of various joints into the robot terminal.Based on the principle of differential transformation of robots, the differential error transformation matrix from the ith joint coordinate system to the terminal coordinate system can be obtained as where n i 7 , o i 7 , a i 7 , and p i 7 are the values of the transformation matrix from the ith joint coordinate system to the terminal coordinate system.Then, we can obtain the terminal errors induced by joint i as The total terminal error induced by all joints can be obtained by where dr denotes the total error obtained by summing the errors of the D-H parameters.

Analysis of Redundancy Kinematic Error of Robots
According to the correlation between the column vectors of the Jacobian matrix, errors can be divided into: independent, correlated, and nonfunctional errors.Independent error refers to an error that can be identified directly, reflected in the fact that the corresponding column of the error is not linearly correlated with other columns of the Jacobian matrix.The correlated error implies that the corresponding error column is linearly correlated with other columns of the Jacobian matrix.Given the interaction of each other during the identification process, only one of them is retained.The nonfunctional error refers to an error that does not affect the terminal error, and since this parameter cannot be identified, it needs to be eliminated.
The transformation matrix from any joint i coordinate system to the base coordinate system of the robot can be expressed as The transformation matrix from the previous joint coordinate system i−1 to the base coordinate system of the robot is Machines 2023, 11, 20

of 23
The relationship between the two coordinate systems can be expressed as Thus, the column vectors in the transformation matrix can be written as The Jacobian matrix J is obtained by taking the total differential of all the kinematic parameters, and J ai−1 , J αi−1 , J di , J θ1 , and J βi indicate the five columns of the Jacobian matrix corresponding to the five D-H parameters, respectively.Then, the Jacobian matrix J can be expressed as According to the differential kinematics of robots, the column vectors of the Jacobian matrix at any joint can be expressed as the cross-product of vectors.
By combining these equations, we get Substituting the theoretical values of the kinematic parameters into the above equations, the following conditions can be derived: 1. If If α i−1 = 0 and a i−1 = 0, δd i−1 * and δd i * are mutually redundant, and either shall be eliminated, and δβ i shall be introduced for identification.

4.
If θ i = 0 and d i * = 0, δa i−1 and δa i are mutually redundant, and some parameters shall be eliminated.

5.
If θ i = 0 and d i * = 0, δα i−1 and δα i are mutually redundant and δa i−1 and δa i are mutually redundant, and some parameters shall be eliminated.

6.
After eliminating redundant parameters, the D-H parameter error still contains 24 errors.

Static Stiffness Model
The overall structure of the robotic arm is simplified to a structure with ten nodes and ten beam units, as shown in Figure 3.
After eliminating redundant parameters, the D-H parameter error still contains 24 errors.

Static Stiffness Model
The overall structure of the robotic arm is simplified to a structure with ten nodes and ten beam units, as shown in Figure 3.The part between joints #1 and #2 is considered the first beam unit, Unit 1; the part between joints #2 and # 3 is regarded as the second beam unit (Unit 2).The main arm between joints #3 and #4 is the third unit (Unit 3).The connecting rod between joint #3 The part between joints #1 and #2 is considered the first beam unit, Unit 1; the part between joints #2 and # 3 is regarded as the second beam unit (Unit 2).The main arm between joints #3 and #4 is the third unit (Unit 3).The connecting rod between joint #3 and the forearm is the fourth unit (Unit 4).The forearm is regarded as the fifth unit (Unit 5).Unit 8 is the telescopic sleeve of the original robotic arm and the hinge mechanism; the connecting rod between joint #7 and the hinge is considered the tenth unit (Unit 10).Unit 7 is perpendicular to both Units 8 and 6, respectively.Units 3, 4, 5, and 6 form a parallelogram mechanism.The mass of Unit i is assumed to be m i (unit: kg), the length is considered to be l i (unit: m), and, in particular, the length of Unit 9 l 9 = l o + d 6 * (l o is the length of the hinge when d 6 * = 0 holds in the D-H coordinate system).As shown in Figure 4, the finite element global coordinate system is defined as follows: the origin of the finite element global coordinate system o-XY coincides with the origin O 1 of O 1 -X 1 Y 1 Z 1 in the robotic arm D-H coordinate system, the x-direction always coincides with the X 2 -direction of O 1 -X 2 Y 2 Z 2 in the robotic arm coordinate system.The y-direction is opposite to the gravity direction.and the forearm is the fourth unit (Unit 4).The forearm is regarded as the fifth unit (Unit 5).Unit 8 is the telescopic sleeve of the original robotic arm and the hinge mechanism; the connecting rod between joint #7 and the hinge is considered the tenth unit (Unit 10).Unit 7 is perpendicular to both Units 8 and 6, respectively.Units 3, 4, 5, and 6 form a parallelogram mechanism.The mass of Unit i is assumed to be mi (unit: kg), the length is considered to be li (unit: m), and, in particular, the length of Unit 9 l9 = lo + d6 * (lo is the length of the hinge when d6 * = 0 holds in the D-H coordinate system).
As shown in Figure 4, the finite element global coordinate system is defined as follows: the origin of the finite element global coordinate system o-XY coincides with the origin O1 of O1-X1Y1Z1 in the robotic arm D-H coordinate system, the x-direction always coincides with the X2-direction of O1-X2Y2Z2 in the robotic arm coordinate system.The ydirection is opposite to the gravity direction.The local coordinate systems are established based on the following principles.Since each unit contains nodes j and k (both satisfying i < j), node j is set as the origin of the local coordinate system, and the direction from node j to node k is set as the positive direction of the x-axis of the local coordinate system of this unit, and the axis obtained by rotating the x-axis π/2 counterclockwise is set as the y-axis.
Since only the robotic arm's deformation under gravity is considered, the theory of small deformation assumption is used here.The uniform load can be converted into the load on the node by the equivalence principle.For Unit i, when the x-direction of the local coordinate system of the unit is the same as the horizontal direction, the equivalent load on the node is as follows (the nodal force is positive in tension, and the bending moment The local coordinate systems are established based on the following principles.Since each unit contains nodes j and k (both satisfying i < j), node j is set as the origin of the local coordinate system, and the direction from node j to node k is set as the positive direction of the x-axis of the local coordinate system of this unit, and the axis obtained by rotating the x-axis π/2 counterclockwise is set as the y-axis.
Since only the robotic arm's deformation under gravity is considered, the theory of small deformation assumption is used here.The uniform load can be converted into the load on the node by the equivalence principle.For Unit i, when the x-direction of the local coordinate system of the unit is the same as the horizontal direction, the equivalent load on the node is as follows (the nodal force is positive in tension, and the bending moment is positive if it is counterclockwise): Herein, F yj = F yk = 1/2 m i g (the concentrated load is positive, which means that the direction is the same as that shown in Figure 5, and the same for the following bending moment), M j = -1/12 m i gl i , and M k = 1/12 m i gl i .Based on the established global coordinate system and the local coordinate systems of different units, we assume the angle between the local coordinate system of each unit and the finite element global coordinate system as where βi * is the angle between the local coordinate system xiyi of Unit i and the finite element global coordinate system o-XY and is positive if the finite element global coordinate system rotates toward the local coordinate system counterclockwise.In vector β * , β1 * = π/2, and β2 * is a non-zero fixed value, which can be obtained by measuring the geometric relationship of the robotic arm structure.Besides, some elements of the vector β * can be described by the variables of the robotic arm in the D-H coordinate system, and some elements of the vector can be expressed as Based on the geometric relationship of each unit in Figure 3, we get Assuming that the stiffness matrix of Unit i in its respective local coordinate system is Ki, then the stiffness matrix of this unit in the finite element global coordinate system is where matrix Ti is Based on the established global coordinate system and the local coordinate systems of different units, we assume the angle between the local coordinate system of each unit and the finite element global coordinate system as where β i * is the angle between the local coordinate system x i y i of Unit i and the finite element global coordinate system o-XY and is positive if the finite element global coordinate system rotates toward the local coordinate system counterclockwise.In vector β * , β 1 * = π/2, and β 2 * is a non-zero fixed value, which can be obtained by measuring the geometric relationship of the robotic arm structure.Besides, some elements of the vector β * can be described by the variables of the robotic arm in the D-H coordinate system, and some elements of the vector can be expressed as Based on the geometric relationship of each unit in Figure 3, we get Assuming that the stiffness matrix of Unit i in its respective local coordinate system is K i , then the stiffness matrix of this unit in the finite element global coordinate system is Machines 2023, 11, 20 where matrix T i is Since the stiffness of Unit 9 changes with the stretching and rotation (i.e., parameters θ 5 and d 6 * ) of the hinge, the stiffness matrix of Unit 9 can be expressed θ 5 about and d 6 * as Herein, the coefficients can be fitted using the stiffness matrix data of Unit 9 extracted in the MySQL Workbench tool.It is assumed that the displacement and load on the ith node under the finite element global coordinate system are as follows: In Equation ( 27), the unit of u i and v i is m, and the unit of ϕ i is rad.Then, the displacement of the whole structure and the load on it under the finite element global coordinate system can be expressed as Assuming that the overall stiffness matrix is K, then Since only the impact of gravity on each unit is considered, the load on each beam is uniform.Additionally, the y-direction of the finite element global coordinate system is the same as the direction of gravity.Therefore, analyzing the nodes under the finite element global coordinate system is more straightforward.Since there is no load in the x-direction of the system under the global coordinate system, the force in the x-direction of each node under the global coordinate system must be zero, i.e., For all nodes except Node 1, the gravity equivalent load on each node is the overall load on the node under the finite element global coordinate system.As seen from the load-equivalent diagram shown in Figure 6, the overall load of all nodes except Node 1 is  Since each node corresponds to a joint of the robotic arm, and the joint is the origin of each coordinate in the D-H coordinate system.Therefore, under the influence of gravity, the origin of the D-H coordinate system of the joint where Node j is located is translated along the translation vector λi = [ui vi 0] under the finite element global coordinate system and rotated along the direction vector p in the z-direction of the finite element global coordinate system with the angle of rotation φi.
It is assumed that the element's coordinate system of the robotic arm, without taking the effect of gravity into account, is OjXjYjZj and the terminal's element coordinate system with the impact of gravity considered is Oj ʹ Xj ʹ Yj ʹ Zj ʹ .Here, the homogeneous transformation matrix from the coordinate system O7X7Y7Z7 to the coordinate system O7 ʹ X7 ʹ Y7 ʹ Z7 ʹ is defined as 7 T7'.
Given that Node 9 is the origin O7 of the coordinate system O7X7Y7Z7, the translation vector of the coordinate system O7X7Y7Z7 under the finite element global coordinate system is [ ] Likewise, λ9 is expressed as 7 λ9 under the coordinate system O7X7Y7Z7.
Then, we can obtain the homogeneous transformation matrix 7 T7ʹ from the coordinate system O7X7Y7Z7 to the coordinate system O7 ʹ X7 ʹ Y7 ʹ Z7 ʹ as R T (36) where O1×3 is a 0 matrix with one row and three columns.For the end displacement p7 = [σx7 σy7 σz7], parameter 7ʹ R7 is derived as follows: Since Node 1 has fixed support, the impact of F 1 on terminal deflection can be ignored.The displacement and angle of rotation of each node can be obtained by eliminating the corresponding rows and columns of Kq = F using boundary conditions.
Since each node corresponds to a joint of the robotic arm, and the joint is the origin of each coordinate in the D-H coordinate system.Therefore, under the influence of gravity, the origin of the D-H coordinate system of the joint where Node j is located is translated along the translation vector λ i = [u i v i 0] under the finite element global coordinate system and rotated along the direction vector p in the z-direction of the finite element global coordinate system with the angle of rotation ϕ i .
It is assumed that the element's coordinate system of the robotic arm, without taking the effect of gravity into account, is O j X j Y j Z j and the terminal's element coordinate system with the impact of gravity considered is O j X j Y j Z j .Here, the homogeneous transformation matrix from the coordinate system O 7 X 7 Y 7 Z 7 to the coordinate system O 7 X 7 Y 7 Z 7 is defined as 7 T 7' .
Given that Node 9 is the origin O 7 of the coordinate system O 7 X 7 Y 7 Z 7 , the translation vector of the coordinate system O 7 X 7 Y 7 Z 7 under the finite element global coordinate system is Likewise, λ 9 is expressed as 7 λ 9 under the coordinate system O 7 X 7 Y 7 Z 7 .

Rigid-Flexible Coupling Kinematic Error Model
In the rigid-flexible coupling kinematic error model, only the influence of the joint's driving direction error on the gravity model is retained, and the effect of other geometric errors on gravitational deformation is neglected.When gravitational deformation is not considered, the Jacobian matrix can be written as e = 7 ∑ i=0 i J 7 G i E i = Jdr, where e is the error transformed into the seventh coordinate system.However, when gravitational deformation exists, since gravitational deformation is a small quantity of first order, and the deviation caused by the mechanism geometric error is also a small quantity of first order, the variation of gravitational deformation caused by the geometric error is a high-order error and thus can be ignored.Therefore, the gravitational deformation of nominal parameters and the parameters compensated in the driving direction are used for calculation in this study.
When the gravity model is considered, the joint coordinate system needs to be transformed using the homogeneous transformation matrix i−1 T i of gravitational deformation of nominal parameters and the parameters compensated in the driving direction.The errors of the terminal coordinate system can be obtained.
Assuming that then, The error vector described in the terminal coordinate system with gravitational deformation taken into consideration is: where J * is the Jacobian matrix of error with gravitational deformation taken into consideration.After the measured terminal error is estimated, one can determine the geometric error from J * , which, in turn, is calculated via the regularized least squares method [29,30].

Error Compensation Strategy for 7-DoF Hybrid Spray-Painting Robots
After kinematic calibration, an error compensation strategy was proposed for the 7-DoF spray-painting robot.The proposed approach solves the inverse motion with errors based on the optimized inverse motion solution.Therefore, the inverse motion with errors also should be optimized in the same direction as the inverse kinematic solution.
Since the inverse motion with errors could not be solved directly using the analytical formula, the inverse kinematic solution with errors was obtained via the Jacobian matrix iteration method.Since the robot contained seven joints and redundant joints resulted in no inverse of the 6 × 7 Jacobian matrix, only six joints were needed to compensate for the identified errors in calibration.Properly selecting the fixed joint is the key to the iteration process.
In solving the inverse kinematics without errors, Wang et al. [45] used inertia volatility as an optimization index, which was found instrumental in solving the inverse kinematics with errors.For kinematic optimization, determination criteria were set, the inertia load sensitivity of the joint corresponding to the unit end error was defined to represent the change of inertia load when the single joint motion compensates for the unit error of the end, and it was taken as the determination index, providing the basis for solving the inverse kinematic solution with errors.The iteration process was as follows.
Herein, the initial value of the Jacobi iteration is the optimized inverse kinematic solution.For a given terminal position p and a given geometric error r, the joint drive variable q = [q 1 q 2 q 3 q 4 q 5 q 6 q 7 ] T can be computed using the numerical method.
When solving the inertial load sensitivity of the joint corresponding to the unit error, the inertial load sensitivity should be solved first.In solving the inertial load sensitivity, to consider the overall effect of a joint on the inertial load, the average inertial load over the entire working space of the remaining joints needs to be computed when each joint of the mechanism takes different values.In this study, given the coordinate of a joint, the inertial loads of the remaining six joints were obtained, integrated about their coordinates, and divided by the length of the integration interval of the remaining six joints to get the average inertia matrix in the entire working space.The inertia matrix was only related to the robot's posture but not to the translational motion of the base, so only Joints #2, #3, #4, #5, and #7 had to be calculated.With joint #6 taken as an example, the variation of the inertia load on joint #6 in extension was as follows: where M m (q) is the inertia matrix of the current joint coordinates, θ u is the upper joint boundary, and θ l is the lower joint boundary.The conformation shows that joints #3, #4, and #6 are more influential than the others, mainly comparing these three joints.The impact of the inertial load on joints #3, #4, and #6 is shown in Figure 7, where the slope is the inertia sensitivity.
DoF spray-painting robot.The proposed approach solves the inverse motion with errors based on the optimized inverse motion solution.Therefore, the inverse motion with errors also should be optimized in the same direction as the inverse kinematic solution.
Since the inverse motion with errors could not be solved directly using the analytical formula, the inverse kinematic solution with errors was obtained via the Jacobian matrix iteration method.Since the robot contained seven joints and redundant joints resulted in no inverse of the 6 × 7 Jacobian matrix, only six joints were needed to compensate for the identified errors in calibration.Properly selecting the fixed joint is the key to the iteration process.
In solving the inverse kinematics without errors, Wang et al. [45] used inertia volatility as an optimization index, which was found instrumental in solving the inverse kinematics with errors.For kinematic optimization, determination criteria were set, the inertia load sensitivity of the joint corresponding to the unit end error was defined to represent the change of inertia load when the single joint motion compensates for the unit error of the end, and it was taken as the determination index, providing the basis for solving the inverse kinematic solution with errors.The iteration process was as follows.
Herein, the initial value of the Jacobi iteration is the optimized inverse kinematic solution.For a given terminal position p and a given geometric error r, the joint drive variable q = [q1 q2 q3 q4 q5 q6 q7] T can be computed using the numerical method.
When solving the inertial load sensitivity of the joint corresponding to the unit error, the inertial load sensitivity should be solved first.In solving the inertial load sensitivity, to consider the overall effect of a joint on the inertial load, the average inertial load over the entire working space of the remaining joints needs to be computed when each joint of the mechanism takes different values.In this study, given the coordinate of a joint, the inertial loads of the remaining six joints were obtained, integrated about their coordinates, and divided by the length of the integration interval of the remaining six joints to get the average inertia matrix in the entire working space.The inertia matrix was only related to the robot's posture but not to the translational motion of the base, so only Joints #2, #3, #4, #5, and #7 had to be calculated.With joint #6 taken as an example, the variation of the inertia load on joint #6 in extension was as follows: ( ) ( ) where Mm(q) is the inertia matrix of the current joint coordinates, θu is the upper joint boundary, and θl is the lower joint boundary.The conformation shows that joints #3, #4, and #6 are more influential than the others, mainly comparing these three joints.The impact of the inertial load on joints #3, #4, and #6 is shown in Figure 7, where the slope is the inertia sensitivity.Secondly, to compare the single joint motion under the unit end error and define the joint sensitivity index.Since this index is affected by the pose of the end, its solution is complicated.In order to simplify the calculation, the joint zero sensitivity index is defined as (with joint #6 taken as an example): where J 6 (q) denotes the column corresponding to d6 in the driving Jacobi matrix, Then, the inertia load sensitivity of the joint corresponding to the unit end error can be expressed as: The inertia sensitivity of joint 6 is relatively high (the slope is relatively high), as shown in Figure 7.The calculated ratio of joint sensitivity of joint 3, 4 and 6 is 2.369:1.647:1.Joint 6 has the minimum joint sensitivity and the maximum inertia sensitivity, so the unit end error of joint 6 corresponds to the maximum inertia load sensitivity of the joint.In this case, according to the optimization direction of inverse motion, Joint #6 was taken as a fixed joint.Then, the coordinates of the seven joints were obtained, and commands were directly entered into the robot for compensation.

Experimental Methods and Data Analysis 7.1. Procedures
The laser tracker model used in the experiment was Leica AT901-B, with an accuracy of (15 + 6.0 × 10 −6 L * ) µm, where L * was the measured distance, in µm.The structural parameters of the robot are listed in Table 2.  complicated.In order to simplify the calculation, the joint zero sensitivity index is defined as (with joint #6 taken as an example): ||J ( ) || S( ) where J6(q) denotes the column corresponding to d6 in the driving Jacobi matrix, Then, the inertia load sensitivity of the joint corresponding to the unit end error can be expressed as: The inertia sensitivity of joint 6 is relatively high (the slope is relatively high), as shown in Figure 7.The calculated ratio of joint sensitivity of joint 3, 4 and 6 is 2.369:1.647:1.Joint 6 has the minimum joint sensitivity and the maximum inertia sensitivity, so the unit end error of joint 6 corresponds to the maximum inertia load sensitivity of the joint.In this case, according to the optimization direction of inverse motion, Joint #6 was taken as a fixed joint.Then, the coordinates of the seven joints were obtained, and commands were directly entered into the robot for compensation.

Procedures
The laser tracker model used in the experiment was Leica AT901-B, with an accuracy of (15 + 6.0 × 10 −6 L * ) μm, where L * was the measured distance, in μm.The structural parameters of the robot are listed in Table 2.As shown in Figure 9, as the terminal position of the robot could not be measured directly, a positioning plate was added to the robot terminal to mount the target sphere base of the laser tracker, and the coordinate system {e} represented the target sphere center coordinate system, that is, the tool's coordinate system.The relationship between the tool's coordinate system {e} and the terminal coordinate system {A7} was described by a transformation matrix A e .Since the terminal adapter panel and laser tracker reflection sphere had high mechanical processing accuracy and were installed at the farthest distance from the robot base, the slight change in the posture had a slight effect on the final measured error.To simplify the calculation, A e only considered the position but not the posture in this study, being was expressed as follows: A e = [p ex p ey p ez 1] T (44) posture in this study, being was expressed as follows: The tool's coordinate system of industrial robots can be converted to the global coordinate system as where w T0 is the robot base coordinate system under the global coordinate system and 7ʹ Te is the target spherical coordinate system under the robot terminal coordinate system.Therefore, the geometric errors in w T0 and 7ʹ Te also had to be identified during error identification.For this reason, these two errors were also included in establishing the error model to reduce the error introduced in the measurement process.In the measurement point selection process, the points were selected from joint space.Since the robot had redundant degrees of freedom, multiple inverse solutions existed.For each joint, the two endpoints of the effective range of motion and the three joint values taken uniformly were selected, giving a total of 5 7 = 78,125 joint input arrangements for the seven joints, which were used as candidates.They were further optimized, according to the measurement posture optimization algorithm proposed in the study [46] to obtain the final 150 sets of measured postures.The tool's coordinate system of industrial robots can be converted to the global coordinate system as w T e = w T 0 0 T 1 1 T 2 . . . 6T 7 7 T 7 7 T e (45) where w T 0 is the robot base coordinate system under the global coordinate system and 7 T e is the target spherical coordinate system under the robot terminal coordinate system.Therefore, the geometric errors in w T 0 and 7 T e also had to be identified during error identification.For this reason, these two errors were also included in establishing the error model to reduce the error introduced in the measurement process.
In the measurement point selection process, the points were selected from joint space.Since the robot had redundant degrees of freedom, multiple inverse solutions existed.For each joint, the two endpoints of the effective range of motion and the three joint values taken uniformly were selected, giving a total of 5 7 = 78,125 joint input arrangements for the seven joints, which were used as candidates.They were further optimized, according to the measurement posture optimization algorithm proposed in the study [46] to obtain the final 150 sets of measured postures.

Validation of Error Elimination Effectiveness
After the elimination of redundancy error, the variation of each error item in the identification process is shownin Figure 10.

Validation of Error Elimination Effectiveness
After the elimination of redundancy error, the variation of each error item in the identification process is shownin Figure 10.As observed, each error term can be accurately identified after eliminating redundancy error, and the identification process is stable and converges in about ten generations.As a comparison, the variation of each error item in the identification process without eliminating redundancy error is shown in Figure 11.As observed, each error term can be accurately identified after eliminating redundancy error, and the identification process is stable and converges in about ten generations.As a comparison, the variation of each error item in the identification process without eliminating redundancy error is shown in Figure 11.

Validation of Error Elimination Effectiveness
After the elimination of redundancy error, the variation of each error item in the identification process is shownin Figure 10.As observed, each error term can be accurately identified after eliminating redundancy error, and the identification process is stable and converges in about ten generations.As a comparison, the variation of each error item in the identification process without eliminating redundancy error is shown in Figure 11.As observed, when redundancy error is not eliminated, the identification process is unstable, there is no convergence in the error term identification, and over-fitting exists, which indicates the correctness of error term elimination.

Experimental Results
One hundred fifty points were selected in the working space for measurement, and As observed, when redundancy error is not eliminated, the identification process is unstable, there is no convergence in the error term identification, and over-fitting exists, which indicates the correctness of error term elimination.

Experimental Results
One hundred fifty points were selected in the working space for measurement, and the measured results were identified, and the results obtained by using the method without calibration, the technique considering only kinematic calibration, and the calibration method considering gravity compensation were compared.
As shown in Figure 12, the errors in the X-, Y-, and Z-directions obtained by using the technique without calibration and the method considering only kinematic calibration were compared to verify the effectiveness of the proposed kinematic calibration As observed, when redundancy error is not eliminated, the identification pro unstable, there is no convergence in the error term identification, and over-fitting which indicates the correctness of error term elimination.

Experimental Results
One hundred fifty points were selected in the working space for measuremen the measured results were identified, and the results obtained by using the method out calibration, the technique considering only kinematic calibration, and the calib method considering gravity compensation were compared.
As shown in Figure 12, the errors in the X-, Y-, and Z-directions obtained by the technique without calibration and the method considering only kinematic calib were compared to verify the effectiveness of the proposed kinematic calibration To illustrate the correctness of the gravity compensation model and the kin calibration method considering gravitational deformation, the errors in the X-, Y-, directions of the measurement points obtained by using the method considerin To illustrate the correctness of the gravity compensation model and the kinematic calibration method considering gravitational deformation, the errors in the X-, Y-, and Z-directions of the measurement points obtained by using the method considering only kinematic calibration and the kinematic calibration method considering gravitational deformation were compared in Figure 13.
The errors obtained after calibration by each method are listed in Table 3.
For calibration alone, the error in each direction shows a reduction trend, indicating the correctness of the kinematic calibration model.Meanwhile, the error obtained by the calibration method considering gravity is smaller than that obtained by the technique considering only calibration, indicating the correctness of the proposed gravitational deformation model and the calibration method considering gravity.
After eliminating redundancy error, the error terms identified using the proposed calibration method were obtained and listed in Table 4.The errors obtained after calibration by each method are listed in Table 3.For calibration alone, the error in each direction shows a reduction trend, indicating the correctness of the kinematic calibration model.Meanwhile, the error obtained by the calibration method considering gravity is smaller than that obtained by the technique considering only calibration, indicating the correctness of the proposed gravitational deformation model and the calibration method considering gravity.
After eliminating redundancy error, the error terms identified using the proposed calibration method were obtained and listed in Table 4.After calibration, the measurement points' errors were identified, as shown in Table 4.To verify the accuracy in the whole working space after calibration, some other measurement points in the entire working space were used for verification.After compensation with the identified data, 100 measurement points were used here to verify the accuracy after calibration.
For calibration only without considering gravitational deformation, the compensated error and the uncalibrated error of each measurement point were compared using the identified data, as shown in Figure 14.
4. To verify the accuracy in the whole working space after calibration, some other measurement points in the entire working space were used for verification.After compensation with the identified data, 100 measurement points were used here to verify the accuracy after calibration.
For calibration only without considering gravitational deformation, the compensated error and the uncalibrated error of each measurement point were compared using the identified data, as shown in Figure 14.After calibration considering gravitational deformation, the compensated error, the error obtained by calibration without considering gravity, and the error obtained by only calibration of each measurement point were compared using the identified parameters, as shown in Figure 15 and tabulated in Table 5.After calibration considering gravitational deformation, the compensated error, the error obtained by calibration without considering gravity, and the error obtained by only calibration of each measurement point were compared using the identified parameters, as shown in Figure 15 and tabulated in Table 5.

Conclusions
A rigid-flexible coupling kinematic modeling method considering gravitational deformation was proposed for the 7-DoF spray-painting robot.The kinematic model was used to identify the robot's geometric errors, which can simultaneously compensate for its geometric errors and gravitational deformation and improve its terminal accuracy.Firstly, the MD-H model constructed the rigid kinematic error model of the 7-DoF spraying robot.The redundancy of the geometric errors was analyzed, and the redundancy error was eliminated.Then, the stiffness matrix of each part was extracted from the finite element software by considering the part's flexibility.The obtained matrices were joined together by MSA, thereby bringing the static stiffness matrix of the whole robot.Based on this static stiffness matrix, the effect of gravitational deformation on the terminal position of the robot was analyzed, and the MD-H rigid kinematic error model was modified to obtain a rigid-flexible coupling kinematic error model considering gravitational deformation.The regularization identification algorithm was used to identify the geometric errors.Finally, the proposed method was verified by kinematic calibration experiments, which compared the results obtained using the technique without calibration, the method considering only kinematic calibration, and the calibration method considering gravity compensation.The experimental results showed that using the rigid-flexible coupling kinematic model for calibration could achieve the highest terminal accuracy, indicating the correctness of the proposed geometric error calibration model and the gravitational deformation model.Additionally, to address the redundancy problem of the 7-DoF robot, an error compensation strategy was developed.This method can be used to improve the accuracy of large-sized 7-DoF H-Ds and is of great reference value for other mechanisms with the characteristics of the robots.

Figure 1 .
Figure 1.3D model of the robot.For the convenience of description, the seven joints are numbered from 1 to 7, from the base to the terminal of the robot.

Figure 1 .
Figure 1.3D model of the robot.For the convenience of description, the seven joints are numbered from 1 to 7, from the base to the terminal of the robot.

Figure 2 .
Figure 2. Schematic of the robot mechanism.

Figure 2 .
Figure 2. Schematic of the robot mechanism.

Figure 3 .
Figure 3. Structure of the robotic arm.1-9 indicates the node numbers and ①-⑨ the numbering of the beam units.

Figure 3 .
Figure 3. Structure of the robotic arm.1-9 indicates the node numbers and 1 -9 the numbering of the beam units.

Figure 4 .
Figure 4. Finite element global coordinate system of the robotic arm.

Figure 4 .
Figure 4. Finite element global coordinate system of the robotic arm.

Figure 8
Figure 8 presents the calibration flowchart.

Figure 8
Figure 8 presents the calibration flowchart.

Figure 9 .
Figure 9. Calibration test using a Leica laser tracker.

Figure 9 .
Figure 9. Calibration test using a Leica laser tracker.

Figure 12 .
Figure 12.Comparison of kinematic calibration and non-calibration in the X-(a) Y-(b) an directions.

Figure 12 .
Figure 12.Comparison of kinematic calibration and non-calibration in the X-(a) Y-(b) and Z-(c) directions.

Machines 2023 ,Figure 13 .
Figure 13.Comparison of gravity compensation and non-compensation in the X-(a), Y-(b) and Z-(c) directions.

Figure 13 .
Figure 13.Comparison of gravity compensation and non-compensation in the X-(a), Y-(b) and Z-(c) directions.

Table 3 .
Maximum and mean errors obtained after calibration by various methods.

Table 4 .
MD-H parameters of the robot after identification combined with gravity compensation.

Table 3 .
Maximum and mean errors obtained after calibration by various methods.
X-Direction/m Y-Direction/m Z-Direction/m

Table 4 .
MD-H parameters of the robot after identification combined with gravity compensation.

Table 5 .
Verification of calibration considering gravity in X-, Y-, and Z-directions.

Table 5 .
Verification of calibration considering gravity in X-, Y-, and Z-directions.