Singular Configuration Analysis and Singularity Avoidance with Application in an Intelligent Robotic Manipulator

Abstract

Recently, robotic sensor systems have gained more attention annually in complex system sense strategies. The robotic sensors sense the information from itself and the environment, and fuse information for the use of perception, decision, planning, and control. As an important supplement to traditional industrial robots, co-bots (short for co-working robots) play an increasingly vital role in helping small and medium-sized enterprises realize intelligent manufacturing. They have high flexibility and safety so that they can assist humans to complete highly repetitive and high-precision work. In order to maintain robot safe operation in the increasing complex working environment and human–computer intelligent interactive control, this paper is concerned with the problem of applicant accuracy analysis and singularity avoidance for co-bots. Based on the dynamic model with load and torque sensors, which is used to detect the external force at the end of the robot, this paper systematically analyzes the causes of singularity phenomenon in the robot motion control. The inverse solution is obtained by analytical method and numerical method, respectively. In order to ensure the smooth and safe operation in the whole workspace, it is necessary for a robot to avoid singularity. Singularity avoidance schemes are utilized for different control tasks, including point-to-point control and continuous path control. Corresponding simulation experiments are designed to verify the effectiveness of different evasion schemes, in which the advantages and disadvantages are compared and analyzed.

1. Introduction

Recently, appliance accuracy analysis and human–computer interaction have undergone rapid development in intelligent assembly machine factories, with typical examples including assembly, polishing, dual-arm coordination, or dexterous hand manipulations [1,2,3]. In order to adapt to small-batch, customized, and short-cycle production, it is necessary for co-bots, a new type of industrial robot that can work with human beings in the same working environment, to develop a hierarchical control algorithm that enables safe and stable cooperative locomotion. With dynamic model and torque sensors, the robot is achieving the drag control by compensations of gravity and friction, which counteract external disturbances. However, it will inevitably lead to the complexity and uncertainty of environment [4,5]. Reaching the singular region makes it easy for a robot to cause various problems, such as instability, poor performance, and so on. Emphasis is placed on the co-bot singularity. Hence, the singularity avoidance is the basis to ensure co-bot stable operation in the man–machine cooperation.

The singular configuration is a phenomenon which hampers the motion of the robot end effector [6,7]. The freedom of the robot end will decrease in a singular configuration, so it cannot be controlled to move in directions. While the robot approaches the singular adjacent region, some joints calculated by the inverse kinematic tend to infinity, which will cause the decline of tracking expected trajectory [8]. The singularity is the inherent characteristic of an articulated robot, which mainly appears in finding the inverse kinematics solution. The actual performance is that the calculated expected joint speed is much greater than the maximum joint speed which the actual motor can provide. Bohigas [9] gave the general algorithm of robot singular configuration with different mechanical structures. Müller [10] analyzed singular points by using Lie groups separated from the specific structure.

For more general singularity avoidance problems, Maciejewski [11] proposed a damped least square method with a damping term near singular region to obtain approximate solution of joint velocity. The essence lies in limiting the joint velocity near the singularity by adding a damping coefficient. Serial robots conduct seam welding, sealant application, spray painting, polishing, deburring, or other tasks requiring uninterrupted continuous motion. Chiaverini [12] introduced a continuous nonlinear function to characterize the damping factor, which is used to ensure the continuity of joint velocity near singular configuration. However, it will increase the tracking error of robot end in all directions. The damping factor is only added to the minimum singular value of a Jacobian matrix. In addition, Megalingam et al. [13] proposed singular separation to determine the singular direction of Kinova robots, and then avoid joint movement to pass through singular points. References [14,15] also proposed similar analysis methods to analyze singular directions, while these methods are not universal. A large number of SVD decomposition operations are required in the above singular point avoidance methods. In order to reduce amount of computation, Xu et al. [16] proposed the method of “singular separation + damping reciprocal” for PUMA robots. It was the first time to separate the singular problem into two parts of position and attitude. The redundant robots have infinite sets of inverse kinematics solutions, hence the optimal solution (i.e., the solution furthest away from the singular region) can be found by designing some optimization objective. Wang et al. [17] summarized some methods of using redundancy to avoid the singular region with series robot. The ultimate goal is to find the solution of the singular configuration of the robot principle in countless redundant solutions. The optimization directions include the determinant and condition number of Jacobian matrixes. However, the calculation using Jacobian matrix is too complex, and an optimization is proposed in reference [18,19].

All the singularity avoidance methods above are at the cost of tracking accuracy. The parameters of most methods are closely related to experience, especially the determination of damping factor and singular region. In addition to singularity avoidance method satisfying trajectory tracking, Taki et al. [20] proposed a method strictly tracking the desired path. Different from damping avoidance method, the avoidance algorithm proposed is to reduce running speed by reducing the expected speed when approaching the singular region, so as to avoid limiting the motion of the robot end effector caused by singular configuration. However, the strange avoidance method will only produce partial errors in tracking directions.

