2.1 Adaptive haptic assistance in telerobotics

Adaptive haptic assistance in telerobotics still is a rather young field of research. A study with different haptic assistance policies in telerobotics is presented in [⁷]. In their teleoperation scenarios, a virtual box has to be guided through a maze without and with obstacles that have to be avoided, and GVFs are implemented to support the user. In the study, one constant assistance policy and two adapting policies (switching and linear continuously adapting) are compared according to task performance, disagreement, perceived workload by the users and efficiency. One conclusion of that study is that adaptive assistance policies have the potential to outperform constant assistance when taking into account multiple performance measures. For known environments with obstacles in particular, adapting policies provide more flexibility. However, the authors also conclude that it cannot be said generally that continuously adapting policies are better than switching policies or vice versa according to their evaluation metrics.

The work in [⁸] also proposes a study to compare constant and adaptive haptic assistance in telerobotics.

The focus of that study is the analysis of how the assistance system behaves when unpredictable disturbances arise. Therefore, scenarios with disturbances forces are established; also, adapting the assistance is developed by changing the stiffness of the GVF linearly with the user grip force measured by force sensors. The authors conclude that adaptive assistance achieves similar performance to high gain (strong guidance) constant assistance. However, adaptive assistance reduces the steering force necessary to overcome situations with unpredictable disturbances.

Corredor et al. [⁶] developed a nonlinear continuously adapting assistance policy with the intention to make the assistance react more “human-like” in partially known environments. The assistance system is based on the drift-diffusion model (DDM), a decision-making model from cognitive science. The DDM represents a decision-making process as a Two-Alternative Force Choice (TAFC) task depending on past decisions. In [⁶], the DDM is integrated into the assistance system to adapt the impedance of a spring-damper GVF in order to improve the performance of the teleoperation task with unknown obstacles. The adaptive assistance proposed by Corredor et al. is more complex than systems with linearly adapting policies, which makes it interesting for further investigation. That study concludes that adaptive strategies that preserve the transparency of the interaction present better performance.

2.2 Sensorless human force estimation for HRI

In the field of sensorless external force estimation in robotics, a lot of research has been done, where human force estimation in HRI can be seen as a special case of sensorless external force estimation. In [⁹], the authors propose a neural network-based contact force observer for haptic applications. They use an MLP neural network to model the inverse dynamics of a 3-DOF Planar Twin-Pantograph haptic device. However, only sensor signals from the joint angles are available, so the velocities and accelerations are obtained using derivative filters. In [¹⁰], the problem of estimating human force, while lifting heavy unknown masses with robotic assistance, is tackled. The human force estimation is formulated as a linear control problem where precise knowledge about the robot dynamics is necessary. The work in [¹¹] proposes a twin robot system called shadow robot to estimate human force. The system consists of two robots of the same type, where one robot is constrained by the human, while the other one can move freely. The robots are controlled in position, velocity and acceleration via a control based on a disturbance observer. By subtracting the disturbance torque in the unconstrained robot from the constrained one, the human force is extracted. In [¹²], the nonlinear disturbance observer proposed by [¹³] is used for human force estimation in a 1 DOF haptic device. To estimate the human force, a parametric dynamic model (including friction terms) of the device is required. This allows the separation of the observed disturbance torque into a friction term and the applied human force. In [¹⁴], a Kalman filter disturbance observer was implemented in order to estimate human force for a delta robot. The disturbance torque is estimated by using measures of the motor currents and the joint angles. That disturbance torque observer is very similar to the one used in [¹¹]. However, this work does not mention how to deal with model uncertainties and unmodeled dynamics in order to separate the torque caused by the operator from those due to disturbances. A sensorless approach to estimate human force for collision detection in HRI is proposed in [¹⁵]. A collision torque observer, based on a nominal dynamic model with generalized momentum, is used. The use of the generalized momentum makes acceleration information unnecessary, which is usually not directly measured in robotic systems. Hence, only the joint angles, angular velocities and motor torques are needed for the proposed approach. Model uncertainties and friction parameters are identified using a genetic algorithm.

In [¹⁶-¹⁸], model-based techniques are proposed in order to estimate external forces under communication time delays. Convergence proofs rely on Lyapunov methods and these works concentrate mostly on enhancing robustness against time-delays of the estimator.

