The optimized algorithm based on machine learning for inverse kinematics of two painting robots with non-spherical wrist

This paper studies the inverse kinematics of two non-spherical wrist configurations of painting robot. The simplest analytical solution of orthogonal wrist configuration is deduced in this paper for the first time. For the oblique wrist configuration, there is no analytical solution for the configuration. So it is necessary to solve by general method, which cannot achieve high precision and high speed as analytic solution. Two general methods are optimized in this paper. Firstly, the elimination method is optimized to reduce the solving speed to 20% of the original one, and the completeness of the method is supplemented. Based on the Gauss damped least squares method, a new optimization method is proposed to improve the solving speed. The enhanced step length coefficient is introduced to conduct studies with the machine learning correlation method. It has been proved that, on the basis of ensuring the stability of motion, the number of iterations can be effectively reduced and the average number of iterations can be less than 5 times, which can effectively improve the speed of solution. In the simulation and experimental environment, it is verified.


Introduction
For the 6R painting robot, non-spherical wrists are usually adopted because of operation requirements. It is mainly divided into orthogonal non-spherical wrist and oblique non-spherical wrist. For the orthogonal non-spherical wrist configuration, three adjacent axes are parallel to each other to meet Pieper criterion, there is an analytical solution, and the simplest analytical solution is given in this paper. For the oblique non-spherical wrist configuration, the Pieper criterion is not satisfied and there is no analytical solution. Among many methods, the elimination method has the most complex calculation process, but it has the highest accuracy and can get all the solutions under the same terminal position. In this paper, the elimination method was optimized to improve the speed of solution and to supplement the integrity of the method. Gaussian damped iteration method is derived by the most widely used damped least squares method. The Gauss damped iteration method is an efficient and stable method. It has the advantages of fast convergence of Newton iteration method, smooth motion of damping method and less error after optimizing damping factor. But in the embedded environment, it cannot achieve the same speed as the analytic solution of general configuration. On the basis of the algorithm, this paper introduced the enhanced step length coefficient, analyzed the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 precision of the solution. An intelligent algorithm based on extreme learning machine and sequential mutation genetic algorithm was proposed in [25]. An online adaptive strategy based on the Lyapunov stability theory was presented in [26]. Radial Basis Function (RBF) Neural Networks (NNs) was employed and Quadratic Programming (QP) method was incorporated in the training algorithm of the NNs. Various soft computing methods can achieve accuracy up to the micron level, but the solution speed is still a problem.
In this paper, the enhanced step length coefficient was studied by machine learning methods. Classification and regression methods commonly used in machine learning mainly include DecisionTree(DT) [27] and NaiveBayes(NB). DT is a prediction model to establish a mapping between object attributes and object values. It has fast speed, and is easy to visualize the model, but easy to overfitting. NB takes the maximum conditional probability based on Bayesian criterion. It needs few parameters to estimate and is insensitive to missing data. However, assuming that attributes are independent of each other, it will affect the accuracy. K-NearestNeighbor(KNN) [28] is based on the limit theorem. Decision-making depends on a very small number of the nearest neighbor samples. It is suitable for multi-classification problems. LogisticRegression(LR) is suitable for continuous and categorical independent variables and is widely used in medical related fields. It is sensitive to multiple collinearity of independent variables in the model. RandomForest(RF) [29] integrates the set of decision trees with control variance by using the idea of bagging. Gradient Boosting Decision Tree(GBDT) [30] generates a weak classifier in each iteration through multiple iterations. Each classifier trains on the basis of the residual of the previous one. Support Vector Machine(SVM) estimates regression based on a series of kernels. In this paper, the commonly used RBF nonlinear kernels are used. Bagging method constructs a series of predictive functions and combines them into a predictive function. ExtraTreeRegressor(ETR) [31] is similar to RF and consists of many decision trees. Each decision tree is obtained by using all training samples, and the bifurcation values are obtained completely randomly. XBGboost combines all the predictions of a group of weak learners to train a strong learner through additive training strategies.

Inverse kinematics of orthogonal non-spherical wrist robot based on analytical method
YASKAWA-EPX2050 (lemma wrist) was taken as an example in this paper. The model is shown in Fig 1, and DH parameters is shown in Table 1.
The transformation matrix between two adjacent coordinate systems of the robot is T i . Positive kinematics can be obtained from Eq (1).
Where T end ¼ l x m x n x p x l y m y n y p y l z m z n z p z By reversible transformation of Eq (1), T 3 T 4 T 5 ¼ ðT 2 Þ À 1 ðT 1 Þ À 1 T end ðT 6 Þ À 1 was obtained.

