A Novel Solution to the Inverse Kinematics Problem of General 7R Robots

Task and Motion Planning for a 7-DOF serial robot containing only revolute joints needs a fast and numerically stable inverse kinematics solution. The analytical methods are real-time, but they are robot dependent. Although the numerical methods based on inverse differential kinematics can solve the inverse kinematics problem (IKP) of 7R robots with arbitrary geometric parameters, their numerical stability is insufficient. Moreover, it is time-consuming for them to converge to a single solution that depends on the initial guess. The main innovation of this article is the development of a practical method for solving the IKP of general 7R robots. The proposed method can find multiple solutions and is also faster than the numerical method. It combines symbolic preprocessing, Sylvester dialytic elimination, and eigendecomposition. First, according to the Denavit-Hartenberg convention of geometric representation, a general kinematics model of a 7R manipulator is established. Secondly, the joint variables are separated by symbolic preprocessing, and then two symbolic matrices are obtained. After numerical substitution, a square matrix of polynomial coefficients is computed through Gaussian and Sylvester dialytic elimination. Finally, utilizing the eigendecomposition of the matrix restructured from the square matrix, the values of some variables are obtained, and then the remaining variables are solved by back-substitution. The simulation and experimental data are analyzed. The comparison results verify that the proposed method applies to solving the IKP of general 7R robots and is almost real-time, effective, and numerically stable.


I. INTRODUCTION
As far as one pose in the three-dimensional workspace for the end-effector of a 7R robot is concerned, an infinite number of inverse kinematics (IK) solutions can make the endeffector reach it. This feature affords a robot with avoiding workspace obstacles [1], refraining from singular configurations [2], overcoming joint limits [3], and improving dexterity [4]. Solving the inverse kinematics problem (IKP) of redundant robots is a nonlinear problem of mapping from the robot workspace to the joint space, seeking the numerical or analytical solutions of a system of nonlinear transcendental equations.
The associate editor coordinating the review of this manuscript and approving it for publication was Yangmin Li .

