Learning the Kinematics of a Manipulator Based on VQTAM

: The kinematics of a robotic manipulator is critical to the real-time performance and robustness of the robot control system. This paper proposes a surrogate model of inverse kinematics for the serial six-degree of freedom (6-DOF) robotic manipulator, based on its kinematics symmetry. Herein, the inverse kinematics model is derived via the training of the Vector-Quantified Temporal Associative Memory (VQTAM) network, which originates from Self-Organized Mapping (SOM). During the processes of training, testing, and estimating of this neural network, a priority K-means tree search algorithm is utilized, thus improving the training efficacy. Furthermore, Local Linear Regression (LLR), Local Weighted Linear Regression (LWR), and Local Linear Embedding (LLE) algorithms are, respectively, combined with VQTAM to obtain three improvement algorithms, all of which aim to further optimize the prediction accuracy of the networks for subsequent comparison and selection. To speed up the solving of the least squared equation, which is common among the three algorithms, Singular Value Decomposition (SVD) is introduced. Finally, data from forward kinematics, in the form of the exponential product of a motion screw, are obtained, and are used for the construction and validation of the VQTAM neural network. Our results show that the prediction effect of the LLE algorithm is better than others, and that the LLE algorithm is a potential surrogate model to estimate the output of inverse kinematics.


Introduction
The robotic manipulator has been commonly used in various fields. Robot control is critical for high-speed and high-precision robot motion, and it is based on the kinematics of robots. Kinematics includes two aspects: forward kinematics and inverse kinematics. Inverse kinematics describes the mapping from the state of effector to the state of actuator [1]. The key to kinematics is establishing the mapping relationship of the manipulator. For the serial robot, the input is the rotation angle of each joint, and the output is the pose of the end effector.
Generally, the inverse kinematics analysis of serial robots requires highly complex calculations, which affect the response speed and work efficiency of the robot controller. The same problem also exists in the forward kinematics of parallel robots. Three types of methods are typically used for the inverse kinematics of serial robots: geometric [2][3][4], algebraic [5][6][7], and iterative algorithms [8][9][10][11][12]. However, algebraic methods cannot always derive a closed-form solution. In addition, the use of geometric methods is limited by the precondition, and the initial position influences the convergence

Research Methodology
A robotic manipulator is a typical open-loop mechanism. Its forward kinematics can be achieved by coordinate transformation, D-H parameters, screw theory, and other mathematical tools. The result can be obtained quickly. However, the inverse solution of serial robots is very complex. The robotic manipulator is a nonlinear system. Generally, solving the inverse kinematics and dynamics accurately and conveniently is difficult.
Thus, the poses corresponding to different joint angles can be obtained using forward kinematics. The data obtained can be used to train the surrogate model for inverse kinematics based on its symmetry. In this paper, a VQTAM neural network is introduced to replace the complex calculation of inverse kinematics. The research content and structure of this paper are shown in Figure 1.
In Section 3, the forward kinematics is studied. Data for training are sampled from forward kinematics, meaning that each pose of the end effector corresponds to a group of joint angles.
In Section 4, the nonlinear system identification method based on VQTAM is studied, and the training method of the VQTAM neural network is established. In the process of VQTAM training and testing, it is critical to search activated neuron nodes. Moreover, the priority search K-means tree algorithm is introduced for this step.
In Section 5, the improved method of local linearization for the mapping process based on the VQTAM neural network is studied.

