Model-based sensorless robot collision detection under model uncertainties with a fast dynamics identification

Abstract

This article presents a novel model-based sensorless collision detection scheme for human–robot interaction. In order to recognize external impacts exerted on the manipulator with sensitivity and robustness without additional exteroceptive sensors, the method based on torque residual, which is the difference between nominal and actual torque, is adopted using only motor-side information. In contrast to classic dynamics identification procedure which requires complicated symbolic derivation, a sequential dynamics identification was proposed by decomposing robot dynamics into gravity and friction item, which is simple in symbolic expression and easy to identify with least squares method, and the remaining structure-complex torque effect. Subsequently, the remaining torque effect was reformulated to overcome the structural complexity of original expression and experimentally recovered using a machine learning approach named Lasso while keeping the involving candidates number reduced to a certain degree. Moreover, a state-dependent dynamic threshold was developed to handle the abnormal peaks in residual due to model uncertainties. The effectiveness of the proposed method was experimentally validated on a conventional industrial manipulator, which illustrates the feasibility and simplicity of the collision detection method.

Introduction

Robot collaborating with humans in both domestic and industrial area is increasingly becoming a research focus, and much effort has been devoted to realize human–robot collaboration in which both counterparts share the same workspace to accomplish a goal via different interaction modes. ^1,2 Compared to conventional robots strictly confined within safe fence to avoid any possibility of physical interaction when robot is performing a predefined task, collaborative robots, also known as co-bots, have been regarded as a new generation of robots with the feasibility to complete the complementary task with close involvement of operator or direct physical contact with human for cooperation through high-level autonomy. The arising concern of safety is a critical challenge when a close distance between humans and robotic apparatus is typically demanded for collaboration scenario since the chance of unexpected collision significantly increases due to one’s frequent intervention into the robot’s operation space. ^{3 –5} To address the issues of safety, several related technical standards have been proposed to provide guidance for collaborative robots in order to minimize the potential injury risks caused by collisions between the robot and human operator or its surrounding environment. ^{6 –8}

Generally speaking, effective schemes dealing with collision can be divided into pre-collision and post-collision strategies based on when the active actions are given with respect to time of impact. On the one hand, the former strategies can usually achieve desirable collision-free performance by actively providing appropriate collision avoidance path using additional sensors to perceive unstructured surrounding environment, which is sometimes sophisticated and probably computational heavy. ^{9 –12} For example, human–robot coexistence can be guaranteed by introducing preventive approaches aiming at monitoring the operator’s distance to robot body with stereo vision ¹⁰ and proximity sensors. ¹² On the other hand, post-collision strategies give a reactive response to the impact without the requirements of real-time locating operator’s position and a complex path planning algorithm to avoid obstacles. ^{2,13 –18} Meanwhile, the algorithmic simplicity of post-collision scheme is always at the cost of the potential collision impact due to its intrinsic property that collision could not be detected until the impact between the robot and its surroundings occurs, which highlights the significance of how to accurately diagnose the abnormal signals resulting from the unexpected collision. However, the post-collision scheme, which succeeds in minimizing the unintended collision threats, is attractive to reduce the costs of external sensing and potential failures due to the integration of additional perceptive components, making the approach preferable for cost-sensitive industrial applications. ^2,19,20 In human–robot collaboration context, a comprehensive risk assessment is particularly required for robotic system adopting the post-collision strategy to embed safety scheme such as collision detection and reaction into the control system. Crash testing in robotics by Haddadin et al. experimentally manifested how the factors like robot mass and velocity affect the injury severity index from both biomechanics and forensics aspects through human–robot impact voluntary testing. ^4,21 It was demonstrated that a lightweight manipulator with compliant actuation is less likely to cause fatal injuries. ^{22 –25} Lightweight designed robots are inherently suitable for close human–robot interaction with the desirable features of low inertia and reduced weight which allows driving torque output to decrease. Thus, the injury caused by collision is also considerably minimized in contrast to conventional heavy and stiff manipulators. ^22,23 In terms of compliant actuation, flexibility is introduced into the actuator, ²⁵ even with adjustable impedance property. ²⁶ Apart from the human-friendly design with compliant and lightweight property which does not necessarily ensure the operator’s safety, control for safe physical interaction allows the realization of collision handling with the primary motivation of detecting the impact-deduced signal as fast as possible before the robot has imposed excessive force on the operator. A simple but feasible approach is to monitor non-model-based signals such as the instantaneous transients in current or the variation of position, which to some extent characterizes the system fault in control system caused by external impact. ²⁷ These approaches suffered from the drawbacks of low sensitivity and system dependence as moderate impact may produce too sight influence to observe and the detection threshold is heavily dependent on control performance. Other accessible methods principally rely on the model-based signal, which requires the complex nonlinear system dynamic description and accurate dynamic parameters provided by robot manufactures or experimentally calibrated. ^13,19,20,28 A straightforward idea is to compare the actual torques with the estimated torques, which is easily calculated using the inverse dynamics providing real-time robot motion information, to generate the residual signal as the difference of these two torques. Thus, the residual signal is attributed to a collision under the assumption of precise dynamic parameters and modeling. This method has been refined by employing an observer of the generalized momentum to avoid double numerical differentiation of position to calculate the joint acceleration, which inherently reduces the measurement noise and improves the accuracy of the residual signal. ^13,15 In a similar perspective, several disturbance observers were also discussed to estimate the exogenous disturbance imposed by the environment. ^{29 –31} In summary, the effectiveness of these model-based methods critically relies on the accuracy of relevant dynamic parameters which always keep confidential and unavailable to users and need parametric calibration by a meticulous identification experiment, since the parametric errors structurally affect the residual signal as external disturbances. The dynamic parameter identification has been extensively researched to explore the issues including symbolic calculation of base parameter, excitation trajectory design, physical feasibility, and estimation method, which requires much effort and attention to derive the sophisticated symbolic expression and enhance the accuracy with a well-conditioned excitation trajectories. ^{32 –35} One practical approach which combines machine learning method with model identification procedure enables the simplification of symbolic computation and a satisfying prediction performance. ^{36 –39} For instance, a data-driven machine learning approach was presented in literature ³⁹ to learn the dynamics of a planar two-link robot using a regressor selection technique while this method is limited to a class of planar robot whose fiction and gravity effect has little influence on the overall dynamics. Regardless of the parameter identification accuracy, the collision is eventually judged with the comparison of the residual signal and a certain threshold. A primary concern in designing threshold is how to achieve a trade-off between prevention of false alarm and detection sensitivity at the same time. It is worthwhile to tune the collision detection threshold rather than select a static threshold to achieve better detection performance. ^21,17 The generation of a dynamic threshold was elaborated in literature ¹⁷ through an adaptive residual evaluation to tackle the uncertainties in robot dynamic parameters for a two-link robot, while the online computational burden would significantly increase for robot with more degrees of freedom (DOFs), and estimation errors would accumulate as a recursive scheme was adopted in the threshold formulation.