A. RELATED WORK
At present, the methods to solve the IKP for redundant degrees of freedom (DOF) robots are mainly numerical, such as pseudo-inverse of Jacobian, extended Jacobian, augmented Jacobian, gradient projection, weighted least-norm, damped least square, weighted least square, sequential quadratic programming, etc.
To resolve motion rate control on the IK of multibody mechanisms, Liegéois [5] proposed the pseudo-inverse method. However, this method has the weakness of causing extraneous motion in null space. Using a composite Jacobian, Oh et al. [6] developed the iterative numerical method of the Newton-Raphson procedure. It can be applied to solve the IKP for redundant or non-redundant manipulators without having to derive beforehand a closed-form solution. Using a modified Newton-Raphson technique for solving a system of nonlinear equations, Goldenberg et al. [7] presented VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ a complete generalized solution to the IK of robots with any arbitrary number of DOF. The method based on the Newton-Raphson algorithm has the disadvantage of uncertain iteration numbers and computing the pseudo-inverse of the Jacobian matrix for singular configurations. To solve the singularity problem, Nakamura and Hanafusa [8] introduced the singularity robust inverse of the Jacobian matrix as an alternative to the pseudoinverse of the Jacobian matrix. Baillieul [9] imposed a constraint task to the mapping between joint displacements and poses of the end-effector and presented the extended Jacobian technique, which was viewed as a generalized inverse technique. By projecting the constraint Jacobian onto the null space of the end-effector Jacobian, Chiacchio et al. [10] developed a closed-loop IK scheme to solve the IK for redundant manipulators when additional constraints were imposed. Employing numerical estimation techniques, Tevatia and Schaal [11] introduced a computationally efficient version of the extended Jacobian that had better performance than the original version. The simplified extended Jacobian algorithm was capable of reducing the irrelevant null-space motion for certain optimization criteria. To solve the approximation problem of the IK algorithms for redundant manipulators, Ratajczak [12] introduced the approximation of dynamically consistent Jacobian. Using the Cholesky decomposition and the Ritz method, he proposed an optimal extended Jacobian IK algorithm. To keep robots far from singularities, Simas and Gregorio [13] proposed the algorithm of adaptive coefficients in constraint functions to obtain a good conditioning index of the extended Jacobian.
The sequential quadratic programming (SQP) method efficiently solves a series of related nonlinear optimization problems. As an iterative algorithm for nonlinear optimization, the SQP method is adequate for solving the problems when the objective function with constraints is twice consecutively differentiable. Leger and Angeles [14]proposed a method to solve the forgoing redundancy of 5-DOF tasks based on SQP. Jiang et al. [15] used the SQP method to improve the efficiency and trajectory tracking accuracy of a robot and eliminate its jerks. After parameterizing the predefined geometric path, the trajectory generation problem can be transformed into a nonlinear optimization problem with a single state variable. Lyu et al. [16] applied the SQP method to generate a time-optimal and energy-efficient trajectory for manipulators. By the SQP method, Taniai and Naniwa [17] computed the optimal arm angles within the constraints of the joint angles. Huang and Chen [18] used the SQP method to obtain the optimal obstacle avoidance pose of a 7-axis robotic manipulator.
Unlike redundancy resolution schemes based on differential kinematics which lack repeatability, the augmented Jacobian method is always cyclic. Using the systematically enlarged direct kinematics, Sciavicco and Siciliano [19] proposed a dynamic solution technique to solve the IKP for constrained redundant manipulators. The kinematics of a robotic manipulator was appropriately augmented by including the constraints of obstacle avoidance and limited joint range, called task-space augmentation. It introduced a constraint task to be fulfilled along with the end-effector task. Then, an augmented Jacobian matrix was set up whose inverse gave the sought joint velocity solution. The concept of taskspace augmentation was also independently proposed by Egeland [20], to avoid singularities and the loss of DOF for generating a position reference of the positioning part online. Incorporating the kinematic constraints and the end-effector motion, Seraji [21] revisited the concept and developed a simple framework for the configuration control of redundant manipulators. By parameterizing the arm angle and constructing its Jacobian matrix, Cao et al. [22] proposed an iteration algorithm based on augmented Jacobian to solve the IKP for 7-DOF manipulators with arbitrary offsets.
Using the gradient projection optimization, Dubey et al. [23] presented a computationally efficient kinematic control scheme for a 7-DOF robot with a spherical wrist. This scheme was modified for general 7-DOF robots [24] and did not need to determine the pseudo-inverse of the Jacobian. Zghal et al. [25] extended the scheme presented by Dubey et al. and developed an efficient gradient projection optimization scheme for robotic manipulators with multiple degrees of redundancy. In order to avoid joint limits for redundant joint manipulators, the gradient projection method provides the optimal direction for the joint velocity vector within the null space. Nevertheless, its magnitude is not unique and is adjusted by a scalar coefficient chosen by trial and error. Using a weighted least-norm solution, Chan and Dubey [26] proposed a manipulation scheme that automatically chose an appropriate magnitude of the selfmotion all through the workspace. Unlike the gradient projection method, this scheme guarantes joint limit avoidance and minimizes unnecessary self-motion.
The IK formulations based on the Jacobian matrix are inefficient and fail near kinematic singularities. To overcome kinematic singularities, Wampler [27] reformulated the exact inverse problem as a damped least-squares problem and presented an efficient vector formulation of the IKP in terms of a partial velocity matrix. Using the Gauss-Newton model of the direct kinematics function with the Levenberg-Marquardt iteration, Deo and Walker [28] proposed an adaptive nonlinear least-squares algorithm to solve the IKP for general robotic manipulators. The algorithm was insensitive to the reachability of the desired position, and convergence was ensured even if the goal position corresponds to a singular configuration. Adopting the damped least-squares technique to provide numerical robustness of the solution in the neighborhood of kinematic singularities, Caccavale et al. [29] utilized closed-loop IK algorithms and a two-stage algorithm to realize real-time kinematic control on a seven-joint manipulator. Buss and Kim [30] introduced the selectively damped least squares method for the IK of multi-body with multiple end-effectors. It adjusted the damping factor separately for each singular vector of the Jacobian singular value decomposition based on the difficulty of reaching the target positions.
It had advantages in convergence with fewer iterations and in not requiring interim damping constants. Le et al. [31] proposed two strategies based on the damped least square method to effectively deal with singularities and minimize error due to the introduction of damping. Unlike other work in which the same damping factor was used for all singular vectors, the genetic algorithm was implemented to search for a different damping coefficient for each singular vector based on the corresponding singular value of the Jacobian. To solve the problem of multi-solution caused by the redundancy of a 7-DOF manipulator, Yang et al. [32] introduced the rule of best compliance based on the weighted least square method. Applying a genetic algorithm to search for all global optimum solutions made the multi-solution a mono-one. Sugihara [33] proposed a robust numerical IK solution based on the Levenberg-Marquardt algorithm. The solution employed the squared norm of residual of the original equation as the damping factor and avoided heavily computing the manipulability and the minimum singular value of the coefficient matrix.
Since redundant robots are highly nonlinear, coupled multi-variable systems, it is challenging to obtain the IK solution in an analytic way that can obtain all the feasible solutions in the global joint space. Analytical methods are robot dependent and only applicable to redundant robots with simple configurations or unique geometric properties. Mostly, analytical solution algorithms employ screw theory and geometric features analysis. Hemami [34] presented solutions to four different classes of robotic manipulators obtained by locking the redundant joints. When the extra joints are randomly locked at arbitrary angles, the resultant will be non-redundant manipulators of various structures for each class. For the 7-DOF robotic manipulators of which all links length is zero, Shimizu et al. [35] derived a closedform inverse kinematic solution based on a parameterization method and developed an analytical method for computing feasible solutions under the joint limits. Singh and Claassens [36] provided an analytical solution to the IKP for the Barrett WAM manipulator and addressed an analytical method to identify a set of feasible poses for some joint-angle constraints. Wang et al. [37] presented an analytic IK solution for a 7-DOF space manipulator with an S-R-S structure and provided an optimization approach to get a near-optimal IK solution.
In order to realize the efficient search of IK solutions for redundant DOF robots, intelligent algorithms are widely used, for example, genetic algorithms, artificial neural networks, and particle swarm optimization algorithms. Parker et al. [38] used genetic algorithms to position a redundant robot at a target location by minimizing the maximum joint displacement in a point-to-point positioning task. Considering the point-to-point motion of redundant manipulators working in environments with obstacles, Nearchou and Aspragathos [39] took the problem as a constrained optimization problem and solved the IKP using a method based on a simple genetic algorithm. Although the method could be successfully applied to any redundant or non-redundant manipulators with obstacles in their workspace, it produced a rather sizeable positional error of the end-effector when the environment was cluttered with obstacles or when some obstacles were too close to the desired location. To solve this problem, Nearchou [4] applied a modified genetic algorithm and presented a new approach. The new approach searched for successive robot configurations in the entire free space so that the robot moved its end-effector from initial placement to a final desired one. It showed better performance than the pseudo-inverse method and simple genetic algorithm. Tabandeh et al. [40] proposed a modified genetic algorithm to find multiple solutions of the IK through an adaptive niching strategy. The algorithm used a few preset parameters and could be generalized to solve the IKP of a robot with an unknown number of DOF and an unknown present configuration.
Artificial Neural network is one of the most important techniques to solve the IKP with the help of different neural network architectures. Since the neural networks work with an acceptable error, the error at the end of IK learning should be minimized. Köker [41] used a simulated annealing algorithm to improve the result obtained from the neural-network-based IK solution of robotic manipulators. A simulated annealing algorithm was used to find the best-fitting ten digits for the decimal part of the best solution selected by three Elman neural networks. Using an artificial neural network, Srivastava and Kumar [42] presented an optimization technique to get the best results for the path of the end-effector with respect to the joint angles. Unlike traditional implementations of neural networks, Feng et al. [43] proposed a hybrid approach to train the neural network using the collected training data from joint subspace and an efficient learning algorithm, the electromagnetism-like method. The hybrid approach could significantly reduce computation time and improve precision.
The particle swarm optimization algorithm is one of the biologically-inspired optimization techniques that can be applied to solve optimization problems with a multimodal function. Huang et al. [44] applied a particle swarm optimization algorithm to solve the IKP of 7-DOF robotic manipulators. The algorithm was capable of finding an optimal configuration and solving the problem more efficiently. By merging genetic algorithms and particle swarm optimization algorithm, Starke et al. [45] developed a novel biologically-inspired approach to solving the IKP efficiently for arbitrary joint chains. The hybridization approach performed more robustly and adaptively than traditional or related methods. Liu et al. [46] proposed an IK calculation method based on an improved particle swarm optimization algorithm. The method is applied to general robots. Compared to the traditional particle swarm optimization algorithm, it solved the local convergence problem and improved stability and convergence accuracy, and operation speed.
Many of the methods mentioned above have been directed to minimize a criterion function. Taking additional VOLUME 10, 2022 constraints into account implies a time-consuming optimization process. Many numerical methods are based on inverse differential kinematics. They involve the computation pseudo-inverse of the Jacobian matrix, which has the problems of large computation amount and long computation time [47]. These methods can only get one approximate solution that depends on the initial guess and has a cumulative error owing to many iterations. The common disadvantage of numerical methods is the lack of repeatability. In other words, the joint space trajectory corresponding to a repeated task space trajectory can not repeat itself, which means that the behavior of a robot is unpredictable when periodic tasks are carried out [48]. Furthermore, numerical methods require the current joint positions of controllers and a well-planned motion path on which the poses of the end-effector in a specific interval are continuously computed. However, robot motion planning and task planning expect the joint configuration of the desired pose to be computed first, and then plan the motion path and trajectory of the end-effector. Therefore, it is difficult for general numerical methods based on inverse differential kinematics to realize them.
Analytical methods have the advantages of fast solving speed, high precision, and obtaining several reliable solutions. However, analytical methods are only applicable to 7R robots with unique geometric configurations, such as redundant humanoid arms with S-R-S structure, a robotic manipulator with all zero of joint offset or link length [49], or a robot with three consecutive parallel joints and spherical wrist joints [50]. They lack universality and portability. Even though the search methods based on intelligent algorithms can obtain multiple solutions or an optimal one, they are timeconsuming, and the obtained solutions are approximate.

