Extended singular robust inverse solution of redundant serial manipulators

We propose the use of a damping matrix β to replace the unified damping coefficient λ in the singular robust inverse method, allowing for adjustment of the manipulator’s singularity for different joints separately. In addition, the correction coefficient matrix α is employed to eliminate the Jacobian matrix approximation error, which usually exists in the iterative control solving process; the extended singular robust inverse method is thereby proposed. This article contains analyses of the singularities caused by three methods and how they result in larger values of the manipulator terminal error and joint angular velocity; additionally, the semi-empirical selection principles of β and α are given. The stability of the proposed method is investigated with the Lyapunov stability criteria, and its effectiveness is verified for spatial linear and curve trajectories with 7-axis and 10-axis redundant serial manipulators. The simulation results show that the proposed method exhibits improved optimization performance and stability, and the terminal error is stable and within the allowable range. The least singular value is larger and the joint angular velocity is the lowest, which is the expected outcome. Furthermore, the simulation results for the 10-axis manipulator also indicate the generality of the proposed method.


Introduction
Kinematics control is one of the most general issues in robot manipulator technology, in which the focus is on the inverse kinematics solution. [1][2][3][4][5][6][7][8] However, it is often difficult to find a suitable solution for redundant manipulators due to high nonlinearity. There have been many studies on classical methods used to find the pseudoinverse solution of various target functions. 2,3,6,[9][10][11][12][13] The singular value decomposition generalized inverse method (SGIM) is a typical and widely used method for solving inverse motion of redundant manipulators, and it provides a defined Jacobian generalized inverse by singular value decomposition. 8,10,14,15 The optimization object of SGIM is the least squares solution with the minimum norm. 10,15 However, when the minimum singular value approaches 0, the norm of angular velocity tends to infinity, and the resulting solution has no practical significance for the control operation.
Therefore, another new scheme, the singular robust inverse method (SRIM), has been proposed and widely used in inverse kinematics of redundant manipulators. [11][12][13][16][17][18][19][20][21][22] Unlike SGIM, SRIM introduces the damping coefficient l to avoid problems caused by small singular values and reconstructs a new optimization target function. However, the optimization effect is largely decided by the selection of the damping factor l and, in recent years, scholars have proposed various definitions of l. Wampler suggested that l was a constant. 11 Nakamura and Hanafusa, 17 Kelmar and Khosla,18 Chan and Lawrence, 19 and others 13,20 have proposed a definition of the piecewise function of l. In addition, intelligent algorithms have also been used in recent years to find better damping factor values. 21 Compared with SGIM, SRIM has achieved better results in solving the inverse motion of trajectory tracking, but there also exist some shortcomings. They are as follows: 1. A uniform value of damping factor l is usually employed in SRIM solutions without considering the optimization requirements of each singular value. As a result, the over-optimization of large singular values results in redundant errors and singularities. 2. In order to effect the change of joint angle Du from t k to t k + 1 during the iteration process, the Jacobian matrix is usually considered unchanged and the Jacobian matrix at t k is considered an approximation substitution value.
Once the manipulator passes a singular point between t k and t k + 1 , the singular values of the Jacobian matrix will change significantly and the approximation will produce a larger terminal error. The error can be eliminated by employing a smaller sampling interval, but more computing resource will then be consumed. Obviously, neither SRIM nor SGIM considers the effect of the approximation.
To eliminate the defects of SRIM, we propose to introduce the damping matrix b instead of the unified damping factor l. A pre-factor matrix a is also employed to eliminate the error caused by approximation of the Jacobian matrix and the extended singular robust inverse method (ESRIM) results. Subsequently, a comparative analysis and simulation examples show that using the proposed method, the stability of the ESRIM is significantly improved and results in the desired end error, smaller joint angular velocity, and larger least singular value.
This article introduces the modified coefficient matrices to overcome the shortcomings of the SRIM. An introduction to the SGIM and SRIM models is shown in section ''Preliminaries'' and the new methodology ESRIM is illustrated in section ''New singularity-avoidance scheme.'' In section ''Simulation analysis verification,'' the modified method is applied to a 7-axis manipulator and a 10-axis manipulator, respectively, for validation and comparison. Section ''Conclusion'' presents the conclusions drawn from this work.