In this article, a novel sensorless collision detection algorithm for human–robot collaboration scenario is presented on a conventional industrial robot without prior knowledge of dynamic parameters. Based on the model-based collision handing framework, the unknown dynamic parameters of the manipulator, which is indispensable to reconstruct the inverse dynamics and consequently generate the residual signal, are estimated with a sequential identification procedure. The first step of the identification experiment is to predict the gravity and friction accounting for the major components of torque when robot is running in collaboration mode. Subsequently, a machine learning approach named Lasso is introduced to identify the remaining torque, avoiding the necessity of symbolic derivation of the sophisticated expression. Furthermore, by synthesizing state-dependent dynamic threshold, the residual signal is evaluated in case of model uncertainties.

The rest of this article is organized as follows. Firstly, the dynamics modeling and model-based collision detection scheme are illustrated. Then, a fast dynamics identification is described with a sequential identification method. The gravity and friction effects are estimated with a classic method, and details are provided on the formulation of the remaining torque prediction using a machine learning approach. The residual signal is analyzed with a dynamic threshold to detect external disturbances. Subsequently, the proposed collision detection scheme is experimentally validated on a robot performing a predefined task. Finally, collusions and future directions are provided.

Model-based collision detection

Robot dynamic model

For n-DOFs serial manipulator, let ${\mathit{\theta }}_m\in {ℝ}^n$ and ${\mathit{\tau }}_m\in {ℝ}^n$ , respectively, denote the motor positions and torques, and $\mathit{\theta }\in {ℝ}^n$ and $\mathit{\tau }\in {ℝ}^n$ , respectively, denote the positions and torques after the gear transmission. These according variables are connected with the transmission matrix $R_{\text{red}}\in {ℝ}^{n\times n}$ with $\mathit{\theta }=R_{\text{red}}^{-1}{\mathit{\theta }}_m$ and $\mathit{\tau }=R_{\text{red}}{\mathit{\tau }}_m$ . However, in the presence of joint elasticity, the link positions $q\in {ℝ}^n$ are not identical to the transmitted motor positions $\mathit{\theta }$ . Stiff effects can be adopted to represent the flexible transmission model, along with additional damping effects. Without loss of generality, only the stiff model will be considered in this study to simplify the analysis. Indeed, the above positions are involved in the dynamic model for robot with elastic joints as follows

${\mathit{M}}_l(\mathit{q})\ddot{\mathit{q}}+\mathit{C}(\mathit{q},\dot{\mathit{q}})\dot{\mathit{q}}+\mathit{G}(\mathit{q})+{\mathit{\tau }}_{fl}(\dot{\mathit{q}})+\mathit{K}(\mathit{q}-\mathit{\theta })={\mathit{\tau }}_{\text{ext}}$

1

$\mathit{B}\ddot{\mathit{\theta }}+{\mathit{\tau }}_{fm}(\dot{\mathit{\theta }})-\mathit{K}(\mathit{q}-\mathit{\theta })=\mathit{\tau }$

2

where ${\mathit{M}}_l(\mathit{q})\in {ℝ}^{n\times n}$ is the positive definite inertia matrix of robot link; $\mathit{C}(\mathit{q},\dot{\mathit{q}})\dot{\mathit{q}}\in {ℝ}^n$ is the centrifugal and Coriolis torque; $\mathit{G}(\mathit{q})\in {ℝ}^n$ is the gravity item; $\mathit{B}\in {ℝ}_{}\mathrm{n\times n}$ is the constant, diagonal matrix of the motor inertia moments; $\mathit{K}\in {ℝ}^{n\times n}$ is the constant, diagonal joint stiffness matrix; ${\mathit{\tau }}_{\text{ext}}\in {ℝ}^n$ is the effective torque at joint level arising from the external impact force exerted on the robot link; and ${\mathit{\tau }}_{fl}(\dot{\mathit{q}})\in {ℝ}^n$ and ${\mathit{\tau }}_{fm}(\dot{\mathit{\theta }})\in {ℝ}^n$ are, respectively, the joint and motor friction torques. The Coulomb–viscous model is the most common model to characterize friction effect and thus the joint and motor friction torques could be depicted as ${\mathit{\tau }}_{fl}={\mathit{F}}_{vl}\dot{\mathit{q}}+{\mathit{F}}_{cl}\text{sign}(\dot{\mathit{q}})$ at joint level and ${\mathit{\tau }}_{fm}={\mathit{F}}_{vm}\dot{\mathit{\theta }}+{\mathit{F}}_{cm}\text{sign}(\dot{\mathit{\theta }})$ at motor level where $({\mathit{F}}_{vl},{\mathit{F}}_{vm})$ are viscous friction coefficients and $({\mathit{F}}_{cl},{\mathit{F}}_{cm})$ are Coulomb friction coefficients. Equations (1) and (2) are referred to as link equation and motor equation, respectively. By summarizing equations (1) and (2), we obtain

$\begin{array}{l}{\mathit{M}}_l(\mathit{q})\ddot{\mathit{q}}+\mathit{C}(\mathit{q},\dot{\mathit{q}})\dot{\mathit{q}}+\mathit{G}(\mathit{q})+{\mathit{\tau }}_{fl}(\dot{\mathit{q}})+\mathit{B}\ddot{\mathit{\theta }}+{\mathit{\tau }}_{fm}(\dot{\mathit{\theta }})\\ =\mathit{\tau }+{\mathit{\tau }}_{\text{ext}}\end{array}$

3

Note that a completely rigid joint is always taken into consideration for simplified analysis ( $\mathit{K}\to \infty$ ), then $\mathit{\theta }=\mathit{q}$ , and we obtain the dynamic characteristic of n-DOFs serial robot with rigid links and joints as follows

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

4

where $\mathit{M}(\mathit{q})={\mathit{M}}_l(\mathit{q})+\mathit{B}\in {ℝ}^{n\times n}$ represents the lumped inertia and ${\mathit{\tau }}_f(\dot{\mathit{q}})\in {ℝ}^n$ captures the combined friction effects at the link level as ${\mathit{\tau }}_f(\dot{\mathit{q}})={\mathit{F}}_v\dot{\mathit{q}}+{\mathit{F}}_c\text{sign}(\dot{\mathit{q}})$ with viscous friction coefficients ${\mathit{F}}_v={\mathit{F}}_{vl}+{\mathit{F}}_{vm}$ and Coulomb friction coefficients ${\mathit{F}}_c={\mathit{F}}_{cl}+{\mathit{F}}_{cm}$ . Although equation (4) is the ordinary form of the rigid robot and has been thoroughly researched theoretically and experimentally, we derive this equation from engineering aspect for further implementation on a manipulator.

The symbolic dynamic model equation in equation (4) can be firstly derived using Lagrange equation or Newton–Euler method, which can be decomposed such that the dynamic model can be reformulated as a linear regression form. The dynamics equation can be written in linear expression with respect to its dynamic parameters as

$\mathit{\tau }=\mathit{H}(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})\mathit{\pi }$

5

where vector $\mathit{\pi }$ represents the dynamic parameter set, and the matrix $\mathit{H}(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})$ is always referred to as regression matrix or regressor. For each link, its dynamic parameters consist of six components of the inertia tensor, three components of the first moment, mass, and the viscous and Coulomb friction coefficients. Instead of including every individual parameter, a base parameter set is preferable in model identification since some parameters have no influence on system dynamics and some can be regrouped in linear combinations of its original set. The base parameter set has a lower dimension than the original one. After analyzing the symbolic expression of robot’s dynamics, a well-conditioned trajectory is delicately designed for the actual robot to execute in order to excite robot’s dynamic characteristic. Given an observation of robot joint motion $\mathit{q}$ and generalized joint forces/torques $\mathit{\tau }$ for N samples at sampling time $t_k$ , the dynamic parameter vector $\mathit{\pi }$ may be determined by least square regression as