Kinematics of the Robotic Manipulator
In screw theory, any motion of a rigid body can be decomposed into rotational motion around a straight line and translational motion along the line. That is to say, the motion of a rigid body can be regarded as spiral motion.
Using screw theory has the following advantages in rigid body kinematics [5,30,31]: (1) The 6-DOF parameter is used to describe the pose relationship between two adjacent coordinate systems completely, thus avoiding the lack of completeness in the D-H model.
(2) The global coordinate system is used to describe the motion state of a rigid body, which overcomes the singularity of the D-H model. (3) The motion characteristics of rigid bodies can be clearly described from a global perspective, thus simplifying the analysis of complex mechanisms and avoiding the abstraction of mathematical symbols.
To sum up, the advantages of screw include a clear geometric concept, clear physical meaning, simple expression, and efficient algebraic operation in the kinematics of the robotic manipulator.
According to Chasles' theorem, the spiral motion of a rigid body can be expressed in the form of exponential coordinates of the screw. For the 6-DOF robotic manipulator (RRRRRR type) shown in Figure 2, the relative inertial coordinate system of each joint is established according to screw theory, as shown in Equation (1). The coordinates of the points on the axes of each joint are shown in Equation (2).
Then, the unit screw of each joint can be obtained as shown in Equation (3).
ωi can be mapped to the rigid body rotation transformation matrix to describe the rigid body rotation transformation in three-dimensional space, as shown in Equation (4). vi is the vector comprising the third, fourth, and fifth components of the screw, ξi.
ξi can be mapped to a rigid body transformation matrix, describing the rigid body transformation in three-dimensional space, as shown in Equation (5).
As a special Euclidean group SE(3), the rigid body transformation matrix in three-dimensional space can be regarded as a six-dimensional manifold, and the points on the manifold correspond to the rigid body transformation. SE(3) comprises all possible rigid body motion mappings. Then, i is the Lie algebra se(3) of SE (3). According to the isomorphism relation shown in Equation (6), the mapping relation between the Plucker coordinate form and matrix i can be defined and expressed using operators ∨ and ∧, as shown in Equations (7) and (8). 6 (3) ; 0 The helical motion of a rigid body g∈se(3) is expressed in the form of exponential coordinates of the motion screw, as shown in Equation (9) [30,31].
Combining the motion of each joint, the exponential product equation of the forward kinematics of the robot is obtained, as shown in Equation (10) [30,31]. 6  , , , , , θ θ θ θ θ θ = θ (10) gst(0) is the pose in the initial state, and gst(θ) is the pose of the manipulator, where θ is the angle vector of each joint relative to the initial state. Given the joint screw and initial pose, the forward kinematics solution of the robotic manipulator can be obtained.
The joint screw parameters of the robotic manipulator are given in Table 1. The pose shown in Figure 1 is determined to be the initial pose. Poses corresponding to different joint rotation angles can be obtained by Equation (10). Some data for computation validation are shown in Table 2 and Figure 3, which were drawn using the robotics MATLAB toolbox by Peter Corke. In Table 2, (x,y,z) represents the position coordinate of the end actuator, and (u,v,w) is the attitude in Euler angle in radians.   Table 2. Table 2. Poses corresponding to different joint rotation angles.