All these observations motivate the current study. In this paper, we concentrate on singular configuration analysis and singularity avoidance of co-bot, taking the JACO2 robot of Kinova Robotics company as an example [4,13]. It has seven DOF (degrees of freedom), with weight of 5.5 kg, a load of 1.5 kg, and a maximum arm span of 98.4 cm. Each joint motor is equipped with a torque sensor. The rotating potentiometer detects robot joint rotation, and the force sensor reflects ground reaction force information. The main contributions can be summarized as follows.

• The D-H modeling method is utilized to build a joint coordinate system of the JACO2 robot. The inverse kinematics is solved by analytical method and numerical method, and the operation speed and accuracy of two schemes are compared.
• The singularity caused by inverse kinematics is analyzed, and the robot singular configuration conditions are based on the block analysis of wrist Jacobian matrix, which are divided into three singular types: internal singularity, external singularity, and wrist singularity.
• According to different task requirements, singularity avoidance schemes are utilized in robot redundancy, damped least squares, and singularity consistency. It is noted that these works focus on the working principle and control methods of robotic sensor system.

2. Materials and Methods

2.1. D-H Model

The kinematics model of the JACO2 robot with D-H model is illustrated in Figure 1, and its parameters are shown in

Table 1. D-H parameters of the JACO2 robot.

$i$	${\alpha }_{i-1}$ (rad)	$a_{i-1}$ (m)	$d_i$ (m)	${\theta }_i$(rad)
1	$\mathsf{\pi }/$2	0	−0.2755	${\theta }_1$
2	$\mathsf{\pi }/$2	0	0	${\theta }_2$
3	$\mathsf{\pi }/$2	0	−0.41	${\theta }_3$
4	$\mathsf{\pi }/$2	0	−0.0098	${\theta }_4$
5	$\mathsf{\pi }/$2	0	−0.3111	${\theta }_5$
6	$\mathsf{\pi }/$2	0	0	${\theta }_6$
7	$\mathsf{\pi }$	0	−0.2638	${\theta }_7$

.

According to D-H parameters, the transformation matrix between each links is

$${}^{\mathrm{base}}T{}_0=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& -1& 0& 1\\ 0& 0& -1& 0\\ 0& 0& 0& 1\end{array}\right],{}^{i}T{}_{i+1}=\left[\begin{array}{cccc}cos{\theta }_{i+1}& 0& sin{\theta }_{i+1}& 0\\ sin{\theta }_{i+1}& 0& -cos{\theta }_{i+1}& 0\\ 0& 1& 0& d_{i+1}\\ 0& 0& 0& 1\end{array}\right], {}^{6}T{}_7=\left[\begin{array}{cccc}cos{\theta }_7& sin{\theta }_7& 0& 0\\ sin{\theta }_7& -cos{\theta }_1& 0& 0\\ 0& 1& -1& d_7\\ 0& 0& 0& 1\end{array}\right],$$

where $i=0,1,2\dots 5$. Thus, the forward kinematics equation of the JACO2 robot can be expressed as

$${}_{}{}^{base}T{}_n={}_{}{}^{base}T{}_0{}_{}{}^{0}T{}_1{}_{}{}^{1}T{}_2{}_{}{}^{2}T{}_3{}_{}{}^{3}T{}_4{}_{}{}^{4}T{}_5{}_{}{}^{5}T{}_6{}_{}{}^{6}T{}_7=f\left(\mathrm{\Theta }\right)$$

(1)

where ${}_{}{}^{i}T{}_j$ denotes the homogeneous transformation matrix from $\left\{i\right\}$ to $\left\{j\right\}$ coordinate system, and $\mathrm{\Theta }=\left\{{\theta }_1,{\theta }_2\cdots {\theta }_n\right\}$ denotes angle information of each joint. The Euler angle (x-y-z)$\left[\alpha ,\beta ,\gamma \right]$ and $p={\left[p_x,p_v,p_z,\alpha ,\beta ,\gamma \right]}^T$ denote attitude of robot end and end pose, respectively, where $\left[p_x,p_y,p_z\right]$ denote robot end position.

2.2. Inverse Kinematics

The analytical solution and numerical solution are generally adopted to resolve robot inverse kinematics. The former uses a robot’s geometric configuration by separating the forward kinematics parameters. The latter is based on numerical iteration, which obtains the optimal solution by setting objective function and controlling joint angle to move to the opposite gradient direction.

The JACO2 robot is a seven-DOF spherical revolution spherical configuration series robot. For a certain end pose, the robot has countless sets of inverse kinematics solutions. The redundancy makes the JACO2 robot operate more flexible than a traditional six-DOF series robot but increases the computational complexity. Two methods are used to solve inverse kinematics, respectively.

• Analytical solution