B. WORK OF THIS PAPER
Although solutions to the IKP have been expressed in closed form for some 7R robots with a particular configuration, it is not easy to directly obtain analytical solutions to the IKP for 7R robots with general geometric parameters due to the highly nonlinear and complex coupling properties. Whereas numerical methods based on inverse differential kinematics can solve the IKP for 7R robots with arbitrary geometric parameters, it is time-consuming to compute one solution governed by the initial guess without numerical stability. In addition, search methods based on intelligent algorithms are time-consuming to obtain approximate solutions.
It is difficult to derive the kinematic and dynamic equations of robotic systems with high DOF. There are some useful works for a wheeled mobile robotic manipulator [51], planner multi-body systems [52], tree-type robotic systems [53], and mechanisms of single-loop or multi-loop linkages [54], flexible link manipulators [55], [56], etc.
To solve the IKP of 7R robots with general geometric parameters, we propose a new method combining symbolic preprocessing, Sylvester dialytic elimination [57] and eigendecomposition to find the IK solutions, called SDE method for short. Given the geometric parameters of a 7R robot and the end-effector's desired position and direction, the displacements of all seven joints should be computed. The IKP boils down to solving a system of multivariate equations, in terms of whose algebraic properties, according to the method proposed by Raghavan and Roth [58], 14 symbolic equations are generated. After symbolic computation to separate variables, we obtain two equal polynomial matrices and then break them down into two systems of polynomial matrices equations. Setting a joint as the discrete variable, inputting numerical values, and thereby applying the Gaussian elimination method, six equations are selected from the 14 reconstructed equations to form a system of linear equations whose coefficient matrix contains only one variable. A 12 × 12 matrix is structured, which contains only one variable by variable substitution of the half-angle tangent function and Sylvester dialytic elimination. Rearranging the 12 × 12 matrix to form a 24 × 24 numerical matrix, we compute its eigenvalues and eigenvectors. The position of each joint is computed successively through back-substitution, and finally, several IK solutions are obtained.
Although the SDE method is not an analytical method, it has the novelty of obtaining several IK solutions at once by almost real-time computation and further obtaining the unique solution under given constraint conditions. The proposed method for the IKP of 7R robots is a combination of algebraic and numerical techniques. It involves symbolic computation, linear algebra, numerical analysis, and matrix operation. It applies to the kinematics model of a 7R robot established according to Denavit-Hartenberg (D-H) parameters convention [59] of Waldron [60] and Paul [61], without considering joint limits. Compared with the works mentioned above, the SDE method has several contributions: • Several solutions at position level are obtained in one call.
• By contrast with numerical methods, solving time of the SDE method is shorter and numerical stability is better. It is beneficial for online motion planning and can be employed in a multi-objective optimization approach [62] for general 7R robots.
• The SDE method can solve the IKP of 7R robots with arbitrary geometric parameters and can be modified to apply to the IK of all 7-DOF serial robotic manipulators.