System Identification of a Nonlinear Dynamical System Using VQTAM
The inverse kinematics of the robotic manipulator aims at solving the corresponding joint rotation once given the pose and initial pose of the manipulator. The joint rotation angle is regarded as output Q and the pose as input p. The inverse kinematics problem of the robot can be transformed into the identification problem of a nonlinear dynamic system. The time-discrete difference equation is established as shown in Equation (11) [27,29]. Q(t) = θ(t) is the joint rotation angle variable of the manipulator. p(t) = [x(t),y(t),z(t);u(t),v(t),w(t)] is the pose of the manipulator at t time. It is expressed by the position coordinates combined with the Euler angle. f(•) represents a nonlinear function reflecting the characteristics of the system. The output Q at t + 1 is determined by the previous nq outputs and np inputs.
The purpose of inverse kinematics is to find the mapping function f(•) between the input and output. As a complex nonlinear problem, on the premise of ensuring the accuracy, simplifying the model can improve efficiency and ensure real-time performance in engineering applications. (13) VQTAM is a variant of SOM, which can be used to realize nonlinear mapping. The concept of nonlinear manifold embedding can be used to establish the nonlinear mapping. SOM is realized through competition and cooperation among feature neurons. In the mapping process, different neurons are activated for different inputs, and the activated neurons affect the output together with the neurons in their neighborhood. Therefore, time series data can be constructed as training samples of VQTAM.
VQTAM adopts training samples to build a topological structure embedded in multidimensional data space, which can construct the mapping relationship between multidimensional input vectors and multidimensional output vectors. VQTAM consists of three layers: input space ω in , output space ω out , and lattice space. The establishment of the mapping process is shown in Figure 4. The definitions of X in and X out are shown in Equations (12) and (13), respectively. The input space and output space are composed of the weight vectors ωi in and ωi out , respectively, where i is the index of the neuron in lattice space that reflects the location of the neuron in the lattice topology structure. Weight vectors ωi in and ωi out have mapping relationships with neurons in lattice space. ωi in and ωi out have the same dimensions as X in and X out , respectively. In the mapping process of VQTAM, the nearest weight vector ωi* in to X in is found in the input space, which corresponds to the activated neuron in the lattice space. In this process, the Euclidean distance is used to measure the distance between two vectors. The definition of i* is shown in Equation (14), where A is the collection of all neuron indexes in lattice space [27].
ωi* out is obtained according to the neuron index i*, and the output X out is estimated, as shown in Equation (15).

Learning Strategy of VQTAM
The main parameters of VQTAM are neuron weight vectors ωi in and ωi out , and lattice topology. The initial lattice topology structure can be determined as a hyperparameter before training, so the learning strategy mainly aims at updating the weight vectors ωi in and ωi out in the training process.
There is competition and cooperation among neurons in SOM. Different neurons and their neighborhoods are activated by inputting to participate in the learning process, and the weight vectors ωi in and ωi out are updated. The index search of the active neurons in the lattice space is consistent with Equation (14).
To improve the convergence rate in the learning process, the influence range parameters σ(t) of the activation neuron and the learning speed α(t) decrease exponentially with the learning epoch, as shown in Equations (16) and (17).
where t is the current learning epoch, T is the total learning epoch.
Complete the updates of ωi in and ωi out , respectively, as shown in Equations (19) and (20).
To sum up, the VQTAM algorithm flow is as follows.

Algorithm 1: VQTAM Algorithm:
Begin (training part) 1 Input: X in and X out in training set 2 Search activating neuron according to X in * argmin{ } Update the weight vector of the neuron Continue execution until termination conditions are met (testing part) 5 Input: Xtest in in testing set 6 Search activating neuron according to

Searching the Activated Neuron by the Priority Search K-Means Tree Algorithm
During the training, testing, and estimation of VQTAM, the weight vector ωi* in to which X in is nearest in the input space is searched, as shown in Equation (14). The search algorithm applied has a great impact on the efficiency of VQTAM. The global traversal method needs to find all the Euclidean distance from the weight vector ωi in to the input X in , and search the minimum value through traversal, which is inefficient. The K-Dtree algorithm can reduce the time complexity of searching, but for d-dimensional data, the time complexity of the K-Dtree data structure search is O (n d/(d-1) ). For a high-dimensional data search, especially when the data dimension d > 20, the K-Dtree data structure is inefficient. ωi in is 30-dimensional. The searching efficiency is severely restricted. The search time complexity of the priority search k-means tree algorithm (PSKMT) is O(Ld logn/logK), where n is the amount of data in the data set, K is the number of nearest neighbors to be searched, and L is the maximum number of data retrieved in the search process. The priority search k-means tree algorithm is suitable for searching high-dimensional data.
The PSKMT algorithm divides the space into B different regions using the k-means clustering algorithm. Then, the points in each region are partitioned by the same operation until the number of data points in the region is no greater than K. In the process of searching, the idea of "divide and conquer" is used to find the nearest point in a smaller area and reduce the amount of data to be searched, so as to improve the efficiency.
The algorithm used for establishing a priority search K-means tree data structure is as follows. if P = Pnew, Pnew is the non-leaf node, and terminate the loop 9 These processes are executed on the sub-regions C until all leaf nodes are created 10 Output: the entire K-means tree data structure In the establishment of the K-means tree, the algorithm Calg is used to select B points from data set D that can be chosen from the random selection, Gonzale, and K-Means++ algorithms. The specific algorithm has little effect on data structure establishment. Generally, the random algorithm is recommended [16].
In the built K-means tree, leaf nodes are the data points in the original data set, and non-leaf nodes are the central points of the regions after segmentation. The search process can only be executed in a few areas, thus avoiding the global search of data sets. The K-means tree data structure is searched from the root node. Sub-nodes are sorted according to the distance between the clustering center and the query data points. The nearest sub-nodes are searched in advance. The priority K-means tree search algorithm is as follows.