$\widehat{\mathit{\pi }}=\text{arg}\underset{\mathit{\pi }}{\text{min}}\Vert {\mathit{H}}_s\mathit{\pi }-{\mathit{\tau }}_s\Vert$

6

where ${\mathit{H}}_s=\left[\begin{array}{c}\mathit{H}(\mathit{q}(t_1),\dot{\mathit{q}}(t_1),\ddot{\mathit{q}}(t_1)))\\ ⋮\\ \mathit{H}(\mathit{q}(t_N),\dot{\mathit{q}}(t_N),\ddot{\mathit{q}}(t_N)))\end{array}\right]$ and ${\mathit{\tau }}_s=\left[\begin{array}{c}\mathit{\tau }(t_1)\\ ⋮\\ \mathit{\tau }(t_N)\end{array}\right]$ .

Since the stacked regression matrix ${\mathit{H}}_s$ is typically not full rank, the Moore–Penrose pseudoinverse is used to find the estimation solution as $\widehat{\mathit{\theta }}={\mathit{H}}_s^{+}{\mathit{\tau }}_s$ , where superscript + denotes the pseudoinverse operation.

Residual generation under model uncertainties

In human–robot collaboration scenario where humans are close to the robot to perform cooperative tasks by sharing workspace, safety is the fundamental issue to be considered such that effective safe regulations must be utilized to ensure human safety in case of unexpected collisions. To address the collision detection issues for safe pHRI, several model-based approaches have been proposed based on residual signal using the difference between nominal and actual values such as torque, energy, and generalized momentum. ^2,13,15 These approaches have been devoted much attention to detect the unexpected physical collision between the manipulator and its environment, since only the proprioceptive sensors were used in the detection procedure. It is not only economical to keep the expenses low for external sensing, but also easy to implement by reducing the integration elements, leading to lower maintenance costs and potential failures. For this type of collision detection approach, a straightforward residual method can be computed as the difference between the estimated torque $\widehat{\mathit{\tau }}$ and the actual torque $\mathit{\tau }$

$\mathit{r}=\widehat{\mathit{\tau }}-\mathit{\tau }$

7

where $\widehat{\mathit{\tau }}$ results from the inverse dynamic model (4) with the estimated parameters and measured robot motion information. The estimated torque is stated as

$\widehat{\mathit{\tau }}=\widehat{\mathit{M}}(\widehat{\mathit{q}})\ddot{\widehat{\mathit{q}}}+\widehat{\mathit{C}}(\widehat{\mathit{q}},\dot{\widehat{\mathit{q}}})\dot{\widehat{\mathit{q}}}+\widehat{\mathit{G}}(\widehat{\mathit{q}})+{\widehat{\mathit{\tau }}}_f(\dot{\widehat{\mathit{q}}})$

8

where $(\widehat{\text{⋅}})$ denotes the measured or estimated value. For rigid robot with perfect system modeling, accurate dynamic parameters, and well-measured robot motion variables, the residual in equation (7) equals the external torque ${\mathit{\tau }}_{\text{ext}}$

${\mathit{r}}_{\text{ideal}}={\mathit{\tau }}_{\text{ext}}$

9

However, the estimated torque value (8) contains modeling and measurement errors, and it yields

$\mathit{r}={\mathit{\tau }}_{\text{ext}}+\delta \mathit{\tau }$

10

where $\delta \mathit{\tau }$ denotes the torque errors induced by uncertainties including model uncertainties $\delta {\mathit{\tau }}_{\text{mdl}}$ and measurement errors $\delta {\mathit{\tau }}_{\text{meas}}$ . Generally, model uncertainties $\delta {\mathit{\tau }}_{\text{mdl}}$ contain modeling errors due to incomplete or simplified expression of robot dynamics and dynamic parametric errors accounting for the inaccuracy of dynamic parameters such as link inertia and mass. Unfortunately, robot dynamic parameters are always confidential to be available from the robot manufactures even though they can be computed based on computer-aided design (CAD) model or by a well-designed parameter identification experiment, while such identification procedure is always complicated for implementation and the parametric accuracy remains to be verified. Moreover, the inverse dynamic model equation containing such dynamic parameters is too complex to be wildly used in the industrial applications as the computation complexity of inverse dynamic model challenges the limited robot controller capacity. Hence, a simplified model is preferred in the robot dynamic modeling with reduced expression by omitting the items having the irrelevant influence on the overall dynamic property while keeping fairly satisfied performance. Another source of the torque errors comes from the measurement. As only the position encoder is available in acquiring robot motion data for most industrial manipulators, thus joint velocity and its derivatives are often approximated using numerical differentiation, consequently amplifying the existing noise at position level. However, for a well-parameterized position-controlled robot running a commanded motor trajectory with smooth derivatives of higher order, $\mathit{q}$ and its derivatives can be replaced with ${\mathit{q}}_d$ and its derivatives or numerically differentiated after a well-tuned low-pass filter to remove noise. It is clear that the torque errors $\delta \mathit{\tau }$ specify the dependence on the acceleration (inertia uncertainty), speed (friction, centrifugal, and Coriolis uncertainty), and position (gravity uncertainty).

Proposed collision detection scheme

In order to formulate an effective model-based collision detection scheme without additional sensors, identification of dynamic parameter is the fundamental basis for the sensorless collision detection. A major difficulty of the classic dynamic parameter identification lies in the derivation of the symbolic expression of robot dynamics and regrouping the linear form, which requires much effort and attention even with the assistance of commercial software for symbolic computation. In contrast to traditional method, an alternative to handle model identification is machine learning approach. ^36,37,39 This method offers a framework for automatically learning the complex dynamics without incorporating knowledge about the symbolic expression while still keeping a satisfactory performance for robot control and planning. In literature, ³⁹ a fast modeling and identification of robot dynamics was proposed based on Lasso ⁴⁰ and a simple example of a two-link direct-drive robot was demonstrated to show the effectiveness of this promising approach. Although this method can skip the analytical derivation and identify the candidate terms that act on robot dynamics by a regressor selection technique, the condition of this experiment is limited to a planar manipulator that has no gravity and friction effects, which is obviously not suitable for most robotic apparatus. In this section, an approach based on Lasso is presented to remedy the defect in ³⁹ to realize dynamic modeling and identification for a generic articulated manipulator so as for further collision detection scheme.

Gravity identification

Robot dynamic model identification, which has been intensively researched, is absolutely complex compared to the gravity and friction identification. Here, we will first examine the problem of gravity identification and design a method for accurate identification. The identification of gravity is only position-dependent and the gravity item can be expressed in the following linear parameterized form

$\mathit{G}(\mathit{q})={\mathit{H}}_G{\mathit{\pi }}_{\mathit{G}}$

11