Preliminaries
The kinematics study of manipulators provides the mapping relationship between joint space and manipulation space 4 and is highly non-linear for redundant manipulators. Therefore, it is necessary to linearize the non-linear motion to effect satisfactory motion control.
The common method of linearization is to establish the equation with the Jacobian matrix. 23 If divided by time, the linear relationship of the corresponding velocity can be obtained as follows where J(u t ) is the m 3 n partial derivative Jacobian matrix, m is the row number of Jacobian matrix and n is the column number. Furthermore, it can be defined as follows where i = 1, 2,..., m and j = 1, 2,..., n.
According to the theory of linear algebra, 24 when _ x t is within the column space C(J(u t )), equation (1) is compatible and has exact solution _ u t . Otherwise, only the least squares solution, which minimizes the residual of the equation is obtained. The residuals result from the unreasonable task planning and the singularity of manipulators. The solution of joint-space velocity is composed of a homogeneous solution and a special solution, and the form is as follows where _ u 0 is the homogeneous solution belonging to the null space N(J(u t )). _ u Ã is a special solution which belongs to the row vector space C(J(u t ) T ), and neither of them is unique. _ u Ã directly determines the position accuracy of the end of the manipulator, and the optimal selection of its value will not produce the final selfmotion; this is often regarded as the priority option for the optimization. The commonly used methods SGIM and SRIM are optimizations for this term and ESRIM also effects the optimization of _ u Ã . The differential form of formula (1) can be written as Then there is where J(u) * (abbreviated J * ) is the inverse of the corresponding Jacobian matrix J(u) (abbreviated J).

SGIM
J is the rectangular matrix, and its inversion is different from that of the general square matrix. In 1955, the British mathematician Penrose 25 gave the definition of generalized matrix inverse which is extended to the rectangular matrix on the basis of Moore's research, 26 and they presented the uniqueness of the M-P generalized inverse (or '' + '' inverse). It can be obtained by singular decomposition 10 where U is an orthogonal matrix of m 3 m order, V is an n 3 m order orthogonal matrix, S + = diag (s À1 1 , s À1 2 , Á Á Á , s À1 r ) (which is a generalized diagonal matrix), and s 1 , s 2 ,..., s r are the singular values of the matrix J (which are decreasing in turn). Finally, r is the rank of J and r ł m.
For the full-row rank matrix, the optimization of the equation (1) with the SGIM is the least squares solution of minimum norm and can be expressed as The least solution norm and the least terminal error norm, corresponding to equation (7), can be expressed as where ½x 1 , . . . , x m T = U T _ x and other parameters are defined as above.

SRIM
The classical SRIM, also known as the damped least squares method, 12 has improved the membership of the minimum norm and least squares in the optimization target in equation (7) to be a parallel relationship. Both can be balanced by employed damping factor l. Therefore, a new optimization target function is constructed as follows Then, the optimization of the solution of equation (1) can be expressed as The solution norm jj _ u l jj 2 and the terminal error norm jjd l jj 2 with damping factor l are expressed respectively as with the parameters defined as above. From the formula (12), it can be seen that the introduction of the damping factor l can keep the norm jj _ u l jj 2 within a reasonable range when the least singular value s m tends to 0. For a full-row rank matrix, the terminal error norm of the M-P inverse solution jjdjj 2 = 0 , while jjd l jj 2 is no longer 0 in formula (13). Therefore, the introduction of damping factor l eliminates the singularity of the manipulator but increases the terminal error.