Algorithm 3: PSKMT Algorithm [32]:
1 Input: K-means tree, query data Q (X in ), the nearest neighbor number K, the maximum number of search data L 2 Initialize a stack and place the root node in the stack 3 Starts loop. The condition for termination is that the stack is not empty and the amount of retrieved data |P| is less than L: 4 If the top node of the stack is a leaf node, add it to the retrieved data array P 5 Otherwise, the top node goes out of the stack, reads the sub-nodes, sorts them according to the distance between the clustering center and the query data points, and pushes them into the stack. 6 Loop terminated 7 Output: K data nearest to query data Q in retrieval data array P.

VQTAM Local Linear Improvement Algorithms
The VQTAM algorithm is used to quantify the input space ω in and output space ω out . The neurons in ω in and ω out correspond to each other through mapping relations. The input X in is approximated to the nearest neuron ωi* in in estimation. The accuracy of the estimation results can be guaranteed when the number of neurons is large enough. However, the increase in number is followed by an increase in network size and a decrease in computing efficiency. Thus, local linearization of the activated node is used, which can balance the number of neurons and prediction accuracy. Based on local linearization, three improved algorithms for VQTAM are proposed: Local Linear Regression (LLR), Local Weighted Linear Regression (LWR), and LLE. The inverse kinematics can be further optimized, ensuring the computational efficiency.
Local linearization is performed in VQTAM networks. The K-means tree algorithm is used to search the nearest n data points ωi*n in and the corresponding output data ωi*n out are mapped, as shown in Figure 5.

Improvement Algorithm of VQTAM with LLR
It is assumed that the relationship between the local input and output can be expressed linearly, so LLR can be performed as shown in Equation (21). = YB Z (21) Y is the nearest neighborhood ωi*n in of X in , and Z is the corresponding data point ωi*n out , as shown in Equations (22) and (23).
The linear relationship between the local input and output can be obtained by solving B. Equation (21) is an overdetermined equation, which is generally treated as a least squares problem and solved by a generalized inverse matrix, as shown in Equation (24).
Y T Y is an n-dimensional square matrix, and the calculation of the inverse operation is huge when the dimension increases. Meanwhile, Y T Y may be a singular matrix or ill-conditioned matrix, so finding the inverse matrix directly is difficult. Singular Value Decomposition (SVD) can be used to solve the least squares problem. Y can be decomposed as shown in Equation (25).
The column vectors of U are the eigenvectors of YY T ; the column vectors of V T are the eigenvectors of Y T Y.
Σ is a diagonal matrix whose value is the singular value υ of matrix Y, i.e., the eigenvalue of YY T is λ = υ 2 , and U and V T are the orthogonal matrix.
U can be disassembled as [Un, Um-n]. For the least squares problem as shown in Equation (21), the solution obtained by SVD is shown in Equation (27).
The input X in is used to estimate the output, as shown in Equation (28).