II. 7R ROBOT KINEMATICS MODEL
In order to solve the IKP for a 7R serial robot, the forward kinematics model should be established first. According to the D-H parameters convention of Waldron and Paul, as shown in Fig. 1, we model the geometric representation of a 7R manipulator. Seven links in a 7R manipulator are numbered from 0 to 7, so the same is the coordinate system attached to each link. The base link is 0, and the outermost link is 7. Seven joints in the robot mechanism are numbered from 1 to 7. The four D-H parameters are joint angle θ i , joint offset d i , link length a i , and link twist α i . For articulated robots, the last three parameters are definite values, and the first one is taken as the joint variable. The general transformation of coordinate system {i} relative to coordinate system {i − 1} is shown as (1), by multiplying four orderly independent elementary transformations based on the above four parameters, where Then, for a 7R manipulator, multiplying the transformation matrix of each link sequentially results in the transformation matrix of the coordinate system {7} attached to the outermost link relative to coordinate system {0} attached to the base link, where, Equation (2) is the formulation of forward kinematics for a 7R manipulator, which expresses the mapping relationship between joint variables and a pose of the end-effector. In other words, given values of all joint variables for a robot, we can use (2) to describe the pose of its end-effector with respect to its base. Conversely, the IKP corresponds to computing the seven joint angles satisfying (2) for the desired pose.

III. SDE METHOD
In this section, we present the process of SDE method to solve the IKP for 7R robots in detail. First, we rearrange (2), as As a result, each entry of the right-hand matrix of (4) is a function of θ 3 , θ 4 , θ 5 and θ 6 , and each entry of the lefthand matrix is a function of θ 1 , θ 2 and θ 7 . We represent a manipulator as a chain of characters that symbolizes its position of revolute joints and general geometry. Equation 4 helps to decrease the degree and the symbolic complexity of the resulting expressions.

A. SYMBOLIC PREPROCESSING FOR VARIABLES SEPARATION
Taking the D-H parameter link length a i , joint offset d i , sines and cosines of link twist α i and 1-3 rows of elements in 0 T 7 as symbolic constants, the entries of 3rd and 4th columns among left-hand matrix in (4) do not contain θ 7 . Thus, comparing the entries of the third and fourth columns among the corresponding right-hand matrix, we can obtain six equations with six joint variables, in which all constant terms are moved to the right side (see equations EQ01-EQ06 in Appendix A).
We rearrange them to structure four vectors, Considering p left , p right , q left and q right to be 3 × 1 vectors, according to the operation rules of dot product and cross product for vectors, the following four expressions are valid [63], From the above four expressions of vectors relationship, we can derive the other eight equations employing symbolic operation (see equations EQ07-EQ14 in Appendix A). Thus, we get 14 equations containing symbolic constants and variables, in which left-hand polynomials are a linear combination of sθ 1 sθ 2 , sθ 1 cθ 2 , cθ 1 sθ 2 , cθ 1 cθ 2 , sθ 1 , cθ 1 , sθ 2 and cθ 2 , and right-hand polynomials are a linear combination of sθ 4 sθ 5 , sθ 4 cθ 5 , cθ 4 sθ 5 , cθ 4 cθ 5 , sθ 4 , cθ 4 , sθ 5 , cθ 5 and 1. We express them as, where, 12 =(sθ 1 sθ 2 ,sθ 1 cθ 2 ,cθ 1 sθ 2 ,cθ 1 cθ 2 ,sθ 1 ,cθ 1 ,sθ 2 ,cθ 2 ) T , 45 = (sθ 4 sθ 5 ,sθ 4 cθ 5 ,cθ 4 sθ 5 ,cθ 4 cθ 5 ,sθ 4 ,cθ 4 ,sθ 5 ,cθ 5 , 1) T , Q is a 14 × 8 symbolic constants matrix and P is a 14 × 9 symbolic polynomials matrix. Q is rearranged from the left-hand polynomials of the 14 equations, and P is derived from the right-hand of the 14 equations. All entries of P are linear functions of sθ 3 , cθ 3 , sθ 6 and cθ 6 . Each entry of Q and P is a symbolic function structured from the D-H parameter, and the relationship in (5) helps to eliminate 4 of the 6 joint variables. Symbolic preprocessing is used to separate variables, and is applicable to 7R robots with arbitrary geometry as a universal method.All entries of P are linear functions of sθ 3 , cθ 3 , sθ 6 and cθ 6 . Each entry of Q and P is a symbolic function structured from the D-H parameter, and the relationship in (5) helps to eliminate 4 of the six joint variables. Symbolic preprocessing is used to separate variables and is applicable to 7R robots with arbitrary geometry as a universal method.

