External force estimation for robotic manipulator base on particle swarm optimization

Abstract

The safe disposal of nuclear waste in radioactive environment urgently needs cost-effective approaches. Toward this goal, this article developed a method to external force estimation based on the identified model without force sensors. Firstly, the mathematical model including joint friction was obtained and transformed into the linear combination of unknown parameter to be estimated. Secondly, the unknown parameters were identified based on the improved particle swarm optimization algorithm, the identification procedure was implemented by optimizing the excitation trajectories to excite joint motion and sampling relevant data. Identified results were compared with the biogeography-based optimization algorithm and the cuckoo search algorithm. Then, the identified dynamic parameter was applied to external force estimation. Finally, the verification of external force estimation has been carried out using the Kinova Jaco² robot manipulator, and the experimental results showed that the external forces by the proposed method could be estimated with an root mean square error of 0.7 N.

Introduction

The safe disposal of nuclear waste is a global problem of great social significance and the global nuclear industry increasingly seeks cost-effective approaches for remote handling techniques to reduce radiation hazard. Remote robot manipulators can take the place of workers to do too hazardous or difficult tasks in radioactive environment where the use of people is not always possible. Many tasks performed in disposal of nuclear waste require remote robot manipulators to interact with the environment, such as grasping, cutting, dismantling, sorting, and segregating nuclear waste. To ensure the accurate execution of these tasks, not only the position/velocity control but the force control is also required to obtain the desired performances. ¹ Therefore, it is important to take the information of interaction forces into account for achieving safe and desired task performance. To obtain the external force information, a force sensor is used; the force sensor is often installed on the wrist of robot manipulator for measuring the external forces applied on the tip of the end effector. However, the force sensor with high accuracy and reliability is very expensive and undoubtedly increases the system cost of the robotic manipulator. Technically, the uncertain disturbances and harsh environment, such as high temperature and the high dose rates of nuclear sites, may affect the accuracy and reliability of the force sensor. Moreover, the use of force sensor not only increase the complexity of mechanical, electrical, and software design of the robot system but also have an influence on system modeling and control algorithm design. ² Considering the abovementioned drawbacks of force sensor, employing the external force estimation method has been a feasible option in radioactive environment where the dose rates limit the use of force sensor.

The external force estimation methods developed to estimate the forces applied on robot end effector require parameter identification of robot for an accurate mathematical model. Different parameter identification methods for the robot model have been suggested. Swevers et al. ³ presented an approach to design the optimized excitation trajectory, this approach was based on a finite Fourier series and guarantees periodic excitation, which allows data averaging in time domain and rejects the negative effects of noise disturbance. Atkeson et al. ⁴ developed the least square algorithm to estimate load and link inertial parameters of a PUMA 600 industrial robot. Wang et al. ⁵ built the Adaline neural network model to estimate the dynamic parameter for the reconfigurable robot. Taking into account the measurement errors, Olsen and Petersen ⁶ presented the maximum likelihood estimation (MLE) for estimating inertial parameter. This scheme can reject the influence of measurement noise remarkably, but it is easy to get stuck at locally optimal result in the global. Jin and Gans ⁷ proposed an intuitive and efficient criterion using Hadamard’s inequality to design the excitation trajectory for dynamic parameter identification. The reliability of the identified parameters has been verified on a six-degree-of-freedom (6-DOF) Staubli TX-90 robot. Wang et al. ⁸ employed deep learning to eliminate the effect of uncertainty to aid parameter estimation. Gautier et al. ⁹ proposed a closed-loop output error approach called the Direct and Inverse Dynamic Identification Models (DIDIM) for offline identification of the dynamic parameters, which only need to collect the joint torque measurement information, instead of calculating the joint velocity and acceleration by the filtered value of position measurement. In recent years, there are several novel methods for identifying the robot dynamic parameter using intelligence algorithms. Xu and Ke ¹⁰ presented an biogeography-based optimization (BBO) method to estimate the dynamic parameters by searching for global optimal solution. Ding et al. ¹¹ developed an artificial bee colony method to identify the unknown parameters of a 6-DOF industrial robot, which employed the measured joint position to derive the joint velocity, acceleration and combined it with the motor current to calculate the joint torque. Ding et al. ¹² proposed an improved cuckoo search (CS) method to obtain the parameter of ER-16 robot manipulator, while it can be easily trapped into local optimum when solving complex multidimensional parameter identification problems. Jia et al. ¹³ developed optimization design method based on multicriteria embedded by optimized excitation trajectories, and the parameter identification is solved by MLE.