The approach we follow is somewhat different in the sense that we opt for a data-driven approach, where the model is trained a priori. In this sense, time delays are assumed to be negligible. Although this assumption is strong, we prefer to concentrate on the mechanics of shared control and the effect that the lack of force sensors has on the overall task execution performance.

3.1 Nominal master robot dynamics Kinematics

The vector of joint variables is defined as , where n _q is the number of the actuated joints and the robot’s position with respect to a fixed workspace coordinate frame is described by the vector , where np is the number of DOF of the end-effector. The Novint Falcon has three DOF and three actuated joints, thus, n _q = n _p = 3 with parameters as presented in [¹⁹].

In [²⁰], the inverse kinematics in the coordinate frames L _i (see Figures 2 and 3 for the position of the coordinate frames) are derived. Because the end-effector’s position is only known with respect to H, the coordinates are transformed to L _i . The position of the coordinate frame H

with respect to L _i can be described by the homogeneous transformation matrix (see Equation (1)]

with the rotation matrix and the displacement vector The values for the rotation angles and displacements are given in Table 1, [¹⁹]. Hence, the coordinates of the end-effector’s position p with respect to L _i can be transformed using (2). From (3), (4) and (5) the equations for the inverse kinematics can be obtained as

from the known end-effector position ^H p. Using trigonometric relations, elements of ^Li p can be expressed as functions of the joint angles θ _ji as

The detailed steps for the derivation of the inverse kinematics can be found in [²⁰]. Due to the constructive design of the Novint Falcon and the restrictions of the joint angles, a unique solution, out of the four theoretical possible solutions, can be obtained (see (6) - (8)], as:

Equation (7) may be used; as an alternative to (8), the following equation may be used instead:

The Jacobian matrix J for the Novint Falcon is derived. With J, velocities and forces can be transformed from the joint space to the workspace and vice versa [²¹] using the relation in (9).

with f being the force vector at the end-effector and τ being the torques at the actuated joints.

A derivation of J for the Novint Falcon can be found in [²²]. The inverse Jacobian matrix is given by:

3.2 Joint velocity estimation

For the application of the NDOB, information about the joint velocities is necessary. However, those velocities cannot be measured directly and thus, they have to be estimated from the computed joint angles. A simple “dirty” differentiator is used (see Equation (13)], i.e., velocities are obtained by applying the low pass filter

Hence, the estimated joint velocities are declared as filtered velocities The cut-off frequency is chosen as ω _c = 100 rad/s. It is assumed that the Novint Falcon end effector does not move with a higher frequency because it is constrained by the human operator during task execution and the maximum tremor frequency of the human hand can be located at approximately 12 Hz [²⁵]. The advantage of the joint velocity estimation with (13) is its computational simplicity, which leads to a low computational effort and thus, it is useful for real-time applications at the expense of being less precise.

3.3 Nonlinear disturbance observer

In this section, the NDOB for robots proposed in [¹³,²⁶] is revisited. Define a new variable that sums the gravity, Coriolis and centrifugal terms as in Equation (14)

The nominal model is not exact, hence the robot dynamics (see Equation (15)] consist of the nominal model and model uncertainties:

Furthermore, τ _c are the actuation torques (control torques). Rigid body movements of a robot are usually described by the ordinary differential Equation (10). Besides rigid body movements, no exogenous excitations such as friction are considered in the model. Inserting (14) into (10), inserting (15) into the resulting equation and then adding a friction τ f and external forces τ _ext term, leads to the extended model:

The uncertainties of the extended models are lumped in a disturbance torque (16)

Thus, the extended model has the form

Chen et al. propose in [¹³] a NDOB for an estimation of the disturbance torque Based on the observer (17)

Where is the observer gain matrix. Alternatively, the auxiliary variable , where is a user-defined function and may be used for estimation.

Taking the time derivative of z as in Equation (18)

and defining (19)

the time derivative (18) may be rewritten as

In contrast to the observer (17), the new observer (20) does not need information about the joint acceleration (see [¹³]). Mohammadi et al. [²⁶] generalize the observer (20) and the way to obtain an observer gain matrix. They define

with X being a constant, invertible n × n matrix and M _no being symmetric and positive definite. By inserting (21) into (19), it can be shown that The convergence may be established (see [²⁶]) via Lyapunov analysis. The disturbance tracking error converges exponentially to zero if there exists positive definite and symmetric matrix X and positive semi-definite matrix Γ such that Equation (22)