where ${\mathit{H}}_G$ is the regression matrix of gravity after regrouping the gravity symbolic expression with respect to the base gravity parameters set ${\mathit{\pi }}_{\mathit{G}}$ . At the same time, it is not easy to extract the pure contribution of gravity from joint torque. Even if the robot is moving in a quasi-static state, the friction effect still plays an unignored role in robot dynamics. In order to obtain the contribution due to gravity effect, the following analysis is conducted by an in-depth study of robot dynamics. Assume two robot configurations ${\mathit{q}}_1$ , ${\mathit{q}}_2$ satisfy

$\begin{array}{l}{\mathit{q}}_1={\mathit{q}}_2=\mathit{\varphi }\\ {\dot{\mathit{q}}}_1=-{\dot{\mathit{q}}}_2=\dot{\mathit{\varphi }}\\ {\ddot{\mathit{q}}}_1={\ddot{\mathit{q}}}_2=\ddot{\mathit{\varphi }}\end{array}$

12

Then adding the resulting torques according to these two configurations yields

${\mathit{\tau }}_1+{\mathit{\tau }}_2=2[\mathit{M}(\mathit{\varphi })\ddot{\mathit{\varphi }}+\mathit{C}(\mathit{\varphi },\dot{\mathit{\varphi }})\dot{\mathit{\varphi }}+\mathit{G}(\mathit{\varphi })]$

13

If $\ddot{\mathit{\varphi }}=0$ and $|\dot{\mathit{\varphi }}|<\epsilon$ (ε is a small positive value), then the previous equation in equation (13) becomes

${\mathit{\tau }}_1+{\mathit{\tau }}_2\approx 2\mathit{G}(\mathit{\varphi })$

14

Hence, the contrition due to gravity in joint torque can be derived if the above motion assumptions are satisfied by moving the robot to a certain position from the opposite direction at a low velocity with a constant value so as to eliminate other torque effects. To estimate the gravity parameters, we extracted the gravity effect from joint torque by gathering data from 100 different robot configurations and applied least squares method on equation (11).

Friction identification

For electromechanical systems including robot, friction effect plays a significant role in system modeling and control, which is usually difficult to predict especially when robot changes its motion direction. Many researchers have proposed different friction models and identification methods for robot control and simulation. ^{41 –44} In this article, we only employ the basic Coulomb and viscous friction model due to its simple description which captures most characteristic of the friction phenomenon yet. Similar to the extraction of gravity effect from joint torque, the friction effect must be processed in order to identify the friction parameters. Suppose two robot motion states satisfy the condition in equation (12), then the difference of these according torques is

${\mathit{\tau }}_1-{\mathit{\tau }}_2=2{\mathit{\tau }}_f(\dot{\mathit{\varphi }})$

15

It is intuitive to design a sinusoidal motion for joint i in the form of $q_i(t)=\text{sin}(\frac{2\pi }{T}t)$ (T is the cycle time) at joint space to meet the requirements of the above-mentioned motion states if two sampling time are chosen as $t_1=kT/4-\Delta T$ for $q_1$ and $t_2=kT/4+\Delta T$ for $q_2$ with $k=1,3\cdots$ and $0<\Delta T<T/4$ . By applying the designed trajectory for the robot to execute and subtracting the joint torques in the according motion states, the equivalent effect due to friction can be acquired based on the analysis in equation (15). The friction torque is characterized in Coulomb–viscous model and is identified in an axis-by-axis manner. The friction of the i th axis could be parameterized in a linear form

${\tau }_{fi}({\dot{q}}_i)=[\text{sign}({\dot{q}}_i){\dot{q}}_i]\left[\begin{array}{c}{F}_{ci}\\ F_{vi}\end{array}\right]$

where $F_{ci}$ and $F_{vi}$ are the Coulomb and viscous parameters of the i th joint, respectively. Gathering several groups of the joint velocity and friction torque and applying the least squares method, the friction parameters could be experimentally calibrated. We gathered 50 pairs of torques at the corresponding motion states in which the absolute values of each component in equation (4) are equal despite of the sign of the friction items.

Remaining torque identification

Since the gravity and friction items have been identified, thus these two components can be removed from robot joint torque for the remaining complex, nonlinear components model learning. Here, the sum of inertial torque $\mathit{M}(\mathit{q})\ddot{\mathit{q}}$ and centrifugal and Coriolis torque $\mathit{C}(\mathit{q},\dot{\mathit{q}})\dot{\mathit{q}}$ is designated as ${\mathit{\tau }}_{\mathit{r}}$ . Using the linearity of the dynamic model, the remaining torque ${\mathit{\tau }}_{\mathit{r}}$ also keeps this linear property with respect to relevant dynamic parameter vector ${\mathit{\pi }}_{\mathit{r}}\in {ℝ}^{\kappa }(\kappa \le 10n)$ in the following form

${\mathit{\tau }}_{\mathit{r}}={\mathit{H}}_{\mathit{r}}(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}}){\mathit{\pi }}_{\mathit{r}}$

16

It is reasonable to consider the elements in regressor ${\mathit{H}}_{\mathit{r}}(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})$ should be linear combinations of certain candidate terms defined by functions of joint position, velocity, and acceleration multiplied by sine and/or cosine value of $\mathit{q}$ . Indeed, each element in ${\mathit{H}}_{\mathit{r}}$ can be calculated as

$H_r^{ij}=\sum _{k=1}^{\mu \rho }{\delta }_{ijk}{\gamma }_k$

17

where $H_r^{ij}$ is element in the $i\text{th}$ row and $j\text{th}$ column of matrix ${\mathit{H}}_{\mathit{r}}$ . Candidate ${\gamma }_k$ in equation (17) is chosen from $\mathbb{A}\otimes \mathbb{B}$ , where $\otimes$ denotes the Kronecker product, $\mathbb{A}=[{\ddot{q}}_1,{\ddot{q}}_2,\cdots ,{\ddot{q}}_n,{\dot{q}}_1{\dot{q}}_1,{\dot{q}}_1{\dot{q}}_2,\cdots ,{\dot{q}}_{n-1}{\dot{q}}_n,{\dot{q}}_n{\dot{q}}_n{]}^T\in {ℝ}^{\mu }$ and $\mathbb{B}=\mathit{\varphi }(q_1)\otimes \mathit{\varphi }(q_2)\otimes \cdots \otimes \mathit{\varphi }(q_n)\in {ℝ}^{\rho },\mathit{\varphi }(q_i)=[1,\text{sin}{q}_i,\text{cos}{q}_i,\text{sin}(2q_i),\text{cos}(2q_i),\cdots ,\text{sin}(m{q}_i),\text{cos}(m{q}_i{)]}^T$ $(i=1,2,\cdots ,n)$ , m is an positive integer and is set to be 2 (see Appendix 1). The dimension of $\mathbb{A}$ is $\mu =(n^2+3n)/2$ and $\rho =(2m+{1)}^n$ for $\mathbb{B}$ . ${\delta }_{ijk}$ can be set to zero when the corresponding ${\gamma }_k$ does not show in ${\mathit{H}}_{\mathit{r}}^{ij}$ or alternatively an unknown nonzero value decided by robot’s structure. Therefore, according to the expression in equation (17), the remaining torque for joint i can be transformed into