The JACO2 robot has one more rotating joint at third joint and one more offset at fourth joint. Taking the third joint angle ${\theta }_3$ as the redundant parameter, other joint angles are

$${\theta }_i=f\left({\theta }_3\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,2,I,7$$

(2)

All inverse kinematics solutions can be obtained by traversing all effective value, and only one set of optimal solution can be selected according to cost function. Therefore, the inverse kinematics optimization problem is simplified as

$$\begin{array}{l}\min f\left(\mathrm{\Theta }\right)\\ s.t.\hspace{0.17em}\hspace{0.17em}{}^{base}T\left(\mathrm{\Theta }\right)={}^{base}T{}^d\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\mathrm{\Theta }\in \left[{\mathrm{\Theta }}^{\min },{\mathrm{\Theta }}^{\max }\right],\end{array}$$

(3)

where $\hspace{0.17em}\mathrm{\Theta }={\left[{\theta }_1\hspace{0.17em}\hspace{0.17em}{\theta }_2\hspace{0.17em}\hspace{0.17em}{\theta }_3\hspace{0.17em}\hspace{0.17em}{\theta }_4\hspace{0.17em}\hspace{0.17em}{\theta }_5\hspace{0.17em}\hspace{0.17em}{\theta }_6\hspace{0.17em}\hspace{0.17em}{\theta }_7\right]}^T$ is angle vectors of each joint, ${}_{}{}^{base}T\left(\mathrm{\Theta }\right)$ is the forward kinematics equation, ${}_{}{}^{base}T{}^d$ is the given end pose, and $\left[{\mathrm{\Theta }}^{\min },{\mathrm{\Theta }}^{\max }\right]$ represents the constraint of each joint.

Substituting (2) into (3), the inverse kinematics optimization problem is turned into a one-dimensional optimization of ${\theta }_3$:

$$\min \hspace{0.17em}f\left({\theta }_3\right)\hspace{0.17em}\hspace{0.17em}s.t.\hspace{0.17em}\mathrm{\Theta }\in \left[{\mathrm{\Theta }}^{\min },{\mathrm{\Theta }}^{\max }\right].$$

(4)

Selecting the motion amplitude of each joint as the cost function,

$$f\left(\mathrm{\Theta }\right)={{\displaystyle \sum }}_{i=1}^7{\lambda }_i{\left({\theta }_i^{t+1}-{\theta }_i^t\right)}^2,$$

(5)

where ${\lambda }_i$ is the weight coefficient, ${\theta }_i^{t+1}$ is the solution of the i-th joint in the current period, and ${\theta }_i^t$ is the solution of the i-th joint in the previous cycle. A greedy search is used to find the optimal solution.

2.

Numerical solution

The numerical method focuses on the equation relationship instead of robot geometric configuration. Given this, $f\left(\mathrm{\Theta }\right)={}_{}{}^{base}T{}_n$, expected position ${}_{}{}^{base}T{}_d$, which are to be solved with constraints $\mathrm{\Theta }\in \left[{\mathrm{\Theta }}^{\min },{\mathrm{\Theta }}^{\max }\right]$. The Jacobian matrix describes the relationship between robot end velocity and joint velocity:

$$\dot{x}=J\dot{q}$$

(6)

where $\dot{x}={\left[v_x\hspace{0.17em}{v}_y\hspace{0.17em}{v}_z\hspace{0.17em}{w}_x\hspace{0.17em}{w}_y\hspace{0.17em}{w}_z\right]}^T$ represents terminal velocity. $R_{3\times 3}$ describes the relationship between Euler angle and terminal angular velocity.

$$\dot{p}=\left[\begin{array}{c}{I}_{3\times 3}\\ R_{3\times 3}{}^{-1}\end{array}\right]\hspace{0.17em}\dot{x}=\left[\begin{array}{c}{I}_{3\times 3}\\ R_{3\times 3}{}^{-1}\end{array}\right]\hspace{0.17em}J\dot{q}=J’\dot{q}, R_{3\times 3}=\left[\begin{array}{c}1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}sin\beta \\ 0\hspace{0.17em}\hspace{0.17em}cos\alpha \hspace{0.17em}\hspace{0.17em}-sin\alpha cos\beta \\ 0\hspace{0.17em}\hspace{0.17em}sin\alpha \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}cos\alpha cos\beta \end{array}\right]$$

(7)

The pseudocodes of analytical and numerical inverse kinematics solutions are shown in Algorithms 1 and 2.