PLOS ONE
The optimized algorithm for inverse kinematics of two painting robots with non-spherical wrist Then we can get Eqs (2)-(4).
Right and left matrices, element (1,4) was correspondingly equal, obtaining the following equations.
Then we can get Eq (7).
Where g ¼ p z s 23 À d 6 ðn z s 23 þ c 1 n x c 23 þ n y s 1 c 23 Þ À a 1 c 23 À a 3 À a 2 c 3 þ c 1 p x c 23 þ p y s 1 c 23 d 5

Inverse kinematics of oblique non-spherical wrist robot based on two optimization methods
YASKAWA-EPX2800 (lemma wrist) was taken as an example in this paper. The model is shown in Fig 2, and DH parameters is shown in Table 2.
Singularity analysis for oblique non-spherical wrist robot.. Differential transform method was used to construct the Jacobian matrix J of the robot in this structure. The robot in a singular configuration had less degree of freedom in operating space, and the determinant of J was 0. Therefore, the determinant could be used to determine the singular configuration. For simplified calculation, set the coordinate system x 5 À y 5 À z 5 to coincide with the coordinate system x 6 À y 6 À z 6 . Via the complicated calculation and derivation, the determinant was obtained as shown in Eq (8).
According to Eq (8), there were three cases when detðJÞ ¼ 0. By decomposing the Jacobian matrix J, J = USV T was obtained, of which S was a diagonal matrix composed of non-zero singular values of J. The diagonal elements were arranged in descending order as s 1 � s 2 � s 3 � s 4 � s 5 � s 6 , of which the σ 6 referred to the minimum singular value. We then analyzed the relations between the singularities of the configuration and the singular value, and the characteristics of the singular values. One thousand positions were selected on a random basis to analyze the singular values, as shown in Fig 3. 10 6 positions were selected at random. When the singular configuration and non-singular configuration, the mean singular value and variable coefficient were in statistics and shown in Table 3. Fig 3 intuitively showed the bigger change in σ 2 , σ 3 , σ 6 , and smaller change in σ 1 , σ 4 , σ 5 , of which the σ 6 had the greatest change. Therefore, in various damping iteration methods for analyzing the inverse kinematics, the method taking all singular values into account had a better effect than the method of determining damping factors only by the minimum singular value. The statistics data in Table 3 showed that the mean value of minimum singular values varied greatly from variable coefficients for singular configuration and non-singular configuration. Therefore, in analyzing the influence factors of the enhanced step length coefficient, the singularity of the configuration could be roughly represented by the minimum singular values. Table 2. DH parameters of oblique non-spherical wrist robot.
https://doi.org/10.1371/journal.pone.0230790.t002 Optimized elimination method.. The specific process of elimination method was as follows. By reversible transformation of Eq (1), it was obtained that: Right and left matrices, elements in column 3, 4 were correspondingly equal, obtaining Eq 6 unrelated to y 6 .
Via eliminating the Eq 14 e 1 e e 14 by θ 1 , θ 2 , we obtained Eq 6. Wrote Eq 6 in the form of a matrix equation, obtaining the matrix equation 1þx 5 2 , and obtained the matrix equation T , (S 00 ) referred to a 12 � 12 matrix with only unknowns x 3 . This was an over-constrained linear system. The coefficient matrix (S 00 ) should be singular for the equation to have a solution, and the determinant of the coefficient matrix was a 16th degree polynomial in one variable regarding x 3 . The roots of the polynomial were corresponded to the solution to inverse kinematics. For improving the solving speed, the problem of solving 16th degree polynomial was simplified to find the matrix eigenvalues and eigenvectors. (S 00 ) can be written as Ax 3 2 þ Bx 3 þ C, and the Eq (12) can be written as A, B, C were a known 12 � 12 matrix, the Eq (13) can be written as The matrix eigenvalues of M were x 3 , and the eigenvectors were V ¼ v In this paper, the elimination method is optimized and supplemented in three points.
1. The paper presented the specific process of Eq 6 obtained by eliminating Eq 14 by e 1 e e 14 .
First, the expression s 1 , c 1 obtained from two equations e 3 , e 6 was substituted into e 7 , e 8 , e 9 , e 14 , getting 4 equations. Two equations would be obtained by simplifying the e 1 , e 2 , e 4 , e 5 , e 10 , e 11 , e 12 , e 13 via Eqs (15)- (16). There were no specific simplified formulas in other literatures. Where In other documents, the final expression of the matrix A, B, C was not deduced, and the solution required to be made in several steps. For improving the solving speed, the paper adopted complex coefficient extraction and simplification to derive the final expression of matrix A, B, C. Maximized the simplification of the whole solution process for achieving one-step solution of the matrix M. This was also the greatest contribution of the paper to eliminate optimization. The expression of A, B, C is shown in the appendix and refer to the uploaded matlab code. Refer to S1 Table. https://github.com/Wangxiaoqi1031/robot 3. In this paper, the solution formula in the simplest form of other joint angles θ 1 , θ 2 , θ 4 , θ 5 , θ 6 was introduced to integrate the whole method. Through a lot of experiments, taking the item with a large absolute value in the M matrix eigenvector as ratio, to get more accurate x 4 , x 5 . Via analysis, the maximum value may be v 1 , v 3 , v 10 , v 12 , and the formula θ 4 , θ 5 obtained was shown as Eq (17).
Where J + is the pseudo inverse of the robot Jacobian matrix J. This method gives the best possible solutions that minimize joint velocity norm k _ yk 2 and the end-effector tracking error k _ X À J _ yk :2 . However, the method is unstable near singularities. SVD decomposing the Jacobian matrix to get S is a diagonal matrix formed by the singular values of J.which are arranged in descending order, s 1 � s 2 � � � � s 6 . Hence, the joint speed could be computed as follows.
While a robot approaching a singular configuration, the smallest singular value reaches 0 and it causes infinite joint velocities as implied in (22). Therefore, damped least squares (DLS) introduces a damping factor, balances the precision of the solution and the increase in joint velocities near singularities by minimizing kJ _ y À _ Xk The joint speed could be computed as follows.
Multiple methods are generated from damping least square method and selecting strategies based on different damping factors. A constant damping factor was proposed in DLS method. The relative relation between the end-effecter position and the target position was considered in SDLS method, and a piecewise function was adopted to determine the damping factor. Gauss damping least square (GDLS) method refers to determine the damping factors with the following functions, based on the Gaussian distribution of the damping factor.
Where λ max is the maximum value of the damping factor and ε is a scalar quantity indicating region of singularities and σ i is singular value. GDLS is more effective than other methods that select constant and piecewise function damping factors. On one hand, each singular value corresponds to a damping factor, so that damping only act on singular vectors, to prevent the introduction of unnecessary damping and reduce the errors. On the other hand, when the robot is in close proximity to singular configuration, the method allows to make the joint velocity continuously change from undamped to damped, for more stable motion. The λ max , ε of the robot configuration in this paper was determined based on the genetic algorithm in the literature [8] and obtained λ max = 0.071, ε = 0.051. The Jacobian matrix pseudo-inverse method J T , DLS,SDLS and GDLS were used to analyze the inverse kinematics when close to the singular configuration and compare the performance of the four methods, as shown in Fig 4.  Fig 4 showed that when the robot approaches a singular configuration, compared with the Jacobian matrix pseudo-inverse method, the damped least square method prevented oscillation and effectively controlled the joint velocity, making the more stable motion. SDLS had similar effects as GDLS had, which reduced unnecessary damping and errors, compared with DLS.
DLS series methods including GDLS can be used to solve the IK problem of robot, and can ensure the stability of motion by controlling the joint speed. However, hundreds of iterations are needed to obtain the solution meeting the accuracy requirement in the solving process, which cannot achieve the goal of high-speed to solve IK. In the iteration process, it is found that the speed of approaching the iteration convergence can be improved through multiplying the iterative increment dθ by a coefficient k, achieving faster iteration convergence. The change of the number of iterations with k can be usually obtained as shown in Fig 5 in the process of solving with GDLS method. Therefore, we introduce the concepts of enhanced step length coefficient and optimal enhanced step length coefficient.
Definition: During the iteration, change θ cur = θ cur + dθ into θ cur = θ cur + kdθ, to make the iteration converge faster, and k is called the enhanced step length coefficient. The number of iterations reached the least, when k was assigned a value, and the value was defined as the optimal enhanced step length coefficient k_optimal.
To study the distribution of k_optimal and the range of the minimum number of iterations. We selected randomly 10 5 points and controlled the initial and target position within a certain