Furthermore, Mohammadi et al. show that the minimum convergence rate of the disturbance tracking error is , where denotes the minimum eigenvalue of a matrix with parameters as in Equation (23)

Using the upper bound on , inequality (22) can then be rewritten as the linear matrix inequality (24) (with a specific choice of parameters (25)]

with Y = X ⁻¹ . An explicit analytical solution of the LMI (24) is obtained for the case

with y and κ being constant scalars. It can be shown that the case (25) leads to with β being the minimum convergence rate of the tracking error. As a trade-off between fast convergence and noise sensitivity, a suitable choice, as in Equation (26) is proposed as

as an analytical solution of the LMI.

Regarding (23), the upper bounds of are necessary. However, the 2-norm of those matrices is hardly feasible in analytical form. Using (11) and the triangle inequality,

where the identity ∥J^T∥₂ = ∥J∥₂ has been used. With (27), a conservative estimation for σ can be obtained. However, the problem of how to obtain max{∥J∥₂} still remains. Due to this, it has been decided to approximate max{∥J∥₂} by computing ∥J∥₂ for a set of 106 uniformly distributed end effector positions. Thus, σ ≤ 0.1116 is computed from (27).

The time derivative of the inertia matrix

can be obtained from the skew-symmetric property of the matrix M˙ (q) − 2C (q, q˙) [²³] and by taking (11) for the Coriolis and centrifugal matrix of the nominal model, the 2-norm can be upper bounded by

The bound is obtained by using the identity . Obtaining an approximation of is not as trivial as with because the dimension of the search space for the upper bound increases from three to six dimensions due to the additional independent joint velocities . Hence, the assumption is made which leads to ζ ≤ 0.6024.

Besides information about the upper bounds of M_no (q) and , also a minimum convergence rate must be determined for (26). Wang et al. [²⁷] determined the minimum convergence rate of an observer for an automotive steering wheel according to the maximum input frequency f _max of the system. They analyzed the step response of the observer and chose a convergence rate, so that within a time of the threshold value is reached close enough. Taking ω _c from (13) as the maximum angular frequency input ω _max . Now β is chosen high enough, so that the observer reaches 95 % of the threshold value within seconds. The value for β is found by trial and error via simulation of [²⁰]. With β ≥ 2.8, a sufficiently high convergence rate is found. Finally, with the values for σ, ζ and β the gain matrix Y = 0.6137 I is obtained.

3.4 Neural network-based estimation of model uncertainties and unmodeled dynamics

According to ⁽¹⁶⁾, the model uncertainties and unmodeled dynamics of the master robot are

The idea is now to estimate τ u with a NN in order to separate τ ext from the τ d. For the generation of training and test data, the Novint Falcon end-effector tracks a reference position trajectory under closed-loop control, where the end-effector positions and control forces are recorded and transformed into joint space using the already proposed kinematic relations. The joint velocities are then estimated from the joint angles using the filter [¹³]. It should be mentioned that during the data generation, no human interaction with the robot takes place; thus, τ _ext = 0.

For the training of the NN, two different methods are proposed and compared. The first method is inspired by the work of Hu and Xiong [²⁸]. In this approach, a semiparametric model is proposed to identify the robot dynamics. The training scheme for Method 1 is shown in Figure 4. The nominal model is simulated with the recorded motion data as input. Additionally to the joint angles and velocities, joint accelerations are also necessary, which are estimated from using [¹³]. The difference between τ _c and τ _no is τ _u, is the training target for a multilayer perceptron (MLP) NN.

The second method assumes that no human interaction with the master robot leads to τ _d = τ _u. That means that the NDOB [²⁰] can be used to produce the target data for the NN. The training scheme for Method 2 is shown in Figure 5.

For both methods, the NN has six inputs: three joint angles θ _1i and their respective joint velocities . The output is the three-dimensional vector The NN itself contains, in each case, one hidden layer with five neurons. Although that configuration is rather simple, no improvements in the estimation with more complex NNs (more hidden layers and/or neurons) were observed. As activation functions, hyperbolic tangent sigmoid functions are used, while the neurons in the output layer have linear activation functions. A set of 10,000 samples is used for training and the NNs are trained with the Levenberg-Marquardt backpropagation algorithm, while 10 % of the training data is used for training validation.