Iterative solving process
The inverse motion solutions based on the Jacobian matrix essentially constitute an iterative process, in which the angular variation can only be solved in a stepwise fashion. 27 The main procedural steps are as follows 1. The pose difference value (Dx) is calculated based on the initial pose (x t ) of the current moment and the desired pose (x t + 1 ) of the next moment; Dx = x t + 1x t . 2. According to the inverse solution formula (5), the corresponding joint angle change Du. is calculated: Du = J * Dx. 3. The joint angle u Ã t + 1 is calculated at the next moment: u Ã t + 1 = u t + Du. 4. The actual position (x Ã t + 1 ) of the end at the next moment is calculated by the positive kinematics formula: The end error jjdjj 2 is calculated by jjdjj 2 = jjx t + 1 À x Ã t + 1 jj 2 , and is used to determine whether the accuracy requirement is met.
If it is, the calculated value is used as the initial value to solve the next point; if not, the parameters are adjusted to calculate a new Jacobian matrix.

Thus, in
Step (2), the following approximation formulas should be used New singularity-avoidance scheme In this work, the inverse solution of equation (1) is obtained with an improved numerical method based on SRIM, named ESRIM. The optimization model in ESRIM can be expressed as follows In equation (15), two modified coefficient matrices a and b are introduced, both of which are diagonal matrices The damping matrix b replaces the unified damping factor l of SRIM. A pre-factor matrix a is also added into the target function, in order to eliminate the error caused by the change of singular values when the manipulator passes singular points. The optimization target function is constructed as In equation (17) Since J T a T aJ + b T b is positive definite in the formula (17), G( _ u) ab is a positive-definite quadratic form of _ u and has a minimum value only if _ u = _ u ab . Then the solution of formula (1) for optimum configuration (15) is as follows The extended singular robust inverse of the Jacobian matrix corresponding to equation (19) is Furthermore, J* ab can be expressed as a simple calculation, 16 as follows Then, for any Jacobian matrix J, the minimum norm jj _ u ab jj 2 and the end error jj _ d ab jj 2 of the extended singularly robust inverse solution are The approximate modified Jacobian matrix J a = aJ(u) is called the truth Jacobian matrix. Equation (4) can be expressed as The pre-factor matrix a is introduced to correct the approximation of the Jacobian matrix at t k + 1 . For the ith axis, it is reasonable to approximate the correction coefficient a i from the change of s i during the time interval, which is the ith singular value of the Jacobian matrix. In this work, the pre-factors are calculated as follows where a 0 is an adjustable critical pre-factor. In order to prevent the singularity caused by the small singular values of the matrix, the truth Jacobian matrix J a is improved by adding the damping matrix b, and the modified Jacobian matrix is denoted as J ab . Introducing b gives rise to an inevitable computational error d(J a ), which is expressed as The relative error d Du ð Þ=Du of the solution of equation (15) satisfies the following formula When the value jjJ Ã a jj jjd(J a )jj \1, the selection of b should minimize min b jjJ Ã a jjjjJ ab À J a jj, and this goal formula can be transformed as Obviously, when b i = 0 (i = 1, 2,..., r), formula (28) reaches its minimum value of 0. A non-zero value of b i helps avoid the singularity caused by small singular values, but it also leads to a larger terminal error. In this work, b i is defined as a piecewise function of s i as follows where i = 1, 2,..., m, and s 0 is the smallest singular value allowed. Equation (29) avoids the terminal error caused by over-optimization when singular values are large enough.
For the full-row rank Jacobian matrix J, the least norm jju ab jj 2 and the terminal error norm jjd ab jj 2 of the extended singular robust inverse solution can be reduced to where ½x a1 , . . . , x am T = U T a À1 _ x, s j satisfies s j ł s 0 and s j + 1 . s 0 . The inverse solution obtained by ESRIM takes the current singular value of the Jacobian matrix into account, and the overlarge value of jju ab jj 2 near the singular point is numerically avoided. In addition, the solution dose not exhibit over-optimization error.
In the process of trajectory tracking, with the solution error, _ u can finally be expressed as The stability of the proposed method is also analyzed. According to the stability criterion, [28][29][30] the Lyapunov function 29 is defined as Then, there are where K i and K p are both positive constant coefficients and are determined by the actual speed. V is always positive definite. Whether positive or not, _ V determines the stability of the proposed method. Obviously, we can choose the appropriate values of a 0 and s 0 , which make up of a and b, to satisfy the stability requirement of _ V as negative definite.