Improvement Algorithm of VQTAM with LWR
Using the LLR algorithm, the impact of all neighbor points on the prediction results is consistent. Considering the different influences of neighbor points, a distance weight is added in the regression process. A local weight linear regression (LWR) is adopted. The cost function is defined as shown in Equation (29).
The weight wk in Equation (29) is determined by the Gauss neighborhood function, as shown in Equation (30), where τ is the bandwidth parameter.
The analytical solution of the cost function minimization is equivalent to the solution of the overdetermined equation shown in Equation (31). (31) where W is an n-dimensional diagonal matrix, and the diagonal element is wk.

WYφ WZ
The SVD for matrix WY is also used to solve the overdetermined equation.

Improvement Algorithm of VQTAM with LLE
It is assumed that an input datum can be represented by a linear combination of several samples in its neighborhood; that is, the input X in to be predicted is represented by a linear combination of n data points in its neighborhood in the input space ωi*n in .
The output has the same linear combination.
The least squares problem, as shown in Equation (35), can also be solved by SVD.

Standard VQTAM Network Test Results
According to the parameters of the robot shown in Table 1, the joint rotation angles for the sample set are generated as shown in Equations (38)-(43) [26,28]. The generated data can effectively cover the whole workspace of the manipulator. Taking time t∈[0,1000] and step t = 0.1, 10,001 joint angles Q(t) = [θ1,θ2,θ3,θ4,θ5,θ6] are calculated as the output of the sample set. The poses p(t) = [x,y,z;u,v,w] of the manipulator are calculated using the exponential coordinates of the motion screw, which are used as the input of the sample set. The VQTAM network is trained with the sample set X in and X out . In each epoch, 1000 pieces of data are selected randomly from the sample set, 80% of which are used as the training set and 20% as the test set. Meanwhile, to test the robustness of the algorithm, random errors are added to the test set data according to N(0,1) distribution.
The hyperparameter setting of the VQTAM network training is shown in Table 3. A VQTAM network with dimensions 60 × 60 (Mx × My) is trained based on the setting of hyperparameters, and the training effect is tested by the test set. Root Mean Squared Error (RMSE), R-squared (R 2 ), and Relative Maximum Absolute Error (RMAE) are used to evaluate the prediction accuracy of VQTAM for inverse kinematics, as shown in Table 4. The value of R 2 ranges from 0 to 1. The closer it is to 1, the better the effect of regression fitting is, which means that the learning VQTAM network has higher accuracy as an approximate model. Generally, R 2 > 0.95 can be applied in engineering. For VQTAM networks with dimensions 60 × 60, the R 2 values after learning are all above 0.98. This shows that this network has value for engineering applications. RMSE and RMAE reflect the prediction error of the VQTAM network, and their small values also indicate the prediction accuracy of the network. The convergence result in the neural network training process is shown in Figure 6.
It can be seen from Figure 6 that the mean square R 2 of the VQTAM neural network tested by the test set in the training process gradually increases and converges near 1. This proves the effectiveness of the VQTAM neural network training algorithm. The consistency between the actual data and the estimated data can be displayed more intuitively through box diagrams, as shown in Figure 7. It can be seen that the box diagrams of actual values and predicted values are very similar, which verifies the excellent prediction effect of the VQTAM network with dimensions 60 × 60.

VQTAM Local Linear Improvement Algorithms
The output of the inverse kinematics problem is estimated using the VQTAM local linear improvement algorithms. The priority K-means tree search algorithm is used to search the neighbor data of the input. The influence of parameter k on the prediction accuracy is analyzed, taking a VQTAM network with dimensions 35 × 35 as an example, as shown in Figure 8. As can be seen from Figure 8, with the increase in k, R 2 increases continuously and approaches 1, which shows an increase in prediction accuracy. When k = 6, the R 2 for all variables is close to 1.
As can be seen from Figure 9, when k = 5, the R 2 for all variables is more than 0.95. When k = 6, R 2 is close to 1. Comparing Figures 8 and 9, the overall prediction accuracy of the LWR algorithm is higher than that of the LLR algorithm.
As can be seen from Figure 8, when k = 5, the R 2 for all variables is greater than 0.99. Comparing Figures 8-10, the overall prediction accuracy of the improvement algorithm of VQTAM with LLE is the highest.  In conclusion, when the number of neurons is small, the prediction accuracy of the network can be improved by combining local linear algorithms.