Algorithm 1: According to the Cost Function (7), the Obtained Optimal Solution  ${\theta }_3$
Input:$T^d$ (End pose to be solved)
${\mathrm{\Theta }}^0$ (Current joint angle vector)
$n$ (Iterative numbers per layer loop)
$\epsilon$ (Minimum allowable error)
$\left[{\theta }_3^{\min },{\theta }_3^{\max }\right]$ (Possible value range)
Output:${\mathrm{\Theta }}^{opt}$ (Optimal joint angle vector)
1:	do
2:	$r=\left({\theta }_3^{\max }-{\theta }_3^{\min }\right)/\left(n-1\right)$ (Search step of ${\theta }_3$)
3:	$S=\left\{{\theta }_3^{\min }+\left(m-1\right)r|m=\mathrm{I}\dots n\right\}$
4:	$C_{\max }=inf$ (Cost)
5:	for $r$ in $S$ do
6:	${\theta }_3=r$
7:	$\mathrm{\Theta }=InvKin\hspace{0.17em}\left(T^d,{\theta }_3\right)$ (Given ${\theta }_3$ to find inverse kinematics solutions. If the
	solution is within the joint limit, $\mathrm{\Theta }$ is returned; otherwise, $\Phi$ is returned.)
8:	if $\mathrm{\Theta }\ne \Phi$ and $f\left(\mathrm{\Theta }\right)<C_{\max }$ then
9:	$C_{\max }=f\left(\mathrm{\Theta }\right)$
10:	${\mathrm{\Theta }}^{opt}=\mathrm{\Theta }$
11:	else
12:	continue
13:	end
14:	${\theta }_3^{\min }={\theta }_3^{opt}-s$
15:	${\theta }_3^{\max }={\theta }_3^{opt}+s$
16:	while $\gamma >\epsilon$

Algorithm 2:$\phi$is the Partial Derivative of the Joint Angle
Input:$P^d$ (Terminal position vector needed to be solved)
${\mathrm{\Theta }}^0$ (Current joint angle vector)
$\lambda$ (Step size per iteration)
$n_{\max }$ (Maximum number of iterations)
$\epsilon$ (Maximum allowable error)
Output:${\mathrm{\Theta }}^{opt}$ (Optimal joint angle vector)
1:	$\mathrm{\Theta }={\mathrm{\Theta }}^0$,$n=1$
2:	$err=P^d-f\left(\mathrm{\Theta }\right)$
3:	do
4:	$\Delta \mathrm{\Theta }=J^{+}\cdot err+\left(I_{7\times 7}-J^{+}J\right)\hspace{0.17em}\phi$ ($\phi$ is used to control the iteration direction)
5:	$\Delta \mathrm{\Theta }=normalize\left(\Delta \mathrm{\Theta }\right)$ (Normalized to [-pi,pi])
6:	$\mathrm{\Theta }=\mathrm{\Theta }+\lambda \cdot \Delta \mathrm{\Theta }$
7:	$err=P^d-f\left(\mathrm{\Theta }\right)$
8:	$n=n+1$
9:	while $\Vert err\Vert <\epsilon$ or $n>n_{\max }$
10:	${\mathrm{\Theta }}^{opt}=\mathrm{\Theta }$

2.3. Singular Configuration Analysis

Both of the two methods illustrated in Section 2.2 will introduce singularity. Whether the robot is in a singular configuration depends on the full rank of the Jacobian matrix,

$$J=\left[\begin{array}{c}{J}_p\\ J_{\omega }\end{array}\right]=\left[\begin{array}{c}v\\ \omega \end{array}\right]=\left[\begin{array}{c}{e}_1\times r_1\hspace{0.17em}\hspace{0.17em}{e}_2\times I\hspace{0.17em}\dots \hspace{0.17em}{e}_7\times r_7\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{e}_1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{e}_2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{e}_7\hspace{0.17em}\end{array}\right]$$

(8)

where $e_i$ is the unit vector of the i-th joint, and $r_i$ is the vector between the origin and end position. The Jacobian matrix at the end and wrist satisfies

$$J_e=\left[\begin{array}{c}{I}_{3\times 3}\hspace{0.17em}\hspace{0.17em}{U}_{3\times 3}\\ 0_{3\times 3}\hspace{0.17em}\hspace{0.17em}{I}_{3\times 3}\hspace{0.17em}\end{array}\right]\hspace{0.17em}{J}_w, U_{3\times 3}=\left[\begin{array}{c}\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-d_7{c}_6\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-d_7{s}_5{s}_6\\ \hspace{0.17em}{d}_7{c}_6\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{d}_7{c}_5{s}_6\\ d_7{s}_5{s}_6\hspace{0.17em}\hspace{0.17em}-d_7{c}_5{s}_6\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\end{array}\right]\hspace{0.17em}.$$

(9)

Let $c_i$ and $s_i$ be abbreviations of $cos{\theta }_i$ and $sin{\theta }_i$, respectively, and that $det(J_e{J}_e^T)=det(J_w{J}_w^T)$, $Rank\left(J_e\right)=Rank\left(J_w\right)$. According to (8), the Jacobian matrix of wrist is