C. NUMERICAL SUBSTITUTION AND MATRIX COMPUTATION
As for the polynomial coefficient matrix P in (5), which contains two joint variables, θ 3 and θ 6 , it is necessary to eliminate one joint variable to form a matrix with only one joint variable. The solution is to choose one of them as a discrete variable and assign an initial value θ ini , such as zero or a nominal angle value, to it within the joint position limit first. If there are no IK solutions found, it is assigned a new value in accordance with a certain increment δθ . In this paper, we choose θ 3 and θ 6 as the discrete fixed joint parameter, respectively, and express it as ξ . Thus, given a specific configuration of a 7R robot, we substitute numerical values of D-H parameters and the sine and cosine values of discrete joint ξ into matrices Q and P, so that (5) is converted into, As a result, we obtain a numerical matrixQ and a numerical polynomial matrixP. All entries ofP are linear combinations of sζ , cζ and 1 (when ξ = θ 3 , ζ = θ 6 , or when ξ = θ 6 , ζ = θ 3 . sζ = sin ζ , cζ = cos ζ ). Based on the description of Section III-B,we can obtain numerical matricesQ 1 and Q 2 , the minors ofQ, and numerical polynomial matricesP 1 andP 2 , the minors ofP. Since numerical matricesQ 1 and Q 2 are derived from symbolic operation, and each entry of them is computed at one time by numerical substitution, they are definite given geometric parameters of a robot and a pose of its end-effector. Also, each entry ofP 1 andP 2 is a linear combination of sζ , cζ and 1, whose numerical coefficients are definite. Thus, matricesQ 1 ,Q 2 ,P 1 andP 2 can be operated according to algorithms of numerical matrix.

1) GAUSSIAN ELIMINATION
We rearrange (8), and obtain two equations. They are expressed asQ As a system of linear equations with respect to variables sθ 1 and cθ 1 , (9) has two variables but six equations, and by Gaussian elimination, four equations can be eliminated. Similarly, two equations can be eliminated in (10). We perform Gaussian elimination with partial pivoting method on (9) and (10) respectively. After employing elementary row operations onQ 1 ,Q 2 ,P 1 andP 2 , we obtain two equations, Since row rank ofQ 1 may be 0, 1, or 2, at least there are 4 rows whose elements are zeros inQ 1 , and the all elements of last four rows are zeros. Since row rank ofQ 2 is at most 6, at least there are 2 rows whose elements are all zeros, and all elements in the last two rows are zeros. Each entry of numerical coefficient polynomial matricesP 1 andP 2 is linear combination of variables sζ , cζ and 1, like the form ηcζ + µsζ + δ.

2) SYLVESTER DIALYTIC ELIMINATION
After performing Gaussian elimination, left-hand polynomials of the last four equations in (11) are equal to 0, whereas the left-hand polynomials of the last two equations of in (12) are equal to 0. We select the last four rows of elements fromP 1 and the last two rows fromP 2 to structure a 6 × 9 matrix . Each entry of is a linear combination of sζ , cζ and 1. As a result, we obtain equation (16), Replacing the trigonometric function in (16)  where y = tan ζ 2 , x 4 = tan θ 4 2 , x 5 = tan θ 5 2 , and multiplying each equation by 1 + y 2 , 1 + x 2 4 and 1 + x 2 5 to clear out the denominators, (16) is converted into the form shown as, where (17), is 6 × 9 matrix, each entry of which is a quadratic polynomial of y. In order to perform eigendecomposition later, it is necessary to construct a new square matrix from using Sylvester dialytic elimination. Multiplying (17) by x 4 , we obtain a square system of the form, where . is a 12×12 matrix, each entry of which is a quadratic polynomial of y. The condition for (18) to exist non-zero solutions is that the determinant of is zero. Therefore, a degree 24 characteristic polynomial in y is equal to 0. Computing the determinant of can introduce significant numerical errors, and the computation of real roots of the polynomial can be ill-conditioned [64], [65]. Furthermore, even though the roots of y can be obtained by computing the determinant of , the next step we back substitute the numerical value of y into in (17) and then we have to compute determinant for x 4 . As a result, the numerical errors accumulate in the intermediate steps of the computing determinant and the roots of polynomials. Therefore, how to compute y from has a significant impact on the computational efficiency and numerical stability. It is difficult to solve the roots of a high degree univariate polynomial, but in this paper, we reduce the problem of computing roots to the problem of eigenvalues computation. Since each entry of is a quadratic polynomial of y, we express it as where, M , N and L are all 12 × 12 numerical matrices.

3) EIGENDECOMPOSITION AND BACK-SUBSTITUTION
The condition for inverse matrix M −1 of matrix M to exist is that its condition number is not infinite. That is to say, M is a well-conditioned matrix. When this condition is satisfied, the unary polynomial matrix Iy 2 + M −1 Ny + M −1 L in y exists. If M is ill-conditioned, we can set a new value to ξ according to ξ = ξ ini + δξ n cnt in Algorithm 1, where n cnt is an integer between the maximum negative and positive increments computed from the lower limit and the upper limit of the free joint. Considering the case that M is wellconditioned, we structure a 24 × 24 matrix from (19), According to the research of Gohberg et al. [66], the 24 eigenvalues of matrix are exactly the roots of Algorithm 1 To Generate a New Guess for the Free Joint Require: ξ ini , the initial guess for the free joint ξ upper , the upper limit of the free joint ξ lower , the lower limit of the free joint δξ , the search discretization increment Ensure: ξ new , a new guess; 1: Initializing: ξ new ← 0, n cnt ← 0 2: n pos ← FLOOR((ξ upper − ξ ini )/δξ ), getting the step number from the initial guess to the upper limit 3: n neg ← CEIL((ξ lower −ξ ini )/δξ ), getting the step number from the lower limit to the initial guess 4: n cnt ← GET_COUNT(n pos , n neg , n cnt ), to generate a new counter between n pos and n neg 5: ξ new ← ξ ini + δξ n cnt , to generate a new guess for the free joint 6: return ξ new determinant( ) = 0. That is to say, the values of variable y correspond exactly to the eigenvalues of matrix . In addition, eigenvectors of matrix compose the following structure [67], where v happens to correspond to X 45 . As a result, we can compute all the roots of y, x 4 and x 5 in (18) by means of the eigendecomposition of . We pick out the real roots to compute ζ , θ 4 and θ 5 . Back substituting sines and cosines of ζ , θ 4 and θ 5 into (8), θ 1 , θ 2 can be solved. At last, θ 7 can be solved using (2). Since V contains several terms consisting of x 4 and x 5 , it is necessary to select a reasonable pair from them to compute x 4 and x 5 . Each term of v has the same bound on the maximum error according to the study of Wilkinson [68]. The term with the maximum magnitude generally has a smaller relative error, and this property is used to compute x 4 and x 5 accurately. Dividing V into two parts, v and y v, let 10 , v 11 , v 12 ) T , which corresponds to X 45 . thus, we can compute x 4 and x 5 from v 1 , whose relative error is low. In vector v 1 , the last element is usually not equal to 1, so we can not directly use the 9th and 11th elements to compute x 4 and x 5 , respectively. To reduce the error, x 4 or x 5 is computed by selecting elements with larger magnitude from v 1 . That is to say, if v 1 has the maximum magnitude, then x 4 is equal to the ratio of v 1 and v 4 , and so on. Two elements are selected, the ratio of which is x 4 . A similar operation is performed to compute x 5 . After obtaining sζ , cζ , sθ 4 , cθ 4 , sθ 5 and cθ 5 , back substituting them into (8) and (2), the values of other variables are solved.