PLOS ONE
The optimized algorithm for inverse kinematics of two painting robots with non-spherical wrist range. The probability distribution of k_optimal and the minimum number of iterations was obtained and shown in Fig 6. Distribution statistics data of k_optimal were shown in Table 4. Fig 6 and Table 4 showed that the k_optimal was mainly distributed between [30,90] and centrally distributed between [55,60). If proper k was selected, the number of iterations can be controlled within 20 times.
To study the influencing factors of k_optimal, the iterative process was analyzed and the supposed possible factors to be: The difference value dif between current position T cur and terminal position T end , iterative convergence error E rr , initial iteration value θ i and singularity. Via the analysis on the singularity in section 3.1, the minimum singular value σ m was adopted to represent the singularity. When other conditions remained unchanged, we got and the variations of k_optimal and the minimum number of iterations with dif, E rr , θ i , σ m were shown in Fig 7.  From Fig 7, when E rr � 0.01, there was almost no influences of E rr on k_optimal. When dif changed over a larger scale, it did affect k_optimal. θ 3 and θ 5 had greater influence on k_optimal when the robot was far away from the singular configuration. When the robot approached the singular configuration, θ 2 , θ 3 , θ 4 and θ 5 had greater influence on k_optimal. After repeated experiments, we get that not every group of data could produce the same change rule and curve, while this section was intended to determine the influence factors of k_optimal, so the curve shape can be neglected. The singularity had no direct effect on k_optimal. Therefore, conclusion comes that the initial iteration value θ i has a great impact on k_optimal.
Based on the analysis on the influencing factors of k_optimal. The initial iteration value θ i was used as input and k_optimal as the output, and the k_optimal of each iteration was determined by the classification and regression models of machine learning. We selected 10 5 random positions as research data and changed the initial iteration value θ i to get k_optimal when other conditions remained unchanged.