$$J_w=\left[\begin{array}{c}{J}_{11\hspace{0.17em}3\times 4}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{0}_{3\times 3}\\ J_{21\hspace{0.17em}3\times 4}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{J}_{22\hspace{0.17em}3\times 3}\end{array}\right].$$

(10)

It is expressed as (11) in the wrist coordinate system:

$${}_{}{}^{w}J{}_w={\left[\begin{array}{c}{}_{}{}^{base}R{}_w^{-1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{0}_{3\times 3}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{0}_{3\times 3}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{}_{}{}^{base}R{}_w^{-1}\end{array}\right]}^{base}{J}_w$$

(11)

where ${}_{}{}^{base}R{}_w$ represents wrist Jacobian matrix.

• Positional singularity

The unfilled rank of the first three rows of the Jacobian matrix will result in position singularity, which can be described as $det({}_{}{}^{w}J{}_{11}{}_{}{}^{w}J{}_{11}^T)=0$,

$$det({}_{}{}^{w}J{}_{11}{}_{}{}^{w}J{}_{11}^T)={{\displaystyle \sum }}_{i=1}^4{M}_i^2$$

(12)

where $M_{i }$ represents the determinant of cofactor of matrix ${}_{}{}^{w}J{}_{11}$.

$$\begin{array}{l}{M}_1=det\left(\left[j_1\hspace{0.17em}{j}_2\hspace{0.17em}{j}_3\right]\right)=0,\\ M_2=det\left(\left[j_1\hspace{0.17em}{j}_2\hspace{0.17em}{j}_4\right]\right)=d_3{d}_5{s}_4\left(d_3{s}_2-d_5{s}_2{c}_4+c_2{d}_5{c}_3{s}_4+c_2{d}_4{s}_3\right),\\ M_3=det\left(\left[j_1\hspace{0.17em}{j}_3\hspace{0.17em}{j}_4\right]\right)=d_3{d}_5{s}_2{s}_4\left(c_3{d}_4-d_5{s}_3{s}_4\right),\\ M_4=det\left(\left[j_2\hspace{0.17em}{j}_3\hspace{0.17em}{j}_4\right]\right)=d_3{d}_5{s}_4\left(d_4{s}_3+c_3{d}_5{s}_4\right).\end{array}$$

The conditions in whether each co-factor equals zero are expressed as

Table 2. The condition in which the determinant of the ${}_{}{}^{w}J{}_{11}$ cofactor of each order is 0.

Number	Condition	M1	M2	M3	M4
1	$s_4=0$	√	√	√	√
2	$d_3{s}_2-d_5{s}_2{c}_4+c_2{d}_5{c}_3{s}_4+c_2{d}_4{s}_3=0$	√	√
3	$s_2=0$	√		√
4	$c_3{d}_4-d_5{s}_3{s}_4=0$	√		√
5	$d_4{s}_3+d_5{s}_3{s}_4=0$	√			√

.

When $s_4=0$, the ${\theta }_4$ joint limitation is 32°~328°, the robot arms stretch out, and the end is at the boundary point, which is called boundary singularity. The ${\theta }_2$ joint limitation is 49°~311°, when $d_4{s}_3+d_5{s}_3{s}_4=0$, the JACO2 robot is in a singular state, which is called internal singularity, which is shown in Figure 2. The singularity condition is illustrated in

Table 3. The position singularity condition.

Number	Condition Combination in Table 2 (Serial Number)	Simplified Conditions
1	1	$s_4=0$
2	2 + 3 + 5	$s_2=0$$, d_4{s}_3+d_5{s}_3{s}_4=0$
3	2 + 4 + 5	Not established (4 and 5 cannot be met at the same time).

.

2.

Attitude singularity

Similarly, the attitude singularity can be judged by analyzing whether the last three rows of the wrist Jacobian matrix are less than the rank. Substituting each item of Jacobian matrix, the robot attitude singularity condition is shown in

Table 4. The robot attitude singularity condition.

Number	Attitude Singularity Satisfying Condition	Note
1	$s_6=0 s_2=0$
2	$s_6=0$$s_4=0$	The position singularity is also satisfied.
3	$s_6=0$$d_4{s}_5+d_3{s}_4{c}_5=0$

, and its wrist singularity condition is shown in Figure 3.

$$det({}_{}{}^{w}J{}_{22}{}_{}{}^{w}J{}_{22}^T)=s_6^2=0, det({}_{}{}^{w}J{}_w{}_{}{}^{w}J{}_w^T)=2d_3^2{d}_5^2{s}_2^2{s}_4^2{\left(d_4{s}_5+d_3{s}_4{c}_5\right)}^2=0.$$

3. Theoretical Analysis

The main purpose of singularity avoidance is to keep the stable, continuous, and bounded running speed of each joint when the robot is in the singular region. The relationship between joint velocity and end velocity is described as

$$\dot{q}=J^{-1}\dot{x}$$

(13)

Since the Jacobian matrix of the JACO2 robot is a nonsquare matrix, (13) is turned into

$$\dot{q}=J^{+}\dot{x}+\alpha \left(I-J^T{J}^{+}\right)\hspace{0.17em}\nabla \phi$$

(14)

where $J^{+}=J^T{\left(J{J}^T\right)}^{-1}$ represents the pseudo-inverse of the Jacobian matrix $J$, $\alpha$ is coefficient, and $\nabla \phi$ is an arbitrary vector. In order to ensure that the joint velocity is bounded, it is necessary to change Jacobian matrix $J$ or end running velocity $\dot{x}$.

3.1. Avoidance Strategy with End Velocity Constant

• Internal singularity

To keep the end velocity constant, according to (15), only the Jacobian matrix J needs to be changed. The cost function (10) is introduced to obtain the optimal solution. The function is designed:

$$f\left(\mathrm{\Theta }\right)=\frac{1}{\sqrt{det\left(J{J}^T\right)}}$$

(15)

where $J$ respresents Jacobian matrix corresponding to each group of possible inverse kinematics solutions. Then, $det\left(J{J}^T\right)$ will be maximized to avoid singularity. The cost function can be designed as

$$f\left(\mathrm{\Theta }\right)=\frac{1}{{\sigma }_r^2}$$

(16)

$$f\left(\mathrm{\Theta }\right)=\frac{{\sigma }_r^2}{{\sigma }_1^2}$$

(17)

where ${\sigma }_r$ and ${\sigma }_1$ represent minimum and maximum singular value, respectively. The convergence direction can be controlled by selecting an appropriate vector $\nabla \phi$, which controls the gradient of cost function with respect to joint angle $\mathrm{\Theta }$. Hence, the equivalent effect of avoiding singularity with the analytical rule is

$$\nabla \phi =\frac{\partial f\left(\mathrm{\Theta }\right)}{\partial \mathrm{\Theta }}$$

(18)

However, considering that all redundant solutions of the JACO2 robot make the manipulator elbow operate on a spatial circle, the above avoidance method can only have a good avoidance effect on the internal singularity and yet have nothing to do with external singularity and pose singularity.

2.

Other singular configurations

Other singular configurations cannot be simply avoided by redundant solutions. When the robot elbow is extended, which belongs to external singularity, the robot end is close to the workspace boundary. The robot cannot find the solution far away from the singularity (the fourth joint angle of all solutions approaches π). For six-DOF series industrial robots (most are PUMA type), the method of “singular separation + damping coefficient” adopted in literature [16,21] can obtain good results and analyze the corresponding errors, while it is not suitable for the JACO2 redundant robot, since singular separation effect cannot be realized directly. A damped least squares method is adopted for singular avoidance in this section, which combines the merit of Newton and steepest descent method. It not only ensures the convergence of iterative calculation, but also speeds up the convergence speed.

$$\dot{q}=J^{\#}\dot{x}$$

(19)

where $J^{\#}$ is the pseudo-inverse of damped Jacobian matrix,

$$J^{+}=J^T{\left(J{J}^T+{\lambda }^{2}I\right)}^{-1}$$

(20)

Then, minimizing the cost function, and regularizing velocity term by coefficient,

$$f\left(\mathrm{\Theta }\right)={{\displaystyle \sum }}_{i=1}^7{\left({\theta }_i^{t+1}-{\theta }_i^t\right)}^2+{{\displaystyle \sum }}_{i=1}^7{\lambda }^2{\left({\dot{\theta }}_i^{t+1}\right)}^2$$

(21)

Then, we regularize of velocity term by coefficient. The minimum singular value of the Jacobian matrix is used as the criterion to judge whether the robot enters the singular region. Equation (20) can be transformed into

$$J={{\displaystyle \sum }}_{i=1}^6{\sigma }_i{u}_i{v}_i$$

(22)

$$J^{\#}={\displaystyle \sum }_{i=1}^6\frac{{\sigma }_i}{{\sigma }_i^2+{\lambda }^2}{v}_i{u}_i^T$$

(23)

$$\dot{\mathrm{\Theta }}=J^{\#}\dot{x}={\displaystyle \sum }_{i=1}^6\frac{{\sigma }_i}{{\sigma }_i^2+{\lambda }^2}{v}_i{u}_i^T$$

(24)

$$\mathrm{where} {\lambda }^2=\{\begin{array}{c}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\sigma }_6>\epsilon \\ \left(1-{\left(\frac{{\sigma }_6}{\epsilon }\right)}^2\right){\lambda }_m^2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\sigma }_6\le \epsilon \end{array}$$

3.2. Avoidance Strategy with End Velocity Changed

The external singularity and wrist singularity are mainly discussed here. The avoidance strategy of end speed change is to change the expected end speed in (15). A path parameter p is introduced to describe the trajectory, and end operation trajectory can be expressed as

$$x=x\left(p\right)$$

(25)

The tracking of desired trajectory is expressed by mathematical expression,

$$h\left(\tilde{q}\right)=k\left(p\right)-x\left(p\right)=0$$

(26)

where $k\left(p\right)$ represents the forward kinematics solution process, $\tilde{q}$ is the set of path parameters p and joint angles, i.e., $\tilde{q}=\left(q,p\right)$. The tracking target will satisfy

$$dh\left(\tilde{q}\right)=H\left(\tilde{q}\right)d\left(\tilde{q}\right)=0$$

(27)

$$H\left(\tilde{q}\right)=\frac{\partial h\left(\tilde{q}\right)}{\partial \tilde{q}}=\left[J\left(q\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-S\right]$$

(28)

where $S$ is a unit vector, which satisfies $dx\left(p\right)=dpS$, and $\Vert dv\Vert =\Vert dv\left(x\right)\Vert$ direction, S is the vector of the end running speed, and parameter P is the end running speed. The general solution of (27) can be expressed as

$$d\tilde{q}=bf\left(\tilde{q}\right)$$

(29)

where b is constant, $f\left(\tilde{q}\right)\in \ker \hspace{0.17em}H\left(\tilde{q}\right)$,

$$f\left(\tilde{q}\right)={\left[f^T\left(q\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}detJ\left(q\right)\right]}^T$$

(30)

where $f\left(\tilde{q}\right)=\left(adjJ\left(q\right)\right)S\left(p\right)$. Equation (29) can be expanded to

$$dq=bf\left(\tilde{q}\right)$$

(31)

$$dp=bdetJ\left(q\right)$$

(32)

As long as the trajectory parameters are designed to satisfy (32) and give the joint speed of formula (31), the robot end trajectory can be always followed. It can be found that when the robot approaches the singular configuration, the determinant of the Jacobian matrix approaches 0, and the given constant and joint velocity both tend to infinity. Therefore, the singular phenomenon can be avoided. The constant b is determined by

$$b^{*}=\frac{d{q}_{\max }}{f\left(\tilde{q}\right)}$$

(33)

Reprogramming end velocity based on (32),

$$d{p}^{*}=b^{*}det\hspace{0.17em}J\left(q\right)$$

(34)

In that the JACO2 robot has redundant joints, (31) is deformed as

$$dq=b_{ep}{n}_{ep}+b_{sf}{n}_{sf}$$

(35)

$$dp=b_{ep}det\hspace{0.17em}{J}_3$$

(36)

where $b_{ep}$ and $b_{sf}$ are all constants, $J_3$ is the Jacobian matrix removing cofactor of column 3, $n_{ep}\in rowJ$, $n_{sf}\in \ker J$. However, when internal singularity occurs, the symbol of $det\hspace{0.17em}{J}_3$ and $dp$ will change, making the expected trajectory reverse, which is equivalent to an obstacle.

$$n_{ep}={{\displaystyle \sum }}_{i=1}^7{b}_i{n}_i$$

(37)

$$dp=b_{ep}{{\displaystyle \sum }}_{i=1}^7{b}_{i}det\hspace{0.17em}{J}_i$$

(38)

where $J_i$ is the Jacobian matrix removing the cofactor of column $i_{}$, the $i\mathrm{th}$ item of $n_i$ is always 0, then remove $n{’}_i$ item, there is $Jn{’}_i=C_{i}det\hspace{0.17em}{J}_i$. The $i\mathrm{th}$ item of $C_i$ is 1, and other items are 0.

Let $b_i=det\hspace{0.17em}{J}_i$, then $dp=b_{ep}{\displaystyle \sum }_{i=1}^7{b}_{i}det\hspace{0.17em}{J}_i=b_{ep}det\left(J{J}^T\right)$, then the symbol $dp$ will not change.

4. Experimental Results and Comparisons

To demonstrate the performance of the proposed algorithm, a series of simulations and comparisons are conducted. The simulations include a variety of avoidance strategies which basically cover most of the actual situations. Four experimental studies including inverse kinematics solutions and singularity analysis are carried out to test the feasibility of practical application. The symbolic operation toolbox in MATLAB is used for dynamic modeling.

The robot home position is selected as the starting point (joint angle [283.24; 162.71; 0.; 43.58; 265.23; 257.52; 288.14], pose position [0.2113; −0.2656; 0.5065; 1.6477; 1.1081; 0.1282]).

4.1. Two Inverse Kinematics Solutions

Let the robot end run at a speed of 0.1 m/s along the y-axis for 3 s, and control period be 0.1 s. A total of 200 trajectory middle points are inserted. The angles of each joint are calculated by two inverse kinematics solution algorithms. Figure 4 shows the trajectories of two kinematics algorithms. The red points represent the desired points, and blue and green describe trajectories with numerical and analytical solutions, respectively. Their solution errors are shown in Figure 5 and Figure 6.

4.2. Singularity Avoidance

To verify the effectiveness proposed in Section 3, singularity avoidance strategies based on robot redundancy, damped least squares, and singularity consistency are realized, respectively. Let the robot end run at the speed of 0.1 m/s along the positive direction along the z-axis for 7.5 s and then turn back, and continue to run along the negative direction along the z-axis for 7.5 s. Set a control cycle as 0.01 s. The joint state and tracking error are plotted with a sampling period of 0.1 s. The compared simulations are as follows.

• Taking no avoiding measures

Figure 7 shows the position change curve of each joint, Figure 8 shows the velocity variation curve (4.5 s~4.9 s), and Figure 9 shows the tracking errors of position and orientation (4.5 s~4.9 s).

2.

Using damping Jacobian matrix to avoid singularity

Assume that $\epsilon =0.1$ is the boundary of the singular region, and ${\lambda }_m^2=0.1$ is the damping coefficient. Figure 10 shows position change curve of each joint, Figure 11 shows the speed change, and Figure 12 shows tracking errors of position and orientation, respectively.

3.

Using singular consistency to avoid singularity

Figure 13 and Figure 14 show the position and speed change curve of each joint, respectively. Figure 15 shows the running speed of the robot end, Figure 16 shows each item changes of constant b, and Figure 17 shows the tracking errors of position and orientation.

5. Results Discussion

From the simulation results in Section 4.1, it is noted that both accuracies can meet the requirement of inverse kinematics. The average calculation time of using the numerical method to calculate the inverse kinematics solution is 0.0021 s, and using the analytical method is 0.0014 s. The calculation time of two methods is in the same order of magnitude. However, when the robot is in the singular configuration, the numerical method takes three times more time than the numerical method due to more iterations needed.

Figure 7, Figure 8 and Figure 9 show variation curves of joint position, velocity, and tracking error with time under no avoiding measures. It is obvious that the robot enters the singular region in about 4.65 s. At this time, the angle value of joint 4 is close to $\mathsf{\pi }$ and meets the condition of external singularity. Due to the ill condition of the Jacobian matrix, the obtained joint angle has a sudden change, and the joint speed is particularly huge, which is impossible for the motor to realize in practice. Meanwhile, the robot tracking error end pose becomes larger after entering the singular region, since the upper motor speed cannot be tracked in practical operation.

From Figure 10, Figure 11 and Figure 12, it can be seen that with damping Jacobian matrix (21), the joint angle changes smoothly, and the whole joint speed can be below 1.5 rad/s. After entering the singular region at about 4.65 s, the pose tracking error will also increase to a certain extent. However, the avoidance method has a certain effect on limiting joint speeds, it is too sensitive to the selection of parameters. Hence, it is necessary to design and select appropriate control parameters for different task requirements.

Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17 show variation curves of joint position, velocity, end position, and tracking error with singular consistency to avoid singularity. It can be seen that the robot enters the singular region at about 4.1 s. By changing constant value $b_i$, which makes the following det J_i smaller, the desired robot end speed dp are forced to change (Figure 15), to limit the speed of each joint (Figure 14). Compared with the singularity avoidance methods in Section 3.1, the joint operation is more stable, and the operation speed of each joint can also transit smoothly when entering into singularity region. Moreover, it will not cause an error deviating from desired trajectory. Therefore, when the robot end velocity is constant, the damping least square avoidance method is selected, while the singular consistency avoidance method is more appropriate when the end velocity changes.

6. Conclusions

The emergence of co-bots is an important supplement to traditional industrial robots, especially for new potential users such as small and medium-sized enterprises and increasingly complex control tasks in 3C electronic products, which determines that co-bots should have characteristics of safe use, simple operation, and convenient deployment. Taking the JACO2 robot as a carrier, this paper provides research on singularity analysis and avoidance strategies to realize safe control and smooth operation of a co-bot in the increasingly complex working environment and man–machine intelligent interactive control. The research results are summarized as follows.

The instability and potential safety hazards of the robot in the singular region are analyzed in detail. The conditions for singularity of the robot are deduced by means of Jacobian matrix separation. Then, the singular configuration conditions are analyzed in detail by means of block analysis of robot wrist Jacobian matrix, which can be divided into three singular types: internal singularity, external singularity, and wrist singularity. For different types of singularity, three singularity avoidance schemes are realized based on robot redundancy and different control tasks. This will ensure safe and stable operation of the co-bot in the whole workspace, and it eliminates safety problems in the man–machine cooperation.

In the future work, it is necessary to study some security issues through other sensors to ensure safety from the perspective of smooth operation, for example, replanning tasks to avoid collision, and adjusting the robot running state according to the working state of the operator. In addition, compliance control can also be studied to improve the robot adaptability with complex industrial tasks.