${\tau }_r^i=\sum _{j=1}^{\kappa }{H}_r^{ij}{\pi }_r^j=\sum _{k=1}^{\mu \rho }{\gamma }_k{\stackrel{\sim }{\pi }}_r^{ik}$

18

where ${\stackrel{\sim }{\pi }}_{\mathit{r}}^{ik}={\sum }_{j=1}^{\kappa }{\delta }_{ijk}{\pi }_{\mathit{r}}^j$ . Furthermore, by designating ${\stackrel{\sim }{\pi }}_{\mathit{r}}^{ik}$ as the $k\text{th}$ component of vector ${\stackrel{\sim }{\mathit{\pi }}}_{\mathit{r}}^i$ , which is the converted dynamic parameter set for $i\text{th}$ joint, ${\gamma }_k$ as the $k\text{th}$ component of vector ${\stackrel{\sim }{\mathit{H}}}_{\mathit{r}}$ , the expression in equation (18) can be restated in a compact form as

${\tau }_{\mathit{r}}^i={\stackrel{\sim }{\mathit{H}}}_{\mathit{r}}{\stackrel{\sim }{\mathit{\pi }}}_{\mathit{r}}^i$

19

In the transformed expression (19), the regressor ${\stackrel{\sim }{\mathit{H}}}_{\mathit{r}}$ no longer contains the information specified by system mechanical structure, which otherwise has to be provided in equation (16) with cautious and tedious symbolic computation, making the dynamic parameter identification procedure more convenient. Compared to the original parameter set defined in equation (16), the converted parameter set ${\stackrel{\sim }{\mathit{\pi }}}_{\mathit{r}}$ mapping from ${\mathit{\pi }}_{\mathit{r}}$ is identified separately for each joint and has sparse feature due to the conversion, while the burden of symbolic computation is significantly reduced. It is worth noting that the regressor ${\stackrel{\sim }{\mathit{H}}}_{\mathit{r}}$ in equation (19) is equal for all joints and the size of parameter vector ${\stackrel{\sim }{\mathit{\pi }}}_{\mathit{r}}^i$ is significantly larger than 10n, leading to more computation time for the regression. However, the regression is done offline to select an adequate size of the parameters. Thus, the computation time for regression analysis is not a concern here.

In order to identify the dynamic parameter while reducing the feature dimension to a certain degree by shrinking some coefficient estimates toward zero, a regressor selection method named Lasso is a privilege to employ. ³⁹ Lasso is a regression analysis method that incorporates the least squares method to minimize regression error and an $l_1$ penalty term to yield sparse model by forcing some coefficients to zero. Generally, Lasso has the advantage of effectively producing simpler and more interpretable models that involve only a subset of the predicator by introducing an $l_1$ regulation term to penalize large coefficients.

Observing N groups of data and utilizing the Lasso method, the parameter identification problem for $i\text{th}$ joint is formulated as

$\underset{{\stackrel{\sim }{\mathit{\pi }}}_{\mathit{r}}^i}{\text{min}}{\Vert {\stackrel{\sim }{\mathit{H}}}_{\mathit{r}s}{\stackrel{\sim }{\mathit{\pi }}}_{\mathit{r}}^i-{\mathit{\tau }}_{\mathit{r}s}^i\Vert }_2^2+\lambda {\Vert {\stackrel{\sim }{\mathit{\pi }}}_{\mathit{r}}^i\Vert }_1$

20

where ${\stackrel{\sim }{\mathit{H}}}_{\mathit{r}s}=\left[\begin{array}{c}{\stackrel{\sim }{\mathit{H}}}_{\mathit{r}}(\mathit{q}(t_1),\dot{\mathit{q}}(t_1),\ddot{\mathit{q}}(t_1))\\ ⋮\\ {\stackrel{\sim }{\mathit{H}}}_{\mathit{r}}(\mathit{q}(t_N),\dot{\mathit{q}}(t_N),\ddot{\mathit{q}}(t_N))\end{array}\right]$ , ${\mathit{\tau }}_{\mathit{r}s}^i=\left[\begin{array}{c}{\mathit{\tau }}_{\mathit{r}}^i(t_1)\\ ⋮\\ {\mathit{\tau }}_{\mathit{r}}^i(t_N)\end{array}\right]$ , and $\lambda >0$ is an adjustable tunning parameter.

It is sufficient to predict that the vector ${\stackrel{\sim }{\mathit{\pi }}}_{\mathit{r}}^i$ has many zero elements which have to be forced into zero in the regression in equation (20). To this end, the choice of the tuning parameter λ critically trades off the regression variance and sparsity which indicates how many variables are involved in the regression. For instance, the solution of equation (20) is exactly equal to the least squares regression when $\lambda =0$ , and when λ increases to some extent, all coefficients estimated are substantially forced to approximate zero in order to regulate the penalty term. A simple way to tackle this problem is cross-validation which can properly maintain the balance of estimation variance and variable selection, which has been implemented in MATLAB Lasso function.

For classic robot dynamic parameter identification, it is essentially significant to design a well-conditioned excitation trajectory in order to reduce estimation error. ³³ When it comes to machine learning method, it also requires a sufficiently large amount of data that can contain the necessary information for modeling as much as possible to improve estimation quality. However, designing a well-conditioned trajectory is not the focus of this article and a practical trajectory generation approach inspired by finite Fourier series trajectory ³³ is briefly introduced here.

The basic form of the excitation trajectory for the $i\text{th}$ joint is as follows

$\begin{array}{l}{q}_i(t)=\sum _{k=1}^W\left(\frac{a_i^k}{{\omega }_{b}k}\text{sin}({\omega }_{b}kt)-\frac{b_i^k}{{\omega }_{b}k}\text{cos}({\omega }_{b}kt)\right)+q_{i0}\\ {\dot{q}}_i(t)=\sum _{k=1}^W\left(a_i^k\text{cos}({\omega }_{b}kt)+b_i^k\text{sin}({\omega }_{b}kt)\right)\\ {\ddot{q}}_i(t)=\sum _{k=1}^W\left(-a_i^{k}k\text{sin}({\omega }_{b}kt)+b_i^{k}k\text{cos}({\omega }_{b}kt)\right)\end{array}$

21

where ${\omega }_b$ represents the fundamental frequency of the motion for all joints and W is the order of the harmonic of the sine and cosine functions. $a_i^k$ and $b_i^k$ are the parameters deciding the patterns of the trajectory which are determined through optimization. $q_{i0}$ is an auxiliary parameter which guarantees the robot to recover to its initial configuration at the end of each period. In order to avoid a sudden change at the start and end of each period when calculating a continuous trajectory, boundary constraints are imposed on the start and end motion states of each period. Additionally, the exciting trajectory for each joint is bounded for demanded specifications. Given the above requirements, the problem of the excitation trajectory design for each period is formulated as