• Classification method
According to the distribution of k_optimal in Table VI,   iterations changes with k continuously and gradually decreases when k approaches k_optimal. Therefore, a slight difference between the final selected k value and k_optimal is acceptable, and the number of iterations will be close to the minimum. The initial iteration value θ i was used as input and the classification of k_optimal used as the output, and the classification models of machine learning were used for classifying. The performance was evaluated by OA (overall accuracy), ECA (each class accuracy), AA (average accuracy) and TT (training time). Seven models with better performance were selected from various models for comparison, as in Table 5. Table 5 showed that the overall and the average accuracy of the DT model were the highest. The ECA obtained from DT and GBDT were relatively high. NB model featured the shortest training time. On the whole, DT model had the best classification result. After classification, in the iteration process, the median of each range is taken for the value of k_optimal of each class, which is k_optimal = 42.5,57.5,62.5,77.5.

• Regression method
The initial iteration value θ i was used as input and the k_optimal was used as the output, and the regression models of machine learning were used for analysis. The paper used three commonly used statistical indicators including root mean square error (RMSE), coefficient of  determination (R 2 ) and mean absolute error (MAE) to evaluate the performance of the model.

RMSE ¼
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Three indicators were calculated according to Eqs (25)- (27). Five models with better performance were selected from various models for comparison, as shown in Table 6, of which four with better results were selected and 1000 groups of data were taken at random to generate the regression fitting error curve, as shown in Fig 8. Table 6 showed that the DT model featured the best results for both training and testing samples. It could also be seen intuitively from

PLOS ONE
The optimized algorithm for inverse kinematics of two painting robots with non-spherical wrist According to the Eqs (2)- (7), we can get θ 1~θ6 . All the results that meet the accuracy requirements are as shown in Table 7.
For similar configurations, the formulas are universal.  (17)- (19). All solutions that meet the accuracy requirements, as shown in Table 8. For analyzing the algorithm performance, we selected 10 6 positions of the robot in the reachable workspace at random, to find a solution by the above optimized elimination method. And the testing environment was (Intel(R) Xeon CPU2.8GHz, python). The statistical results showed that the average solving speed was reduced from 11-13ms before optimization to 1.48ms, and the solving speed was reduced to 20% of the original, the average accuracy of solution reached 2.3205e −13 . During solving by elimination, the proportion of the points that cannot be solved accurately since the A matrix was singular or the robot approached the singular configuration was 0.086%, which can be solved by the following solutions.