D. CONVERGENCE ANALYSIS
In the previous sections, we have presented a practical method for the IK of general 7R serial manipulators. Currently, the method is restricted to manipulators consisting of seven revolute joints, and it does not assume the geometry of the manipulator. Based on the algebraic formulation in Section III-A, we reduce the problem to computing zeros of matrix polynomials and matrix pencils. The matrix polynomial is regular if its symbolic determinant is non-zero. Otherwise, it is a singular matrix polynomial. However, it is difficult to check whether is singular. As an alternative approach, we compute the condition number of numerical matrix M . If the matrix is singular, its condition number is infinity. We set the rational maximum condition number as n con and take it as a prerequisite for the eigendecomposition of . If the condition number of M is less than n con , the next computation will carry on. Otherwise, according to Algorithm 1, a new M is generated corresponding to the new value of the free joint. Because the free joint has the upper and lower limit, the amount of total iteration between them is determined, and it is the sum of n pos and n neg . The sum relates to the search discretization increment, δξ . For a smaller δξ , the SDE method needs more iterations, but it obtains better performance in success. In addition, if the solution result does not meet the final error tolerance, the iterative process will also continue. In any case, the SDE method ends after a finite number of iterations, regardless of whether the solving operation is successful. But in most cases, this method succeeds in computing the IK solution, as shown in Table 5 and Table 6.

IV. SDE METHOD IMPLEMENTATION
In this section, combining the Robotics Toolbox for MATLAB developed by Corke [69], we realize the SDE method proposed in Section III by programming in the mathematical computing software of MathWorks, MATLAB R2019B.
We use Baxter, a dual-arm robot in our laboratory, as an example to illustrate the solution to the IKP with the SDE method. Analyzing the URDF files in the control system of Baxter, we can obtain its D-H parameters. When the upper shoulder link is taken as the base link, its left arm and right arm have the same D-H parameters, and their installation positions have the same coordinate transformation, transl(0.055695, 0, 0.011038)rotz(0)roty(0)rotx(0). In this case, the SI unit of angle is rad and the SI unit of length is m. We take the left arm as the research object and list its D-H parameters in Table 1. The corresponding D-H parameter kinematics model is shown in Fig. 2. We

A. SPECIFIC EXAMPLE
Let the free joint variable be θ 6 and θ 6 = 1.83, ζ = θ 3 . Substituting D-H parameters and the elements of 0 T 7 into symbolic matrices Q and P, after operation according to Section III, we can obtain eight real eigenvalues about θ 3 . The results is shown in Table 2. And the 8 solutions of joint variables for the specific pose of the left end-effector and the corresponding pose errors are computed and given in Table 3.

B. GENERAL EXAMPLE
Let the free joint variable be θ 3 and θ 6 , respectively. The initial guess θ ini is zero, and the discretization value δθ for searching the valid free joint displacement is 5 degrees. In Robotics Toolbox for MATLAB, tr2delta(T 0 , T 1 ) is the function to compute the differential motion corresponding to infinitesimal motion from pose T 0 to T 1 which are homogeneous 4 × 4 transformations. We can use it to compute the error of the pose corresponding to the computed joint values relative to the desired pose T eef . The error is expressed as a vector (δx, δy, δz, δR x , δR y , δR z ) × 10 −12 , which represents  infinitesimal translation and rotation. The final error tolerance for the pose of a valid IK solution is expressed as . It is computed on the norm of the error between the current and desired pose, and = 10 −10 . Failure may happen when the SDE solver runs at the initial guess for the free joint. In this case, the value of the free joint can be set according to Algorithm 1. In this way, the success rate will improve markedly. We substitute D-H parameters and the pose data into the computation program. When ζ = θ 6 , 8 valid eigenvalues are computed after one iteration, which are 4.4597,  Table 4. We can choose the No.5 solution as the best one of the 12 solutions, which has the max manipulability index. In addition, for several special configurations, such as zero joint angle configuration, vertical configuration, stretched configuration, and nominal configuration, the SDE method can find the IK solution, too.