Simulation analysis verification
Taking the 7-axis manipulator ''YL101'' (designed in our laboratory) as the object, the effects of ESRIM are studied and compared with those of SGIM and SRIM; the latter are well known as traditional methods providing good inverse solutions of trajectory tracking. ''DH'' method was first proposed by Denavit and Hartenberg. 31 It is the most common and mature method to describe the structure of chain manipulator. 32 The structure of the 7-axis manipulator ''YL101'' is shown in Figure 1 and the DH structural parameters are listed in Table 1.

Example 1: linear trajectory
In the first simulation example, the tracking trajectory of the manipulator is a spatial straight line L S , ranging As a result of the accuracy requirement, the step size is set as 0.05 mm. The SGIM has no specific parameter settings. In SRIM, the damping coefficient is set as l = l 0 jjdjj 2 according to the recommendation of Chan and Lawrence, 19 and l 0 = 2:22 to achieve the best results. In the proposed ESRIM calculation, the critical pre-factor and singular value are a 0 = 2.30 and s 0 = 2.75, respectively. The simulation results are shown in Figures 2-7.
The terminal errors resulting from the three different methods are compared in Figure 2. Before the manipulator moves close to the first singular point, all three methods obtain good trajectories with acceptable small terminal errors. However, only ESRIM smoothly finds the solution over the entire process. Figure 2(a) shows the vicinity of the first singular point S P1 . When the manipulator approaches S P1 , the terminal error of the SGIM solution oscillates rapidly and then diverges until it exceeds the allowable error value; this leads to the failure of the whole optimization process. Adding the damping factor l allows SRIM to find the solution around S P1 . However, large oscillation occurs in the subsequent SRIM solution, and this violates the accuracy requirement around point S P2 (Point b in Figure  2). Compared with the former two methods, ESRIM obtains the solutions of S P1 and S P2 simultaneously. The terminal error does not exhibit oscillation with long duration or large amplitude and shows better optimization performance as the optimization process continues. Figure 3 shows the Lyapunov objective function derivative values _ V for the different optimizing methods. _ V should be negative, which is necessary for the stability of the optimization process. Before the first singular point, the _ V values of the three methods are all negative definite. However, the values of SGIM and SRIM oscillate and turn positive and eventually lead to the failure of the subsequent iteration process. Only the proposed ESRIM maintains the negative _ V and satisfies the stability requirement.
In addition to the optimization of terminal error, the norm of joint angular velocity should be small enough to be practical. As shown in Figure 4, the joint angular velocity for SGIM increases rapidly at the singular point S P1 , and this is consistent with the variation in its error.   The same relationship can be observed for joint angular velocity and terminal error with SRIM. Compared with SGIM and SRIM, the joint angular velocity of the proposed ESRIM is much more stable and remains a low level even at both of the singular points.
Equation (8) indicates that smaller singular values will lead to larger joint angular velocity and reduced flexibility of the manipulator. Therefore, for a certain terminal position, the least singular value s m of the Jacobian matrix should be as large as possible. Figure  5 shows the variation of the least singular values with the three optimization methods, and ESRIM produces the largest s m during the entire iteration process. These results show that ESRIM improves the flexibility of the manipulator. Figure 6 shows the track of the manipulator when different optimization methods are used to produce a spatial linear trajectory. The three tracks are generally consistent, but only the ESRIM successfully passes both singular points and completes the whole trajectory.
In addition, Figure 7 shows the trends in the angles of each joint during the optimization with ESRIM. All the joint angles change smoothly without sharp variations, and they remain within the range limits.