On the other hand, various methods for estimating the interaction between the robot and its environment have been proposed. A robust force estimation observer using semiparametric model and disturbance Kalman filter is proposed to estimate the contact force of multilink manipulator. ¹⁴ A high-order finite-time observer is proposed to estimate the time-varying and nonconstant interaction forces in a finite time horizon. ¹⁵ An external force observer-based admittance method with time-varying delay compensation is developed for hydraulic series elastic actuators. ¹⁶ Sensorless contact force estimation method based on the weighted moving average with variable span is developed for serial manipulator. ¹⁷ Disturbance observer-based method is proposed to estimate the unknown interaction force. ^18,19 Interaction force estimation that utilizes extended state observer was presented to estimate the interaction time-varying forces exerted on a robotic manipulator. ²⁰ However, the above method may suffer from decreased estimation accuracy due to model errors.

Motivated by the need to perform tasks in disposal of nuclear waste in radioactive environment, we proposed a method to estimate the external forces applied on robot end effector without force sensors. The main contributions of this article are listed as follows: (a) Compared with two different intelligence algorithms, the superiority of the particle swarm optimization (PSO) algorithm in the application of parameter identification is validated. (b) A new method to estimate the external forces without force sensor based on the identified model together with the built-in torque sensor. (c) Experimental results carried on the Kinova Jaco² robot manipulator verify the effectiveness of the proposed method.

The remaining part of this article is organized as follows. The linear dynamical model of robot manipulator with friction is described in the second section. Parameters identification based on the PSO algorithm with excitation trajectory design and optimization is elaborated in the third section. In the fourth section, the identification procedure as well as experimental validation of external force estimation was implemented on the 6-DOF Kinova Jaco² robot manipulator. Additionally, the performance of the PSO algorithm is compared with some other intelligence algorithms. Finally, a brief conclusion is given in the final section.

Dynamic modeling

Applying the recursive Newton–Euler ²¹ or the Lagrange formulation, ²² the dynamic model of a robot manipulator consisted of n rigid links connected by n rotational joints can be formulated as

$\mathit{D}(\mathit{q})\ddot{\mathit{q}}+\mathit{C}(\mathit{q},\dot{\mathit{q}})+\mathit{G}(\mathit{q})={\mathit{\tau }}_d$

1

${\mathit{\tau }}_d=\left({{\displaystyle \mathit{\tau }}}_{d1},{{\displaystyle \mathit{\tau }}}_{d2},{{\displaystyle \mathit{\tau }}}_{d3},\dots ,{{\displaystyle \mathit{\tau }}}_{dn}\right)$

2

where $\mathit{D}(\mathit{q})$ refers to the symmetric inertial matrix, $\mathit{C}(\mathit{q},\dot{\mathit{q}})$ denotes the nonlinear Coriolis and centrifugal force vector, and $\mathit{G}(\mathit{q})$ denotes the gravity loading force vector; $\stackrel{ }{\mathit{q}}$ , $\dot{\mathit{q}}$ , and $\ddot{\mathit{q}}$ refer to the joint angle vector, velocity vector, and acceleration vector, respectively; ${\mathit{\tau }}_d$ is the actuator torque acting on joints; ${{\displaystyle \mathit{\tau }}}_{di}$ is the actuator torque of each link i.

To effectively identify the dynamic parameters, equation (1) can rearrange into the linear combination of identifiable parameters by the modified Newton–Euler method ²³

${\mathit{\tau }}_d={\mathit{\Phi }}_d(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}}){\mathit{\theta }}_d$

3