C. MORE EXAMPLES
We test the SDE method on seven types of 7R robots which are selected from the robots models library of RoboDK V5.2.2, a simulation software. The seven robots are Daihen OTC FD-V6LS, Franka Emika Panda, Kinova Gen3, KUKA LBR iiWA 14 R820, Motoman SIA20D, Siasun SCR3 and Yumi IRB 14050. The standard D-H parameters of these seven robots are imported into the MATLAB program of the SDE method, and the pose of the end-effector is computed according to the configuration data of the parameter panel for each robot in Fig. 6 to Fig. 12 (see Appendix B). For each robot, we solve the IKP corresponding to the pose of its endeffector, and we list the IK solutions in Table 12 to Table 18 (see Appendix B), respectively.
Generally, the wrist structures of 7-DOF manipulators are designed in two different types: the Euler wrist and the offset wrist. Offset wrists are used with the robot manipulators that are needed to get long horizontal reach while maintaining the appropriate angle in their workspace, while the Euler wrists provide regularly shaped workspace [70]. Solving the IKP is very difficult for the robot manipulators that have an offset wrist. We test the SDE method on the YASKAWA VA1400II robot [22] that has the offset wrist (see Appendix C).

V. EXPERIMENTS AND COMPARATIVE ANALYSIS
In this section, we carry out simulation experiments and one practical experiment. The simulation experiment is a comparative study on the performance of the SDE method and the IK tools in the Robotics Toolbox for MATLAB (version 10.4).

A. PERFORMANCE
Robotics Toolbox for MATLAB is a MATLAB-based robot modeling and simulation toolbox developed by Peter Corke from Australia. It contains three built-in functions, ikine, ikunc, and ikcon, to solve general IKP. Function ikine applies the Gauss-Newton iterative algorithm [71] optimized by the Levenberg-Marquardt method [72] to compute IK solutions. Function ikcon searches the minimum of nonlinear multivariate objective function considering joint limits, whereas function ikunc does the same but does not consider joint limits. Function ikunc applies the Quasi-Newton method [73], whereas function ikcon applies the SQP method [74]. We take the kinematics model of the left arm in Section IV as the simulation object. At first, 1000 sets of joint positions are generated by random function within the joint limits. Each configuration corresponds to a pose of the end-effector, so we have 1000 configurations and 1000 corresponding poses. The final error tolerance are = 10 −2 , = 10 −3 , = 10 −4 , = 10 −5 , = 10 −6 , = 10 −7 , respectively. We compute the solutions corresponding to each pose using the SDE method and the three IK functions, respectively. The configuration of zero joint displacements is considered the VOLUME 10, 2022   initial guess for these four IK solution methods. Computing results are listed in Table 5. From Table 5 we can see that the numbers of successes for ikunc and ikcon decrease when the convergence condition gets stricter. For = 10 −10 ikunc and ikcon will fail, so Table 6 only compares the SDE method to ikine for 10000 random poses. When ξ = θ 3 , the average number of iterations is 1.0004 before the SDE method finds the IK solution, whereas the average number of iterations is 13.8042 when ξ = θ 6 .
To study the data in Table 5, we can see that functions ikine, ikunc, and ikcon cost more time to search IK solutions. Compared to the numerical methods in Robotics Toolbox for MATLAB, the running time of the SDE method is very short, and meantime, it ensures a good success rate for the IKP. On a personal computer, configuring with CPU of Intel Core TM i7-4712MQ@2.3GHz and 16G memory, the average running time is reduced to 20ms, whereas the other three inverse methods take ten to twenty times as long as SDE method. Fig. 3 shows the experimental platform, and it is composed of a Baxter robot, a Leica laser tracker, a PC workstation with the Robot Operating System (ROS) to control the robot, and a PC station to collect the data generated by the laser tracker.

B. EXPERIMENTAL VERIFICATION
At first, in the MoveIt tool, we choose 15 groups of different joint positions according to Fig. 4. Each group of joint positions corresponds to a pose of the end-effector, and the SDE method is used to compute the IK solutions. We select a reasonable one conforming to the joint limit conditions from all IK solutions for each pose. Then, the workstation configured with ROS and Baxter SDK is connected to Baxter through the network protocol, and the reflector ball of the laser tracker is mounted at the end of the left arm. The terminal executes two Python scripts containing the joint positions generated by the MoveIt tool and computed IK solutions by the SDE method. The robot is controlled online to reach 15 configuration states successively. The laser tracker captures the Cartesian space positions of the end-effector, and two sets of point cloud data are obtained. The former set of data is taken as the reference data, and the deviation of the latter set of data is measured.
The data in Table 7 are joint positions generated in the MoveIt tool, and the spatial position of the end-effector corresponding to each configuration measured by the laser tracker is listed in Table 9. The data in Table 8 are selected from IK solutions solved by the SDE method, and the measured spatial position of the end-effector for each configuration is listed in Table 10. We regard the data in Table 9 as reference data, and compute deviation of the data in Table 10. The deviation results are listed in Table 11. Fig. 5 shows the positions of two sets of data in the same coordinate system. The red points are   reference data, and the blue points are the data generated by the SDE method. In Fig. 5, the measured data are very close to the observational data.