The resulting NN’s do not approximate fast varying disturbances and peaks very well. Peaks in the torques arise on the one hand, due to peaks in joint velocities, on the other hand due to inauspicious combinations of joint positions and velocities in the computation of Due to noticeably lower RMSE and MAE, the NN trained by Method 2 is used in the human force estimation in this work.

4.1 Master robot

The Novint Falcon is a low-cost haptic control device for video games, developed by the former company Novint Technologies, Inc. It is capable to produce forces of up to 2 lbs and its workspace is roughly 4 × 4 × 4 inches. The control torques are updated with a frequency of 1 kHz, which is a typical value for a realistic haptic sensation [²⁹]. In the teleoperation setup, the Novint Falcon communicates with MATLAB/Simulink via the Haptik Library [³⁰].

Novint Falcon communicates with MATLAB/Simulink via the Haptik Library [³⁰].

4.2 Slave robot control

The virtual mobile robot is simplified to a holonomic point moving in a 2D workspace. The human operator controls the velocity of the slave robot by moving the Novint Falcon end-effector. There, the x-y-elements of Hp are proportional to the velocity vector of the slave robot. By integrating the slave robot’s velocity with respect to time, the actual position of the mobile robot is obtained.

4.3 Adaptive assistance

In this subsection, the adaptive haptic assistance proposed in [⁶] is revisited. The idea of the adaptive assistance is to schedule the level of assistance depending on the task objective, which can change over time due to uncertain environments. A decision-making model decides between two different objectives: to improve task performance or interaction performance. Depending on the decisions, the properties of the GVF are changed. The GVF in this work is designed as a system that forces the robot to return to the desired path with a control action given by

with fa⊥ being the assistance force perpendicular to the desired path, d⊥ the damping coefficient, k⊥ the stiffness and the deviation of the slave robot from the desired path defined as