$\begin{array}{l}\hspace{1em}\hspace{1em}\underset{a_i^k,b_i^k}{\text{min}}J\\ \text{s.t}{\text{.}}_{\hspace{1em}\hspace{1em}}^{}q(0)=q(t_f)=0\\ \hspace{1em}\hspace{1em}\dot{q}(0)=\dot{q}(t_f)=0\\ \hspace{1em}\hspace{1em}\ddot{q}(0)=\ddot{q}(t_f)=0\\ \hspace{1em}\hspace{1em}{q}_{\text{min}}\le q(t)\le q_{\text{max}}\\ \hspace{1em}\hspace{1em}{\dot{q}}_{\text{min}}\le \dot{q}(t)\le {\dot{q}}_{\text{max}}\\ \hspace{1em}\hspace{1em}{\ddot{q}}_{\text{min}}\le \ddot{q}(t)\le {\ddot{q}}_{\text{max}}\end{array}$

22

where J is an optimization objective to be determined afterwards; $q_{\text{max}},{\dot{q}}_{\text{max}}$ , and ${\ddot{q}}_{\text{max}}$ denote the upper physical limits of joint position, velocity, and acceleration, respectively, with $q_{\text{min}},{\dot{q}}_{\text{min}}$ , and ${\ddot{q}}_{\text{min}}$ as the corresponding lower limits of relevant joint variables. Usually, the optimization objective can be selected as the condition number of the regression matrix in order to plan a well-conditioned trajectory for accurate parameter estimation. Instead of optimizing a complex nonlinear objective function, a practical optimization function which can generate the necessary information required in the regression is adopted as

$J=\sum _{i=1}^n\sum _{k=1}^{\mathit{M}}{\Vert {\zeta }^k{a}_i^k+{\eta }_k{b}_i^k\Vert }_2^2$

23

where both ${\zeta }_k$ and ${\eta }_k$ are randomly generated between $-1$ and 1. This selection of optimization objective can not only simplify the symbolic derivation of the nonlinear function used in the study by Swevers et al., ³³ but also can include a variety of motion patterns for a good data fitting hereafter.

Residual threshold design

As the elaboration in previous section, the residual signal consists of two parts including the external torque effects indicating impacts and the unfavorable disturbances due to parametric uncertainties and model inaccuracies which corrupt the collision detection. It is observed that even under the free motion circumstance, the residual signal has the trajectory-dependent variation in which the peaks due to model uncertainties surprisingly larger than the peak produced by the collision, revealing a static threshold above the maximal residual amplitude is not suitable to obtain a good collision detection performance (see Figure 1). In order to overtake the inappropriate static threshold, the peak caused by collision may have to rise so excessively that collision may be critically harmful to some extent.

Figure 1.

Residual signals for free motion (blue, solid line) and collision (red, dashed line) for joint 1.

A reliable and achievable alternative to distinguish the meaningful peaks caused by collision from false alarm peaks provoked by model uncertainties is to adopt a time-varying dynamic threshold. Instead of the method proposed in literature, ¹⁷ the residual signal is considered as moderate when there is no collision; otherwise, further model refinement and identification has to be implemented to modify the residual. Thus, for a relatively accurate model identification, the residual can be bounded with a dynamic threshold dependent on trajectory. The inner structure of the residual is not modeled in detail but deemed to be a black box. When there is no collision, the dynamic threshold for each joint holds

$|r_i|<r_{\text{det}}^i={\mathit{\Psi }}_i{\mathit{u}}_i+{\upsilon }_i,i=1,\cdots ,n$

24

The threshold $r_{\text{det}}^i$ is built on input ${\mathit{u}}_i$ composed of different functions of robot states. Taking the residual and motion states into consideration, it is intuitive to choose the dependence of the residual on acceleration and speed by setting

${\mathit{u}}_i=[|{\ddot{q}}_i|\hspace{1em}|{\dot{q}}_i|\hspace{1em}{e}^{-{\rho }_i{\dot{q}}_i^2}{]}^{\text{T}}$

25

The last item in ${\mathit{u}}_i$ is to compensate for the uncertainty when robot changes its motion direction, and the corresponding parameter ${\rho }_i$ determines the slope of the peaks which is manually tuned for each link. Moreover, the threshold is characterized by a diagonal matrix ${\mathit{\Psi }}_i=\text{diag}({\psi }_{i1},{\psi }_{i2},{\psi }_{i3})\in {ℝ}^{3\times 3}$ which is subsequently determined by optimization. Additional parameter ${\upsilon }_i$ is added to ensure a certain safety margin, which guarantees the robust detection. In order to ascertain the related unknown parameters to construct a functional threshold, an optimization approach is proposed to compute these parameters. For compactness, we only present how the $i\text{th}$ threshold for joint i is practicably designed by solving the following problem

$\begin{array}{l}\hspace{1em}\hspace{1em}\hspace{1em}\underset{{\mathit{\Psi }}_i,{\upsilon }_i}{\text{min}}\alpha \\ \text{s.t.}\hspace{1em}\hspace{1em}\\ \hspace{1em}\hspace{1em}\hspace{1em}\alpha =r_{\text{det}}^i-|r_i|\\ \hspace{1em}\hspace{1em}\hspace{1em}\alpha \ge {\xi }_i\end{array}$

26

where ${\xi }_i>0$ . The definition in equation (26) can not only keep the threshold from the residual with a secure range larger than ${\xi }_i$ to insure detection robustness but also force the threshold to approximate the residual as close as possible to realize accurate detection of tiny external disturbances. The problem in equation (26) is convex and thus can be solved rapidly with an optimal solution. Here, we use the information of robot motion and residual value acquired from the previous training trajectory. Eventually, the binary detection output CD which denotes the occurrence of an external collision is computed according to the comparison between absolute value of residual signal and the proposed threshold