where ${\mathit{\Phi }}_d(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})$ denotes the (n × 10n) identification matrix, which can be treated as a nonlinear function depending only on the motion data; ${\mathit{\theta }}_d$ denotes the unknown dynamic parameters. The dynamical parameters of link i contained 10 parameters are given by

${\mathit{\theta }}_d=\left[{\mathit{\theta }}_{d1},{\mathit{\theta }}_{d2},\dots ,{\mathit{\theta }}_{dn}\right]$

4

$\begin{array}{c}{\mathit{\theta }}_{di}=\left[l_{xxi},l_{xyi},l_{xzi},l_{yyi},l_{yzi},l_{zzi},\right.\\ {\left.m_i{r}_{xi},m_i{r}_{yi},m_i{r}_{zi},m_i\right]}^{\text{T}}\end{array}$

5

where $\left(l_{xxi},l_{xyi},l_{xzi},l_{yyi},l_{xzi},l_{zzi}\right)$ is the inertial tensor of link i, $\left(m_i{r}_{xi},m_i{r}_{yi},m_i{r}_{zi}\right)$ denotes the first mass moments of link i; m_i is the mass of each link i. Noting that m_i is contained within $\mathit{G}(\mathit{q})$ in equation (1), and other components of dynamical parameters of link i except for m_i are contained within $\mathit{D}(\mathit{q})$ and $\mathit{C}(\mathit{q},\dot{\mathit{q}})$ in equation (1).

Due to the friction comprising a large proportion of the joint output torque, the torques caused by joint friction should be also included in dynamic model of robots in equation (3). To simplify the friction model, the robot joint friction can be regarded as the linear form of viscous friction and Coulomb friction in robotic applications. ²⁴ The friction torques are described as

${\mathit{\tau }}_f={\mathit{f}}_c\text{\hspace{0.17em} sign}(\dot{\mathit{q}})+{\mathit{f}}_v\dot{\mathit{q}}=[\text{sign}(\dot{\mathit{q}}),\dot{\mathit{q}}]\cdot [{\mathit{f}}_c,{\mathit{f}}_v]$

6

where ${\mathit{\tau }}_f$ denotes the vector of joint friction torques; f _c is the Coulomb coefficients vector; f _v is the viscous friction coefficients vector. Furthermore, in the presence of the external torques, the total joint torque vector can be expressed as

$\mathit{\tau }={\mathit{\tau }}_d+{\mathit{\tau }}_f+{\mathit{\tau }}_{\text{ext}}$

7

where ${\mathit{\tau }}_{\text{ext}}$ denotes the external torque on joints generated by the external force/torques, the integrated dynamic model can be rewritten as a linear parameterized form in the absence of external force/torques, which depends upon a set of parameters to be identified. ²⁵

$\mathit{\tau }=\mathit{\Phi }(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})\mathit{\theta }$

8

where τ refers to joint actuator torques, $\mathit{\Phi }(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})$ denotes an (n × 12n) identification matrix as described in equation (1), and θ is the (12n × 1) unknown inertial parameters vector. Thus, the unknown parameters of link i can be rewritten as ${\mathit{\theta }}_i={[{\mathit{\theta }}_{di},f_{ci},f_{vi}]}^{\text{T}}$ . Indeed, the identification matrix Φ in equation (8) is not full rank, it means only part of the inertial parameters have made a dynamic contribution. Therefore, the inertial parameters to be estimated can be simplified to the minimal set of inertial parameters with m identifiable parameters, that is, the vector of base parameters. It can be obtained by the numerical QR decomposition (QRD) and singular value decomposition (SVD). ²⁶ Equation (8) can be rewritten as

$\mathit{\tau }={\mathit{\Phi }}_m(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}}){\mathit{\theta }}_m(m\le 10n)$

9

where ${\mathit{\Phi }}_m(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})$ denotes the n × (m + 2n) vector identification matrix. ${\mathit{\theta }}_m$ represents the unknown dynamic parameters with (m + 2n) vector, which contain the vector of base parameters and the friction parameters.

Once the dynamic parameters vector ${\mathit{\theta }}_m$ was estimated by collecting the measured position and joint torque. The external torque on joints can be estimated by employing equation (3), which results from the reorganization of equations (1) and equation (7)