C. ANALYSIS AND DISCUSSION
In Section V-A, the simulation shows that the SDE method can solve the IKP in a shorter time than general numerical methods in Robotics Toolbox for MATLAB, and it can successfully solve the IKP of any poses in a workspace. Moreover, it has good numerical stability, better accuracy, and computational efficiency. Due to Algorithm 1, the SDE has a nearly 100% success rate, and its numerical stability is improved. For the same initial guess, the SDE method spends several iterations before it finds the solutions, and VOLUME 10, 2022 its computational complexity is considered to be O(n). Even though the error tolerance is getting smaller, the average time to compute the joint angles does not change significantly. By contrast, ikine has an almost 40% probability of failing because it reaches the max iteration number, and ikunc and ikcon have the worse performance in searching for a higher precision solution. The proposed method is more efficient than the three numerical methods due to the following reasons: 1) Convergence The SDE method has a certain number of iterations. Although the three numerical methods set a maximum number of iterations, they can not guarantee convergence before meeting the end condition. 2) Time Cost In Table 5, the SDE method spends the shortest time among the four IK methods. 3) Tollerance The SDE method can compute the IK solution at higher precision and does not affect the solving success rate. The numerical methods fail probably to obtain higher precision. 4) Success Rate The SDE method has a better performance in solving the success rate than the three numerical methods.
To compare the configuration data provided by RoboDK and the corresponding IK solutions, we can see that the SDE method is suitable for 7R robots with different geometry. Given the desired pose of the end-effector, the SDE method can obtain several IK solutions, some of which accord to the joint limits and some do not. Imposing limit constraints can reject the IK solutions that do not satisfy the robot structure. In addition, the value of joint 6 in Table 12 to 18 has an offset plus the joint position provided by RoboDK, respectively. If the offset is zero, there must be one among computed IK solutions almost consistent with the data in RoboDK. Thus, it proves that, under the condition of the existence of IK solutions for a 7R robot, the SDE method can solve the IKP successfully.
Considering that the position accuracy of the Baxter robot is about ±0.5 mm and the manipulator is not stationary when   it has been in the desired state because of the series elastic actuators, the result in Table 10 is acceptable. The average deviation is 7.942 mm, and we can see that the IK solution obtained by the SDE method can make the robot reach a desired spatial position. The main interest of our work is on 7-DOF serial robotic manipulators with only revolute joints. We hope our work can contribute to the IK solution for general 7R robots.

VI. CONCLUSION
1) To solve the IKP of a 7R robot with arbitrary geometry, analytical methods are not applicable. In contrast, the numerical methods based on inverse differential kinematics can only find an approximate solution given an initial guess, and they are numerically unstable. This paper proposes a general method called SDE, applicable to 7R robots with arbitrary geometry. The SDE method can compute several IK solutions by combining symbolic preprocessing, Sylvester dialytic elimination, and eigendecomposition. It employs symbolic preprocessing to separate variables, Sylvester dialytic elimination to construct a square matrix of coefficients, and eigendecomposition to compute roots of the specific variables. Finally, the IKP is solved by matrix operation and back substitution. 2) Comparing to the analytical methods only capable of solving the IKP of 7R robots with a particular configuration, the SDE method has general applicability to 7R robots with arbitrary geometric parameters. 3) Comparing to the solution methods based on inverse differential kinematics, the SDE method is almost realtime to obtain solutions and is more effective and numerically stable. Given the pose of the end-effector and an initial discrete value of the free joint, the SDE method can compute several IK solutions that are suitable to be taken as initial states for motion planning. 4) SDE method is proposed for general 7R robots without considering joint limits. Although it applies to the kinematics model of 7R robots established according to the D-H parameter convention of Waldron and Paul, after parameters converting, it is also applicable to the kinematics model of 7R robots established according to the D-H parameter convention of Craig [75]. This method has the potential to provide a feasible way to solve the IKP of hyper redundant robots and redundant robots with prismatic joints. 5) According to the basic D-H parameter description, we know that each joint has four parameters, namely θ, d, a, and α. So for a 7R robot, there are 28 design variables in total. Seven of them vary, and the rest fix and determine the robot's geometry. However, due to manufacturing errors, the number of geometry parameters of each joint will be expanded to six. The two added are the offset in theŷ axis direction and the rotation around theŷ axis. In fact, after the calibration process, the values of the two extra parameters can be determined. In this way, although 42 design variables are generated, in the symbol processing stage, by considering six θs as variations and the rest as constants, 14 equations similar to those in Appendix VI can still be derived. Of course, these equations are a bit more complicated, but the derivation process is similar. Then, in the following procedure of symbolic matrix analysis, numerical substitution, and matrix calculation, at least six equations with the right side of zeros can be obtained. Thus, an algebraic matrix with only one variable y is formed. Finally, is transformed to a 12 × 12 matrix through the Sylvester dialytic elimination, and its dimension corresponds to X 45 , neither increasing nor decreasing. 6) The algebraic matrix contains only one variable y, the equation whose determinant is zero is of a high order, and it is not easy to find all the roots. Even if the value of a joint angle variable can be found by constructing a univariate higher-order equation, substitute it back to (17), and then construct a higher-order equation about θ 4 or θ 5 , and calculate the value of each joint angle in turn. The whole calculation progress will generate a more significant calculation error. By calculating the condition number of the matrix M in (19), M is judged whether ill-conditioned or not. If wellconditioned, the matrix is constructed, and three joint variables are calculated at one time through eigendecomposition. Otherwise, a new matrix M is generated by changing the angle of the free joint without complicating the calculation process. Therefore, although there are unreasonable coefficients, the calculation of the three joint angle variables with significant errors can be avoided by the condition number calculation. In addition, the SDE method can filter out unnecessary inverse solutions by setting an error threshold after obtaining multiple inverse solutions. 7) For the Baxter robot, the SDE method can find the IK solution after fewer iterations when ξ = θ 3 . Nevertheless, a failure situation happens when the free joint variable is θ 3 for robots Daihen OTC FD-V6LS, Kinova Gen3, KUKA LBR iiWA 14 R820, Motoman SIA20D, and Siasun SCR3. In many cases, when the free joint variable is θ 6 , the SDE method can find the IK solution faster but with more iterations. Further research is needed before we can explain why the above situation occurs.
By the forward kinematics, the pose of the end-effector for the given configuration is Let the free joint is the third and the sixth joint, respectively. The IK solutions computed by the SDE method are listed in Table 20.