Examples of oblique non-spherical wrist robot
1. In the case of A is singular, finding inverse is changed to find pseudo-inverse.
2. According to the singularity analysis in Section 3.1, there is no singular regions in trajectory planning.
3. Points that cannot be solved accurately can be processed by iteration method.
Examples and properties of optimized Gaussian damped least square method. For comparing the solving performance of the classification and regression, we selected 10 5 positions at random and adopted DT model to predict the k value separately, and solved it in accordance with the above steps. The resulting statistics, as shown in Table 9. Table 9 gave that the regression method featured less average number of iterations. The two methods had almost the same average prediction time, and the average prediction time was negligible compared with the solution time. The regression method had a lower proportion of iterations over 10 times and over 20 times, less than 0.1%. On the whole, the regression method had better effects than classification.
For comparing the solving performance of the methods before and after optimization, the Gauss damped least square method and the optimized method in the paper were used for solving separately. A complete experimental process was: we randomly selected 10 5 data points as the original data set, each data point contains the initial iteration value and k_optimal. The collected data set was used to train a DecisionTree model through supervised learning to obtain a model that can be used to predict k_optimal. Then 10 4 data points were randomly selected in the robot workspace, that is, the initial iteration values were randomly selected, and the IK was solved by the GDLS method before and after optimization. Then we counted the number of iterations and the solution time. The above experiment was repeated 10 times and the statistics were obtained. The best, worst, and average cases in the statistical repeat experiments were shown in Table 10. Table 10 showed that the average iteration times of the optimized method was obviously reduced, and the average solution time was reduced to 5% of the original, the average solution time of optimized method was close to the solution time of the analytic solution of the general configuration. The optimized method in this paper was very effective for improving the solving speed. A simulation platform was developed, based on Qt integrated opengl, QML and C ++ language. It was superior in cross-platform, high flexibility and high reusability, etc. Combined with the interpolation algorithm, a trajectory was planned. The algorithms in the paper were adopted, the robot model moved the planned trajectory on a curved surface in the simulation platform, and the collecting data were analyzed, as shown in Fig 10. According to Fig 10, the planned trajectory of a robot model moving on a curved surface was accurate and the robot moved smoothly. The joint angle changed smoothly in the joint space, and the trajectory tracking errors were up to grade. The optimized method could ensure the stability and accuracy of the motion.
There was the experimental platform of oblique non-spherical wrist robot in the lab as shown in Fig 11. Same as the configuration of robot studied in the paper, the structure parameters were a 1 = 0.135, a 2 = 0.7, a 3 = 0.1, d 4 = 0.433, d 5 = 0.115, d 6 = 0.174, β = 60˚. The two methods optimized in this paper were encoded and implemented on the DSP end of the controller, to control the motion of the robot at high speed and smoothly. The effectiveness of the methods were verified.

Conclusions
For the orthogonal non-spherical wrist configuration, the simplest analytical solution was given in this paper. For the oblique non-spherical wrist configuration, elimination method and Gauss damped least squares method were selected to optimize the solution. For the elimination method, this paper mainly made three points optimization, gave the expression of equality simplification, deduced the expressions of final A, B and C matrix to improve the speed of solution, and gave the simplest expression of all joint variables to complement the integrity of elimination method. Elimination method was mainly time-consuming in solving

PLOS ONE
The optimized algorithm for inverse kinematics of two painting robots with non-spherical wrist eigenvalues and eigenvectors of M matrix. Increasing the speed of this step implies the total cost of the algorithm will be decreased. For the Gauss damped least squares method, the enhanced step length coefficient was proposed. By analyzing the influencing factors of the coefficient, the correlation between the coefficient and the initial iteration value was determined. The concept of optimal enhanced step length coefficient was introduced, and studied using the classification and regression methods of machine learning. By comparing the performance indicators, it was determined that the Decision Tree algorithm had better performance for the study of this problem. Follow-up research expects to obtain more accurate coefficients approaching the optimal enhanced step length coefficients through other methods, so that the number of iterations for each solution can be minimized.
Supporting information S1