${\mathit{\tau }}_{\text{ext}}=\mathit{\tau }-\mathit{D}(\mathit{q})\ddot{\mathit{q}}-\mathit{C}(\mathit{q},\dot{\mathit{q}})-\mathit{G}(\mathit{q})-{\mathit{\tau }}_f$

10

The estimated external force/torques on the robot end effector in Cartesian space can be calculated by

${\mathit{\tau }}_{\text{ext}}={\mathit{J}}^{\text{T}}{\mathit{F}}_{\text{ext}}$

11

where ${\mathit{J}}^{\text{T}}$ denotes the transpose of the robot’s Jacobian matrix. ${\mathit{F}}_{\text{ext}}$ is defined as follows

${\mathit{F}}_{\text{ext}}={\left[F_x,F_y,F_z,{\tau }_x,{\tau }_y,{\tau }_z\right]}^{\text{T}}$

12

The estimated external force/torques on the robot end effector are obtained by multiplying the inverse of ${\mathit{J}}^{\text{T}}$

${\mathit{F}}_{\text{ext}}={({\mathit{J}}^{\text{T}})}^{-1}{\mathit{\tau }}_{\text{ext}}$

13

Parameters identification procedure

PSO for parameters identification

PSO is a heuristic optimization algorithm, which was inspired from birds’ flocking for exploring global optimal solution. ²⁷ PSO can be applied to estimate the dynamic parameters by initializing a set of particles randomly over the searching space to find the optimal solution with iteration. Particles update themselves by tracking two “extremum” in each iteration. The first one is individual extremes (p _Best), and the another one is global extremes (g _Best) found by the whole PSO. Figure 1 shows the illustration of the PSO algorithm to obtain the optimal solution, the particles are searching in the n-dimensional solution space, where each particle adjusts its velocity and new position in accordance with its own and companions’ experience.

Figure 1.

Geometrical illustration for PSO algorithm. PSO: particle swarm optimization.

The velocity and position of the particles are updated in each iteration according to the following equations

$\begin{array}{c}{\overrightarrow{v}}_j(t+1)=\omega \cdot {\overrightarrow{v}}_j(t)+C_1\cdot R_1\cdot \left(p_{\text{Best}}-{\overrightarrow{x}}_j(t)\right)\\ +C_2\cdot R_2\cdot \left(g_{\text{Best}}-{\overrightarrow{x}}_j(t)\right)\end{array}$

14

${\overrightarrow{x}}_j(t+1)={\overrightarrow{x}}_j(t)+{\overrightarrow{v}}_j(t+1)$

15

where ${\overrightarrow{x}}_{{ }_j}(t)$ and ${\overrightarrow{v}}_{{ }_j}(t)$ stand for separately the position of the particle and the speed of the particle, respectively; ${\overrightarrow{x}}_{{ }_j}(t+1)$ and ${\overrightarrow{v}}_{{ }_j}(t+1)$ are the particle position and the velocity for the next moment, respectively; coefficients C ₁ and C ₂ are the positive constant; $R_1,R_2\in [0,1]$ are independently random number. ω is the inertia weight coefficient in a range of.[0,1]

The inertia weight coefficient ω in equation (14) can be expressed by ω(t + 1) = αω(t) (0 < α < 1) which decrease linearly with the increase of iteration. This coefficient enables the particle swarm to search in the whole solution space with high probability at the initial stage and quickly converge to the local optimal solution. Then, the optimal solution can be searched in the local detail with the decrease of inertia weight.

To avoid the PSO algorithm falling into the local optimum, the mutation operation is carried out on the particles by referring to the concept of mutation in the genetic algorithm