$\text{CD}=\{\begin{array}{lll}0\hfill & \text{if}\hfill & \forall i:|r_i|<r_{\text{det}}^i\hfill \\ 1\hfill & \text{else}\hfill & \hfill \end{array}$

27

Experiment and validation

Experiment setup

The proposed collision detection scheme has been performed on Efort ER3A industrial robot in our lab (see Figure 2). The 6-DOF manipulator has a spherical wrist, 3 kg of payload capability, and maximal stretch of 0.63 m from the base. The link frames of the Efort ER3A robot are illustrated in Figure 3 and geometry information of the manipulator in modified Denavit and Hartenberg (MDH) representation is given in Table 1. The robot is powered with AC servo motors, which are connected with harmonic gears to increase the effective torque. The robot system uses an industrial computer equipped with 2G RAM and a Hilscher cifX bus communication card, running VxWorks RTOS to realize real-time control of the robot servo drivers. The communication card enables the EtherCAT communication between the servo drivers as the slave devices and the industrial computer as the master device, which allows sending motion or torque commands to control robot and receiving real-time robot’s information such as joint positions and currents at every 4 ms. Since the control architecture is open, the users can program the robot through a teach pendant or use the C++ language with a graphical user interface on a personal computer. Hence, with appropriate programming, the robot control system can implement motion control on both joint and Cartesian space and collect data from both proprioceptive sensors (e.g. motor position encoder and motor current sensor) and exteroceptive sensors (e.g. force/torque sensor, vision sensor). The overall overview of robotic control system is shown in Figure 4.

Figure 2.

Experimental platform of Efort ER3A robot.

Figure 3.

Link frames of Efort ER3A robot.

Table 1.

The MDH parameters of Efort ER3A.

Link i	${\alpha }_{i-1}$ (rad)	$d_i$ (m)	$a_{i-1}$ (rad)	${\theta }_i$ (rad)
1	0	0	0	0
2	$\pi /2$	0	0.050	$\pi /2$
3	0	0	0.270	0
4	$\pi /2$	0.299	0.08	0
5	$-\pi /2$	0	0	0
6	$\pi /2$	0	0	0

MDH: modified Denavit and Hartenberg.

Figure 4.

Control framework of Efort ER3A robot.

The robot is position-controlled with a well-tuned proportional-derivative control law embedded in the servo drivers and only the first three joints are activated for simplicity throughout the article. The input provided to the controller is typically reference position ${\mathit{q}}_d$ at the joint level and the available output from the servo drivers including the joint positions $\mathit{q}$ measured by motor-side encoders and the motor currents I. Meanwhile, the useful joint toque $\mathit{\tau }$ can be obtained using the motor currents with the following transformation

$\mathit{\tau }=\mathit{K}I=\text{diag}(k_1\mathrm{,...},k_n)\mathit{I}$

28

where $k_i$ is the current-to-torque gain of the $i\text{th}$ motor provided by the manufacturer as a constant. The acquired motor torque is, however, too noisy to be directly used in distinguishing the external torque as the random noise signal will annihilate the torque transient caused by collision. On the other hand, the high-frequency noise can be greatly reduced using an appropriate low-pass filter, and a smooth torque is retained while an effective filtering can prevent the effect of hard collisions from overly filtering and avoid corrupting the sensitivity of the collision detection scheme in principle. It is worth noting that the implementation of suitable filtering of the torque requires that the main frequency distribution of external impact torque is within the filter cutoff frequency and the filter has a stable structure and no accumulated filtering errors in digital implementation. In order to solve this problem in filer design, a scheme of distinguishing between intended contact and unexpected collision is introduced in the study by Cho et al. ¹⁸ The low-pass filter $H(z)$ was finally chosen as

$H(z)=\sum _{i=0}^{5}h(i)z^{-i}$

29

with $h_0=h_5=0.1118$ , $h_1=h_4=0.1677$ , and $h_2=h_3=0.2205$ , which represent a fifth-order finite impulse response low-pass using a cutoff frequency of 10 Hz and a Kaiser window. The cutoff frequency was selected to remove the noise signal as much as possible while preserve the meaningful torque data. At the same time, the order of the filter was limited to 5 to reduce the resulting time delay.

Dynamics identification and residual threshold design

In the case of the Efort Ertext 3A robot, gravity effect only presents at the second and third joints and the base gravity parameters can be parameterized with four grouped expressions of mass and center mass. Based on the analysis in the gravity and friction identification, the identification procedure was conducted on the industrial robot sequentially. The corresponding parameters of the first three joints are summarized in Table 2.

Table 2.

Dynamic parameters of gravity and friction effects.

Parameters	Unit	Value
$m_3{l}_{3y}$	kg ⋅ m	−0.6595
$m_3{l}_{3x}$	kg ⋅ m	0.2278
$m_2{l}_{2y}$	kg ⋅ m	0.2139
$m_2{l}_{2x}+m_3{a}_2$	kg ⋅ m	0.8865
${\mathit{F}}_{v1}$	Nm ⋅ s/rad	6.2073
${\mathit{F}}_{c1}$	Nm	3.1127
${\mathit{F}}_{v2}$	Nm ⋅ s/rad	7.0136
${\mathit{F}}_{c2}$	Nm	5.2570
${\mathit{F}}_{v3}$	Nm ⋅ s/rad	4.0979
${\mathit{F}}_{c3}$	Nm	4.5163

In our experiment of training trajectory generation, some basic parameters were designed with $W=5$ , fundamental frequency ${\omega }_b=2\pi /5$ . The physical bounds on joint position, velocity, and acceleration are given in Table 3. According to the training trajectory generation scheme based on the optimization problem in equation (22), we finally produced a trajectory composed of 50 periods of motion data which is continuous at the beginning and end of the individual trajectory for training the remaining torque. Figure 5 shows part of joint position, velocity, and acceleration of the training trajectory. Consequently, the demanded motion was executed by the robot to acquire the relevant joint information including actual joint position and torque. For joint actual position, a low-pass, zero phase filter and a central difference algorithm were used to calculate joint velocity and acceleration. Meanwhile, the joint torque effect due to gravity and friction was removed to obtain the remaining torque ${\mathit{\tau }}_r$ for further data process. Once the income and outcome given in equation (20) are available, the Lasso algorithm automatically computes the parameter estimation with variable selection using the package provided by MATLAB.

Table 3.

Motion constraints for training trajectories.

Joint	1	2	3
Position range (°)	$\pm 90$	$\pm 60$	$\pm 60$
Maximum velocity (°/s)	60	60	60
Maximum acceleration (°/s2)	90	90	90

Figure 5.

Reference trajectories for training. (a) Joint position. (b) Joint velocity. (c) Joint acceleration.

In the context of a three joint articulated manipulator, because the robot’s motion is centrosymmetric with respect to joint 1, $\mathit{\varphi }({\mathit{q}}_1)$ can be excluded from set $\mathbb{B}$ . Hence, $\mathbb{B}$ has a reduced number of candidates to 25, while $\mathbb{A}$ has the dimension of 9; the total number of elements in ${\stackrel{\sim }{\mathit{H}}}_r$ is calculated as 225, leading to the dimension of the converted dynamic parameter set ${\stackrel{\sim }{\mathit{\pi }}}_r^{\mathit{i}}$ shares the same dimension value which is larger than the dimension of ${\mathit{\pi }}_r$ .

The regularization progressively strengthens as the value of λ increases, resulting in fewer nonzero regression coefficients and generally larger cross-validated mean square error suggesting inadequate fit of the responses. However, the parameter λ was automatically selected with the optimal value provided by Lasso function. Practically, for coefficient ${\stackrel{\sim }{\pi }}_r^{ik}$ which is less than 0.2, this coefficient is regarded as having negligible influence on the entire dynamics and always discarded in recovering the remaining torque. For a group of λ value, the cross-validation error for each value of λ is automatically computed. At the same time, the relevant coefficient estimates could also be obtained for corresponding λ value. With the trace plot of coefficients and cross-validated mean square error plot of Lasso fit, one could visually observe the difference of λ selection making on the coefficient dimension reduction and estimation error. If the zero components in ${\stackrel{\sim }{\mathit{\pi }}}_r^i$ are selected more, that is, higher value of λ parameter in the Lasso algorithm is preferred, then the threshold margin ${\xi }_i$ may be manually set more conservative. In our implementation, we chose the λ which reduces the dimension of the coefficients larger than 0.2 to 10 and keeps the estimation error satisfactory as well. Eventually, the recovered robot inverse dynamics composed of the abovementioned three parts and the actual toque is illustrated in Figure 6 corresponding to the motion in Figure 5.

Figure 6.

Actual joint torque and estimated torque based on the proposed method for trajectories in Figure 5.

Figure 6 compares the measured and predicted torques, revealing the inverse dynamics can be reconstructed given the demanded motion. It is worth noting that the predicted torque roughly fits the actual torque well expect some derivation and there exists slight torque variation during the motion because a simplified dynamics identification process is adopted in this article. Meanwhile, the absolute residual value under free motion in Figure 5 can be calculated once the inverse dynamics is available, and the torque residual for each joint is depicted in Figure 7.

Figure 7.

Residual signals for free motion (blue, solid line), collision (red, dashed line), and threshold (green, dashed line) for training trajectory in Figure 5.

It is observed that the residual signal fluctuates near zero except an abrupt peak arises when velocity changes its sign, which implies the reliability and practicality of the proposed dynamics identification. The abrupt peak arises when velocity changes its direction since it is difficult to accurately model the friction in low-velocity range. When velocity traverses the low-velocity regime, several phenomena like Stribeck effect, predisplacement, rate dependence, and hysteresis will significantly dominate the dynamics of friction. To address these dynamic behaviors in low-velocity regime, several friction models incorporating more coefficients and complicated structure were proposed such as LuGre model, Maxwell-slip model, and Leuven model. These models could to some extent illustrate the property of friction dynamics in low velocity, but the identification procedure is not trivial and the accuracy cannot be guaranteed. Indeed, the Coulomb–viscous model is still the most popular one in constructing the friction in robotics due to its linearity in parameters and satisfactory description of friction dynamics in most velocity regime.

In addition, collisions were imposed by operator on the manipulator when it executed the same trajectory and corresponding residual signals were computed as presented in Figure 7. In contrast to collision-free scenario, the peaks due to collisions have the influence on the residual as the model uncertainties do. Hence, the torque errors induced by the dynamics uncertainties including unmodeled dynamics and estimation bias will cause adverse effects on accurate collision detection, and this is tackled by employing the dynamic threshold.

The parameters for threshold is given in Table 4, and the corresponding reconstructed threshold is compared with the torque residual, as shown in Figure 7, illustrating the acquired threshold can bound the residual signal to ensure the detection mechanism not be activated by false alarms due to measurement errors and dynamic model uncertainties as well as small enough to allow fast detection of collision. The effectiveness of the proposed threshold was validated for a different trajectory as detailed in the following section.

Table 4.

Parameters of the proposed dynamic threshold.

Joint i	1	2	3
${\psi }_{i1}$	−0.1433	−0.335	−0.7733
${\psi }_{i2}$	0.2295	0.7350	1.3615
${\psi }_{i3}$	5.1588	5.5866	9.1100
${\upsilon }_i$	3.4021	5.6258	3.6457
${\xi }_i$	2.0000	2.0000	2.0000
${\rho }_i$	150	200	150

Validation

The proposed scheme for sensorless robot collision detection was experimentally validated on the industrial robot ER3A performing a continuous cyclic motion task with a triangle curve in Cartesian space. The triangle motion in Cartesian space is with the maximal velocity of 0.2 m/s and trajectories in the joint space are calculated using the inverse kinematics.

Figure 8 shows the residual signal for each joint under different collision conditions. For example, a collision near 10 and 14 s is relatively moderate and these joint peaks may not be detected if static threshold was adopted. Thanks to the proposed dynamic threshold, the tiny torque peaks resulting from collision could be sensitively identified. For example, a collision occurred near 14 s and the effective torque caused by collision at joint 1 is smaller than the maximal dynamics error induced by model uncertainties. The dynamic residual threshold can also effectively recognize the source of the torque peaks. For the collision at t = 18 s, the peaks in residual signals were recognized at all joints, suggesting the unexpected collision with slight impact occurred at the body before the third joint. The effective external torque at joint level satisfies $\mathit{\tau }={\mathit{J}}_c^T{\mathit{F}}_{\text{ext}}$ , where ${\mathit{F}}_{\text{ext}}$ is the Cartesian contact force and ${\mathit{J}}_c(\mathit{q})$ is the Jacobian matrix at the collision point. When the external force ${\mathit{F}}_{\text{ext}}$ belongs to the kernel of ${\mathit{J}}_c^T(\mathit{q})$ , the external force will be badly reflected at the joint level. The collision peak at t = 14 s is only detected at the first joint, while the slight peak at the second joint could not be sensed. This may be owing to that the matrix ${\mathit{J}}_c^T(\mathit{q})$ is ill-conditioned at that configuration. Generally speaking, the collision detection scheme could locate the probable collision position, which may be utilized for further collision reaction scheme. In our article, a simple collision reaction was adopted by stopping the manipulator for 1 s to prevent harm to operators as soon as the collision was detected.

Figure 8.

Residual signals for free motion (blue, solid line), collision (red, dashed line), and threshold (green, dashed) for validation trajectory.

Conclusion

In this article, a current-based whole body collision detection is presented on a conventional industrial manipulator based on a quick dynamic model identification using Lasso and a residual evaluation method to overcome the dynamic model uncertainties. The presented method uses only proprioceptive sensors providing motor-side information without additional exteroceptive sensors. The collision effect is diagnosed as system fault through comparison between the actual torques related to the internal motor currents and the reconstructed torque from the joint position and its derivatives. In this framework, a sequential identification of dynamic parameters is introduced where the major components of the dynamics (i.e. gravity and friction) in pHRI scenario are identified in the first step and the remaining torque is formulated with a simplified model. On the basis of this model, Lasso algorithm is adopted to predict the relevant parameters along with parameter selection in order to capture the characteristic of the remaining torque. The residual signal which is defined as the difference of the actual and predicted torque is used to recognize whether the external impact is imposed on the robot. To tackle the irregular peaks due to dynamic uncertainties, dynamic threshold is imposed on the residual which ensures the sensitivity of the detection scheme as well as the robustness to false alarm caused by system errors.

The method has been implemented on a conventional industrial robot ER3A with EtherCAT protocol for real-time data communication. Instead of the classic parameter identification with complicated symbolic derivation and heavy computation burden for real-time control, the proposed method can fast and effectively reconstruct the inverse dynamic model. Consequently, the residual signal for collision detection is generated through the simplified model, and threshold is conceived to remedy the peaks due to system uncertainties. Compared with a static threshold, the dynamic threshold reduces the collision torque sensitivity from 10 Nm to 5 Nm while maintaining the robustness. Although the joint sensitivity does not share linear relationship with respect to the impact sensitivity in the Cartesian space, the impact can even be impossibly sensed if the impact force belongs to transposed Jacobian matrix kernel space. However, the robot always holds favorable configurations when robot shares workspace when collaborating with humans; external impact is reflected on the joint space with high fidelity which guarantees the detection of collisions occurs on the whole body of the robot. In our experiment, the collisions can be detected with the sensitivity less than 5 Nm in joint space with a short time delay less than 50 ms.

Despite a fast and simplified model identification, a dynamic threshold for residual signal is proposed to achieve a satisfactory whole body collision detection on an industrial manipulator; our future plan is to evaluate the sensitivity of forces in Cartesian space of the collision classification using a 6D F/T sensor. Additionally, the presented idea in recovering dynamics and residual threshold generation can be brought into the external force estimation at the end effector so as to realize a low-cost virtual force sensor which benefits the force-control scenario like assembly and human–robot cooperation. Further investigation directions in the collision detection also include collision detection scheme of a robot with payload fixed at the end effector and how the factors like speed and stiffness affect the accuracy of the collision and corresponding strategies.