Overall Examples and Test Results
According to the VQTAM network hyperparameters shown in Table 3, the prediction accuracy of the four algorithms under different network dimensions (5 × 5, 10 × 10, 15 × 15, 20 × 20, 25 × 25, 30 × 30, 35 × 35, 40 × 40) is investigated, as shown in Figure 9. The number of the nearest nodes of the local linear algorithm is set to k = 6.
As can be seen from Figure 11, with the increase in the number of neurons, the prediction accuracy of the four algorithms significantly increases. The overall prediction accuracy of the improvement algorithm of VQTAM with LLE is the highest. The local linearization of VQTAM achieved remarkable results. Using the improved VQTAM algorithm, the estimation accuracy level of a network with dimensions of 35 × 35 is close to that of the standard VQTAM with dimensions of 60 × 60. Taking the improvement algorithm of VQTAM with LLE as an example, the output of the inverse kinematics, i.e., the rotation angles of six joints, is estimated, as shown in Figure 12. The estimated rotation angles of joints are compared with actual angles. It is shown that the VQTAM neural network can effectively estimate the output.  Substituting the joint rotation angle into the forward kinematics equation for calculation, the error of VQTAM can be obtained by comparing the result with the actual pose of the robot. The maximum error in the tool space is 3.9686 × 10 −4 mm, and the average error is 4.6857 × 10 −5 mm. This precision has met the control requirements of general robots in industrial applications.
The multilayer perceptron (MLP) can also provide the estimated results for the inverse kinematics. Thus, an MLP NARX network is trained for the comparison of prediction accuracy with VQTAM and its improvement algorithms. The MLP network has a hidden layer, and the number of hidden neurons is 3600. The sample set p(t), Q(t) for training MLP is the same as that for VQTAM. The Levenberg-Marquardt algorithm is chosen as the training algorithm. The comparison of prediction accuracy is shown in Table 5. The prediction accuracy of the standard VQTAM network is significantly higher than that of the MLP neuron network. Meanwhile, the prediction accuracy of the improvement algorithm of VQTAM also has better performance. The VQTAM is less sensitive to weight initialization than the MLP. Training multiple times will generate different results due to different initial conditions. The variation in the range of results of VQTAM is smaller. Compared with the results of a one-dimensional VQTAM in a previous study [27], the increase in dimensions improves the performance of the neural network.

Conclusions
The kinematics of the robotic manipulator is studied, and the forward kinematics solution method in the form of the exponential product of the screw is derived. Given the input angle of each joint of the robot, this method can be used to calculate the pose of the manipulator conveniently and quickly. It can also be used as the computational basis of a VQTAM network to establish sample sets.
A fast inverse kinematics mapping method for a robotic manipulator is proposed by using VQTAM for identification of nonlinear dynamic systems. The VQTAM algorithm is constructed to train and test the network. Based on the priority K-means tree search algorithm, the training, testing, and estimating processes of the VQTAM network are improved to enhance the search efficiency of nearest neighbors.
According to local neuron activation, the VQTAM network is improved based on local linearization. To further optimize the prediction accuracy of the network and reduce the dimensions of the network in application, VQTAM local linear improvement algorithms are proposed.
In this study, a sample data set for VQTAM training is constructed and the VQTAM algorithm and its improved algorithms are tested. The test results show that the increase in network dimensions can improve the prediction accuracy of VQTAM, but the computational efficiency is affected. By using VQTAM local linearization improvement algorithms, the estimation accuracy can be optimized for a low-dimensional network. When the search range is k = 6, the above algorithm can show good prediction accuracy, and the prediction effect of the LLE algorithm is the best among the three algorithms.