Example 2: curve trajectory
In the second example, the tracking trajectory is set as a spatial arc L C1 , ranging from the initial point    Because of the accuracy requirements, the arc-length step size is set as 0.05 mm. SGIM has no specific parameter settings. In SRIM, the parameters are set as l = l 0 jjdjj 2 , where l 0 = 1.8 can achieve the best results. The parameters of ESRIM are set as a 0 = 2.1 and s 0 = 0.5. The simulation results are shown in Figures 8  and 9.
As with straight lines, tracking trajectories of spatial curves are also among the most common basic tasks.
Taking arc trajectory L C1 as an example, the optimization results with SGIM, SRIM, and ESRIM are studied. As compared with straight lines, tracking spatial curves are more difficult. From Figure 8, it can be seen that the singularity of SGIM occurs shortly after the simulation begins. At the singular point, the terminal error breaks through the allowed limit and the tracking fails. SRIM with a damping factor passes through the initial singularity, but the terminal error oscillates severely  during the iteration process. When approaching the next singularity point, the terminal error experiences a sharp divergent oscillation and this leads to tracking failure. However, ESRIM optimization exhibits the best result. Generally, ESRIM obtains lower joint angular velocities and larger singular values, corresponding to better practicability and flexibility. As with the results in Example 1, the terminal error of ESRIM is sometimes higher than SGIM and SRIM, but the error still meets the accuracy requirement. Moreover, only ESRIM passes through both two singular points successfully and obtains joint angle solutions during the whole curve-tracking process. As shown in Figure 9, only the Lyapunov derivative value of ESRIM is negative and satisfies the stability requirement. It can be seen that ESRIM again exhibits better optimization performance in this curve trajectory tracking for the 7-axis manipulator.

Example 3: 10-axis manipulator
In this example, the 7-axis is replaced with a multi-axis redundant manipulator to validate the performances of the different methods. The DH structure parameters of the 10-axis manipulator are listed in Table 2.
The optimization is applied to spatial arc curve L C2 , the RPY attitude angle extends from S a (84.5°, 84.0°, and 90.4°) to E a (120.6°, 40.9°, and 140.8°), and other parameters are the same as those described in Example 2. The arc-length step size is also set to 0.05 mm. The parameters of SRIM are set as l 0 = 0.01 and the parameters of ESRIM are a 0 = 1.7 and s 0 = 1.3. The simulation results are shown in Figures 10 and 11.
As with the 7-axis manipulator, the SGIM solution for the 10-axis manipulator also becomes singular near the beginning of the simulation. SRIM fails to perform better than does SGIM in the tracking process for trajectory L C2 , while the proposed ESRIM performs well. In the limited range, the terminal error is smooth and continuous, and no overall oscillation occurs. The least singular value at the same position is the largest while the joint angular velocity is the smallest. Therefore, the ESRIM shows better optimization performance for the multi-axis redundant manipulator. In addition, the Lyapunov derivative value of ESRIM is negative, which is as expected.

Conclusion
The conclusions of this study are as follows: 1. To address the problem that the manipulator's singularity of different joints cannot be adjusted separately by the classical SRIM, a scheme in which the damping factor matrix b replaces the   unified damping coefficient l is proposed. To address the Jacobian matrix approximation error, which usually exists in the iterative control optimization process, the correction coefficient matrix a is employed and the ESRIM is formed. 2. Large values of manipulator terminal error and joint angular velocity, caused by the singularities in three optimization methods are analyzed, and the semi-empirical selection principle of b and a is given. Furthermore, the stability of the proposed method is investigated with reference to the Lyapunov stability criteria. 3. Three simulation methods are compared in the tracking of spatial straight lines and curved trajectories. The simulation results show that the ESRIM exhibits better optimization performance and stability as the process continues. In the prescribed range, the terminal error does not show oscillation with longer duration or large amplitude. Meanwhile, under the same conditions, the least singular value of ESRIM is the largest and the joint angular velocity is the smallest, which is as expected. 4. The different optimized methods have also been applied with a 10-axis manipulator, and positive results are obtained with ESRIM; this indicates that the extended robust inverse method for multi-axis redundant manipulators has broad applicability and can avoid singularity. However, it should be noted that all three methods based on singular value decomposition of the Jacobian matrix require long calculation times, which is deserved to research and this should be investigated and improved in future work.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.