$y_{jk}=\left\{\begin{array}{c}{r}_k\text{, rand()}<R\\ x_{jk}\text{, otherwise}\end{array}\right.$

16

where $x_{jk}$ , $y_{jk}$ represent the k-dimensional element of particle j, r_k is a random number in the range of the k-dimensional element, $\text{rand}()$ is a random number in the range of [0,1]; R ∈ (0,1) denote the mutation probability.

The fitness function determines whether the particle is the optimal condition in the PSO algorithm. In this work, the fitness function $\text{fit}(x_j(t))$ is defined as

$\begin{array}{c}\text{fit}(x_j(t))=\frac{1}{N}{\displaystyle \sum _{i=1}^N({\delta }_1{({\tau }_{1i}-{\tau }_{p1i})}^2+{\delta }_2{({\tau }_{2i}-{\tau }_{p2i})}^2}\\ +\cdots +{\delta }_n{({\tau }_{1n}-{\tau }_{p1n})}^2\end{array}$

where N is the length of the measured data, ${\tau }_{\xi i}$ (ξ = 1, 2,…, n) is the actual torque of the joints that need to be identified. Similarly, ${\tau }_{p\xi i}$ (ξ = 1, 2…, n) denotes the predicted torque according to equation (9). ${\delta }_{\xi }$ (ξ = 1, 2…, n) is a weighting coefficient in the range of [0,1]. The goal of the search is to get a minimum of the function.

Based on PSO for estimation of dynamical parameter, 12 parameters should be identified for each joint in the dynamic model of the manipulator, the detailed procedure for executing the PSO algorithm to estimate the dynamic parameters is described as follows:

Step 1: Randomly initializing the particle swarm. Initializing the swarm size according to the number of parameters needed to identify, randomly initializing the particle position and velocity in the range of kinetic parameters, initialize $p_{\text{Best}}$ and $g_{\text{Best}}$ .

Step 2: Updating the particle swarm according to equations (14) and (15).

Step 3: Doing the following for each particle:

Step 3.1: Obtaining the variant particle swarm $y_{(t)}$ according to equation (16).

Step 3.2: Selecting the favorable variation and carrying out the variation on the particle swarm $x\text{(}t\text{)}$ :

$x_j(t)=\left\{\begin{array}{ll}{y}_j(t),\hfill & f\left(y_j(t)\right)<f\left(x_j(t)\right)\hfill \\ x_j(t),\hfill & \text{otherwise}\hfill \end{array}\right.$

18

Step 3.3: Carrying out fitness comparison between each individual in the current particle swarm and $p_{\text{Best}}$ according to equation (17), and update $p_{\text{Best}}$ .

Step 4: Comparing the fitness of $p_{\text{Best}}$ for each particle, and update $g_{\text{Best}}$ .

Step 5: Checking the stop criterion. If the stop criterion is not matched, go to Step 2, otherwise, go to Step 6.

Step 6: The particle corresponding to the global extreme value is the optimal solution of the population, that is, the identified parameters of the robot manipulator.

The detailed PSO algorithm procedure for parameter identification is depicted in Figure 2.

Figure 2.

Flowchart of PSO for parameter identification. PSO: particle swarm optimization.

Excitation trajectory design

To decrease the influence of measurement noise in estimating ${\mathit{\theta }}_m$ from equation (9), the identification matrix ${\mathit{\Phi }}_m(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})$ should be well-conditioned. The movement of the robot manipulator along the designed excitation trajectory needs to ensure the sufficiency of excitation. Several trajectory parameterization methods have been applied to the parameter identification of robot manipulator, such as the interpolated trajectory ²⁸ and the fifth-order polynomials. ²⁹ We employed a finite Fourier series as the exciting trajectory in this work, for improving the signal-to-noise ratio of the measurement and reducing the effect on robot flexibility. ³⁰ The trajectories for joint i are parameterized as follows

${\mathit{q}}_i(t)=q_{i,0}+{\displaystyle \sum _{k=1}^N{a}_{i,k}\text{sin}\left(k{\omega }_{f}t\right)}+{\displaystyle \sum _{k=1}^N{b}_{i,k}\text{cos}\left(k{\omega }_{f}t\right)}$

19

where $q_{i,0}$ denotes the constant offset term, ${\omega }_f$ denotes the fundamental frequency of excitation trajectories; N is the number of harmonics and exciting trajectory for each joint contains $2N+1$ parameters. The coefficients $a_{i,k}$ , $b_{i,k}$ are the amplitude of the trigonometric functions.

Convergence rate of the parameter estimation and the noise rejection depend on the condition number of identification matrix, and it should be ensured that the objective function of optimization is to minimize the number of conditions. Thus, the condition number of identification matrix is taken as the fitness function in PSO algorithm, and the initial population is designed according to motion constraints of robot joints. The excitation trajectory optimization with the motion constraint can be described as ¹¹

$\left\{\begin{array}{l}\text{min}(\text{cond}(\mathit{\Phi }))\\ q_{i\text{min}}\le q_i(t)\le q_{i\text{max}}\\ {\dot{q}}_{i\text{min}}\le {\dot{q}}_i(t)\le {\dot{q}}_{i\text{max}}\\ {\ddot{q}}_{i\text{min}}\le {\ddot{q}}_i(t)\le {\ddot{q}}_{i\text{max}}\\ \left\{\psi (q(t)\right\}\subset \Omega \end{array}\right.$

20

where Φ is the identification matrix, $q_{i\text{min}}$ , $q_{i\text{max}}$ , ${\dot{q}}_{i\text{min}}$ , ${\dot{q}}_{i\text{max}}$ , and ${\ddot{q}}_{i\text{min}}$ , ${\ddot{q}}_{i\text{max}}$ are the lower and upper constraints of the joint position, velocity, and acceleration on each joints, respectively. $\left\{\psi (q(t)\right\}$ represents the set trajectories of robot end effector, Ω denotes the workspace of the robot manipulator. The excitation trajectory optimization could be solved using optimization toolbox in MATLAB.

Experimental results

The identification procedure was implemented on the 6-DOF Kinova Jaco² robot manipulator as shown in Figure 3. This robot manipulator has six rotation joints, including two base joints, one elbow joint, and three wrist joints. The joints of the robot manipulator have speed limitations at 48° s⁻¹. The position encoders and the built-in torque sensors are fixed on all joints, all of six joints can provide the joint torque as well as the measured information of joint position, velocity while working. The classic Denavit-Hartenberg (DH) parameters for the 6-DOF Kinova Jaco² robot manipulator are given in Table 1.

Figure 3.

6-DOF Kinova Jaco² manipulator. DOF: degrees of freedom.

Table 1.

The classic DH parameters for Kinova Jaco² manipulator.

Link	α i (rad)	ai (m)	di (m)	θi(rad)
1	π/2	0	0.2755	q1
2	π	0.41	0	q2
3	π/2	0	−0.0098	q3
4	π/2	0	−0.2814	q4
5	π/2	0	0	q5
6	π	0	−0.2341	q6

Figure 4 shows the optimized excitation trajectories with the base frequency of 0.1 Hz. The excitation trajectories are the parameterized five-term Fourier series given in equation (19), containing 11 optimum parameters for each Fourier series. All the optimum parameters are presented in Appendix 1.

Figure 4.

Optimized trajectory of the joints.

The optimal excitation trajectories $q_i(t)$ (i = 1, 2,…,6) for the six joints were controlled by the Kinovarobotics/matlab_Kinova-api_wrapper for the robot manipulator. The sampling frequency was set to 100 Hz, and the joint torque, position, and velocity data are recorded over 10 periods. It should preprocess the recorded data before identification due to measurement noises. A Savitzky–Golay filter was employed to smooth the raw data of joint position and velocity due to its excellent properties for smooth differentiation, and the conventional low-pass IIR Butterworth filter was chosen to filter the measured and averaged joint torques owing to no differentiation is needed. The end-effector path in Cartesian for Kinova Jaco² manipulator during the optimal trajectory is shown in Figure 5.

Figure 5.

The end-effector path in Cartesian during the excitation trajectory.

To highlight the performance of the proposed method, comparison with the BBO method ¹⁰ and the CS algorithm ¹² is performed for the same excitation trajectory. Parameter identification procedure using the above three algorithms is performed in MATLAB 2016b environment. The parameter settings for the three algorithms are given in Table 2.

Table 2.

Parameters for PSO, BBO, and CS algorithm.

Algorithm	Parameters
PSO	Population size N = 72
	Particle number = 100
	Learning factor C1 = C2 = 2.0
	Mutation probability =0.8
	Inertia weight ω = [0.4, 0.9]
BBO	Population size N = 72
	Mutation probability $m_{\text{max}}=0.00\text{5}$
	Maximal immigration rate I = 1
	Maximal immigration rateE = 1
	Habitat modification probability $P_{\text{mod}}=1$
CS	Population size N = 72
	Switching factor $P_a\hspace{1em}=0.\text{4}$
	Step size scaling factor α = 0.25
	Lévy exponent λ = 1.5

PSO: particle swarm optimization; BBO: biogeography-based optimization; CS: cuckoo search.

The max iteration number is set to 100 for all the algorithms and the results of parameter identification are obtained based on 10 independent runs. Figure 6 shows the solving process of the three algorithms, the result of fitness function decreases with the increase of iteration. When iterations reach 55, the optimization process of the PSO algorithm tends to be stable and the optimal result of fitness function is 2.7273, whereas those of BBO and CS were 27.8158 and 40.0971, respectively. The final optimization results which are calculated by equation (17) are 0.9293, 1.7815, and 1.9622 corresponding to the PSO, BBO, and CS algorithm, respectively. The results turned out that the PSO algorithm has a faster convergence rate and a better result than the BBO and CS algorithm. The identified dynamic parameters obtained from the PSO algorithm using the optimal excitation trajectories are shown in Table 3.

Figure 6.

Optimization process for the three algorithms.

Table 3.

The identified dynamic parameters of 6-DOF Kinova Jaco² manipulator.

Parameter	Joint 1	Joint 2	Joint 3	Joint 4	Joint 5	Joint 6
mi (kg)	0.7677	0.9957	0.7967	0.4174	0.4174	1.0754
mirxi (kg·m)	0.0015	0.0051	−0.0020	0.0002	0.0050	0.0049
miryi (kg·m)	0.0016	0.0669	0.0035	−0.0526	0.0023	0.0021
mirzi (kg·m)	0.0030	−0.0035	0.0026	0.0038	0.0211	−0.0007
Ixxi (kg·m2)	0.0022	0.0041	0.0028	0.7085	0.8275	0.0048
Ixyi (kg·m2)	0	0	0	0.0082	0.0083	0
Ixzi (kg·m2)	0	−0.0015	0	0.1128	0.1016	0.0003
Iyyi [kg·m2]	0.0006	0.0255	0.0027	0.7404	0.8520	0.0048
Iyzi (kg·m2)	−0.0003	0	−0.003	0.0005	0	0
Izzi (kg·m2)	0.0023	0.0217	0.003	0.1782	0.1707	0
fci	0.2415	0.3732	0.1102	0.1000	0.1200	0.1000
fvi	0.3190	0.1813	0.0730	0.1767	0.0346	0.0871

Figure 7 shows the comparison of the estimated joint torque based on the identified parameters for each joint with the actual torque (obtained with built-in torque sensors). All the estimated joint torque of the three algorithms can track the actual torque in all the joints, which implied that the estimated torques can be remodeled based on the identified parameters and the data of joint motion. It can be also observed that the estimated torques obtained with the PSO algorithm match better with the actual torque than the BBO and the CS algorithm, which highlights that the PSO algorithm has a good performance of global searching and obtained the optimal parameters for robot manipulator. It should be noted that the estimated torque not match well with the actual torque and the estimated torque errors increase during joint velocity reversal. It is mainly due to the bias of parameter estimation, the unmodeled dynamics, and the nonlinear friction which is not included in the model.

Figure 7.

The comparison of the estimated torque by the three algorithms with the actual torque: (a) torque of the joint 1, (b) torque of the joint 2, (c) torque of the joint 3, (d) torque of the joint 4, (e) torque of the joint 5, (f) torque of the joint 6.

To quantitatively evaluate the performance of the identified model by using the PSO, BBO, and CS algorithm, the root mean square (RMS) error between the estimated joint torque and the actual torque from tracking the optimal excitation trajectories are compared. The RMS error is defined as ³¹

${\epsilon }_{\text{RMS}}=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{({\tau }_{ai}-{\tau }_{ei})}^2}}$

21

where ${\tau }_{ei}$ is the estimated torque value of the ith point, ${\tau }_{ai}$ is the ith sample of the actual torque. N is the total number of sampling points. The RMS error between the estimated joint torque and the actual torque is compared in Table 4. It can be seen from the results that the identified parameters using the PSO algorithm are more accurate than those of other algorithms.

Table 4.

A comparison of RMS error between the estimated joint torque and the actual torque.

Joint	PSO (nm)	BBO (nm)	CS (nm)
1	0.57	0.66	0.72
2	1.36	1.46	1.48
3	0.57	0.81	0.83
4	0.51	0.81	0.96
5	0.37	0.48	0.44
6	0.46	0.42	0.61
Sum	3.84	4.64	5.04

PSO: particle swarm optimization; BBO: biogeography-based optimization; CS: cuckoo search.

To validate the efficiency of the identified parameters applying to external force estimation by the proposed method, the elastic impact experiment was carried out. The experimental setup was composed by the Robotiq FT 300 Force Torque Sensor and the spring with elastic constant of 1.12 N/mm, the spring was attached to Robotiq Force Sensor, as shown in Figure 8. It should make sure to keep the posture of the end effector approximately perpendicular to the horizontal plane, so that the spring maintain perpendicular to force sensor as much as possible. The manipulator was commanded to accomplish an up–down movement with the period of 5 s.

Figure 8.

Experimental setup for an elastic impact.

The comparison of the actual torques with the estimated torques based on the identified parameters by the three algorithms is shown in Figure 9. The known models can only track the movement induced by inner torque of robot manipulator during the spring compression, and yet the actual torque included the inner torque and external torque acting on the joints. Thus, the external torque can be derived from the difference of the actual torque and the estimated torque. It should be note that the main effort for the movement of end effector was realized by joints 2, 3, and 4 during compression. The external torques exerted on joints shown in Figure 9 were given by equation (10), which considered the external torque as the unknown disturbance to the robot manipulator.

Figure 9.

The comparison of the actual torques with the estimated torques: (a) the joint 2, (b) the joint 3, and (c) the joint 4.

Based on the estimated external torques of joints (see Figure 9), the estimated external forces acting on the robot end effector can be calculated by equation (13). Figure 10 shows estimated external forces and the measured forces (obtained with Robotiq FT 300 Force Torque Sensor) in the Z direction. The difference between the measured force and the estimated force in the Z direction is shown in Figure 11. As can be seen, the external force estimated by the three algorithms has the same trend as the measured force in terms of magnitude and frequency, which revealed the estimated external force has high linearity to the spring response. However, the external force estimated by the PSO algorithm matches better with the measured force in magnitude. The RMS error of the estimated external forces acting on the robot end effector by the three algorithms is compared in Table 5.

Figure 10.

The comparison of the measured force with the estimated force in the Z direction.

Figure 11.

The difference between the measured force and the estimated force in the Z direction.

Table 5.

RMS error of the estimated external forces for the three algorithms.

	PSO	BBO	CS
RMS (N)	0.73	1.10	1.69

From these results we can observe that the RMS error of the external force estimated by the PSO algorithm is smaller than those of other algorithms, which highlights the precision of the external force estimation by the proposed method.

Conclusion

In this article, we developed a new method to estimate the external forces without force sensor based on the identified model together with the built-in torque sensor. The optimal periodic excitation trajectories are designed for parameter identification, and all the unknown parameters for the model of a 6-DOF robot manipulator are identified by the proposed algorithm. The identified results showed that the proposed algorithm has fast convergence rate and higher identification accuracy compared with the BBO algorithm and the CS algorithm. Further, we conducted the elastic impact experiments to confirm the validity of the identified parameters applying to external force estimation, and the experimental results demonstrated that the proposed external force estimation method can successfully estimate the external force acting on the end effector of the robot manipulator with an RMS error of 0.7 N. This method can be not only a cost-effective alternative to the conventional force sensors but also a feasible solution for redundant force sensing in radioactive environment. In the future, we would focus on the controller design based on the estimated external force.