The stiffness of the GVF, which is the only parameter varied by the adaptive assistance here, is mapped as a linear homotopy function k_⊥ = αk _low + (1 − α)k _high , with assistance level α ∈ [0; 1[ and upper and the lower bounds for the stiffness k _high and k _low . It can be seen that a low α leads to a more stiff GVF and a high α to a less stiff GVF.

The performance of the teleoperation task is evaluated by two different metrics: Task performance TP and interaction performance IP . TP is based on the position error ϵ and is normalized by

In this work, the position error is only the deviation from the desired path in y-direction and thus, ϵ is a scalar. For mean (|ϵ|), the RMS value is given by (28)

IP is calculated using the internal force fi between the human operator and the assistance. According to [³¹], the internal force is defined in Equations (29) and (30)

with f_⊥ being the human force and f_a⊥ the assistance force. For the performance metric, f _i is normalized by

The forces perpendicular to the desired path are assumed one-dimensional, hence mean (|f _i |) may be defined as the RMS value. A high internal force indicates a low agreement with the assistance; hence the interaction performance is defined as I_P,n = 1 − f _i,n. The idea of the DDM is that every choice z of the TAFC task is rewarded according to the function,

The rewards depend on the performance measures and are computed by reward functions. As a reward structure, amatching shoulder with the linear reward functions

is used. There, k _I and k _T are the slopes of the reward functions, I _P0 is the minimum reward for interaction performance and T _P0 is the maximum reward for task performance. Reward functions can be designed with or without any preference for task or interaction performance, depending on the teleoperation task. For a neutral reward structure, the intersection point of both reward functions is in the middle of the performance measures, which is 0.5 for normalized performance measures. By moving the intersection point to the right or left half, either preference on interaction or task performance is obtained. An elementary part of the DDM is taking into account the actual and past rewards in order to accumulate evidence for future decisions. The evidence in favor of one decision is computed as

with w _z being the evidence for the actual choice being the evidence for the option which was not chosen z̅, and λ being the learning rate with λ ∈ [0; 1[. If λ is low, past decisions are taken more into account, and vice-versa. With the accumulated evidence, the probability of decision making in order to improve interaction performance p _I is computed using the soft-max model

For a large value for µ, the user might feel sudden changes in the stiffness of the GVF. If pI exceeds a certain value (0.5 in this work), future decisions are made in favor of improving interaction performance (see Equation (31)],

This means that more authority is given to the human operator by increasing α. It is now possible to define the adaptation rule as α(pI ) = pI . Tests with adaptive assistance have shown that some modifications in the task performance measure and the reward structure are necessary. It has been observed that T _P,n increases rapidly for very small path errors, therefore T _P,n is weighted with

In (32), ρ is a design variable that must be determined individually for a teleoperation task. For large ρ, small path errors have less impact on the task performance measure.

Furthermore, it has been observed that reward functions with an intersection point at T _P,n = I _P,n = 0.5 do not practically lead to a neutral reward structure in this teleoperation setup, but lead to a reward structure with preference on task performance. The parameters for that neutral reward structure are defined as k _I = 1, I _P0 = 0, k _T = −0.5, T _P0 = 0.5

Figure 8 shows the DDM response using the parameters of neutral reward structure for a path following task with an obstacle avoidance event. It can be observed that the stiffness is high when following the path and decreases when moving around the obstacle.

Besides the adapting assistance force, which acts in the perpendicular direction of the desired path (y-direction), a virtual fixture for the x and z-direction of the Novint Falcon is also considered. The virtual fixture for the x-direction has a spring model f _a,x = −k _x ^H p _x , _master . For the z-direction, a virtual fixture with the damper model f _a,z = −d _z ^H p˙ _z,master is used. Finally, the assistance force vector f _a, which is the force vector sent to the master robot, consists of f _a = [f _a,x f _a⊥ f _a,z ] ^T .

5.1 Experiment

The teleoperation task is to follow a desired horizontal path from left to right and to avoid the unknown obstacle during task execution (see Figure 7]. The experimental setup does not include communications or any other sources of time delays. As soon as the x-coordinate of the desired end position is reached, the task ends. Five subjects (one female and four male) have participated in the experiment as human operators. Four of them did not have any experience in telerobotics, while one person already had some experience. Each person was allowed to do a trial run without obstacle or assistance in order to get familiar with the control of the virtual slave robot by the Novint Falcon before the experiment started. The users were told to complete the task within 15 seconds, but the task completion time was not considered in the evaluation of the experiment.

The user study is executed in two sessions. In the first session, the users perform the task with four different constant assistance levels (one assistance level per run). These assistance levels are classified as no, low, medium and high assistance. For every assistance level, two runs are performed, but in a random order. As a result, eight runs in total are performed for the first session. The second session is executed with adaptive assistance. There are three different reward structures (no preference, preference on task performance and preference on interaction performance) and two different values for the learning rate λ (high and low). This leads to six different configurations for the adaptive assistance. In addition, in the second session, the users perform the teleoperation task twice for each configuration (in random order). The parameters for each assistance configuration are given in Table 3.

A decision in the DDM is made every 100 ms; a measurement of the path tracking error, as well as an estimation of the human force, is made every 10 ms. The experiment has been executed with MATLAB R2015b (32-bit) on a computer running with Windows 7 Ultimate. The computer is equipped with 12 GB RAM, an Intel Core i5-3330 (4×3 GHz) and an Intel HD Graphics graphic board. The teleoperation task was visualized on a 19-inch screen for the human operator.

5.2 Evaluation metrics

Due to the presence of obstacles that must be avoided during the task execution, the workspace of the slave robot is divided into different evaluation regions. While moving on the path, the mobile robot stays in the on-path region. That obstacle-free region is divided into two segments: one in front of and one behind the obstacle. The area between these two on-path regions is the off-path region, where the obstacle is located. In order to prevent learning effects during the experiment, the obstacle size and position are randomly chosen for every run and thus, the threshold value is determined as 25 % of the obstacle edge length. The obstacle size varies between 20 and 40 length units.

5.3 Results

A Kolmogorov-Smirnov test with a significance level of 5 % is applied to the data sets in order to check if the data has a normal distribution. Next, an ANOVA test is made in order to check if the data of at least two assistance configurations differ significantly from each other in a statistical sense.

If the p-value from the ANOVA is lower than 0.05, the null hypothesis (i.e., all assistance configurations are the same) is rejected. If this is the case, Tukey’s HSD test is applied in order to find the assistance configurations which differ from each other.

5.4 Discussion

In this subsection, the results of the user study are discussed referring to the questions which have been formulated in the introduction of this section.