3.1. Decision layer

There are two command modes generated by the decision layer for the control of the mobile robot: force and position commands modes. Force command mode: This mode corresponds to a completely restricted mode, that is, the human brain converts the task into a certain desired force output according to the situation of the object in the environment. The desired position is set to be equal to the actual position of the human hand, that is,

(1) \begin{align}{X_d} = {X_h}\end{align}

where ${X_d}$ and ${X_h}$ are the desired and current human hand’s positions, respectively.

Position command mode: This mode corresponds to the free motion mode, that is, the human brain converts the task into a certain desired position output according to the situation of the object in the environment, and the expected force is equal to the actual force of the human hand, that is,

(2) \begin{align}{F_d} = {F_h}\end{align}

where ${F_d}$ and ${F_h}$ are the desired and actual human hand’s forces, respectively.

3.2. Fusion layer

In the fusion layer, the operator’s positioning process and force determination process are fused. The output force at the fusion layer, ${F_s}$ c, can be obtained in s domain as follows:

(3) \begin{align}{F_s}(s) = \frac{{{T_p}s + 1}}{{{T_{pl}}s + 1}}{k_p}{e_p}(s) + \frac{{{T_f}s + 1}}{{{T_{fl}}s + 1}}{k_f}{e_f}(s)\end{align}

where ${e_p}(s) = {X_d}(s) - {X_h}(s) + {r_p}(n)$ , ${e_f}(s) = {F_d}(s) - {F_h}(s) + {r_f}(n)$ , ${r_p}(n)$ and ${r_f}(n)$ are the white noises generated from the contacted object’s position and force errors, respectively, which are very small to be eliminated without any problem [Reference Talukdar and Barman24].

3.3. Conductive layer

In the conductive layer, a value obtained by applying the delay time of nerve transmission is transmitted to the movement of the human arm. The output force at the conductive layer, $F$ c, can be described as

(4) \begin{align}F(s) = {F_s}(s){e^{ - \sigma s}}\end{align}

where $\sigma $ [second] represents the time delay of the transmission nerves.

3.4. Human arm’s motion

The human arm’s motion after the gravity compensation in Fig. 2 can be described using the impedance modeling as follows:

(5) \begin{align}F(s) + r(n) & = ({M_h}{s^2} + {B_h}s + {K_h}) \cdot {X_h}(s)\nonumber\\[4pt]& = {Z_h} \cdot {X_h}(s), \end{align}

where $r(n)$ is the white noise caused by muscle fatigue, which can be negligible; ${M_h},{B_h},$ and ${K_h}$ care the mass, damping, and stiffness coefficients of the human arm, respectively.

From Fig. 3, $F(s)$ can be obtained after the blocks of ${G_p}(s)$ and ${G_f}(s)$ . As a result Equation (5) can be represented as follows:

(6) \begin{align}{G_p}(s)({X_d}(s) - {X_h}(s)) + {G_f}(s)({F_d}(s) - {F_h}(s)) = {Z_h} \cdot {X_h}(s)\end{align}

where ${G_p}(s)$ and ${G_f}(s)$ are obtained from Fig. 3 as follows:

(7) \begin{align}{G_p}(s) = \frac{{{T_p}s + 1}}{{{T_{pl}}s + 1}}{k_p}{e^{ - \sigma s}}\end{align}

(8) \begin{align}{G_f}(s) = \frac{{{T_f}s + 1}}{{{T_{fl}}s + 1}}{k_f}{e^{ - \sigma s}}\end{align}

${G_P}(s)$ is a position factor and ${G_f}(s)$ is a force factor. When the operator works in position command mode, the position prediction factor becomes the dominant factor in determining human motion and other factors can be eliminated. That is, ${G_P}(s)$ and ${G_f}(s)$ can be simplified ignoring the delay time and the lag factors, ${T_{pl}}$ and ${T_{fl}}$ , which can be approximated to be zero. Also the force prediction factor, ${T_f}$ , can be eliminated. Equations (7) and (8) can be simplified as follows:

(9) \begin{align}{G_p}(s) = ({T_p}s + 1){k_p}\end{align}

(10) \begin{align}{G_f}(s) = {k_f}\end{align}

As the result, the Eq. (6) can be presented as follows:

(11) \begin{align} \left(T_p s + 1\right)k_p(X_d(s)-X_h(s)) + (F_d(s)-F_h(s)) = Z_h X_h.\end{align}

According to the neurophysiology research results [Reference Merfeld and Clark25], the response time distribution of the operator in the processes of perception, decision, and control is as follows: the process by which the receptor transforms the stimulus into a nerve impulse needs 1 to 38 ms, and the nerve-conduction process needs 2 to 100 ms. The neural decision-making process requires 70 to 300 ms, and the process by which the efferent nerve transmits the command to the neuro muscle requires 10 to 20 ms.

As the operating frequency of the human robot interaction system is affected by the operating bandwidth of the actuator, it generally works in the low frequency region. Some factors can be ignored to simplify the operator arm’s dynamic modeling and the modeling can be represented as follows:

(12) \begin{align}{F_h}(s) & = \left({M_h}{s^2} + {B_h}s + {K_h}\right) \cdot \Delta X(s)\nonumber\\[4pt]& = {Z_h} \cdot \Delta X(s).\end{align}

This modeling can also be called an impedance modeling. The joystick occupies a 3D space, and thus ${F_h}$ and $\Delta X \in\square^{3\times 1}$ , ${Z_h},{M_h},{B_h},$ and $K_h \in\square^{3\times 3}$ .

4.1. Human arm impedance estimation

To determine the impedance coefficients of the human arm, the operator works in the position command mode. At this time, ${X_d}$ is a preset reference path and only the operator arm’s impedance coefficients in this reference path can be estimated because the arm’s impedance changes with muscle contractions [Reference Atsushi and Shohei26]. F/T sensor is used to measure the interaction force between the human hand and the joystick. Based on the interaction force between the hand and joystick measured by the F/T sensor and the position increment calculated by the angles measured by the encoders, the impedance coefficients of the human arm can be obtained as follows:

(13) \begin{align}{Z_h} = \frac{F_h}{\Delta X}\end{align}

If ${F_h}(x,y,z)$ can be expressed as a differentiable nonlinear function of the position of the end point of the human arm as follows:

(14) \begin{align}d{F_i} = \left(\frac{{\partial {F_i}}}{{\partial x}}\right)dx + \left(\frac{{\partial {F_i}}}{{\partial y}}\right)dy + \left(\frac{{\partial {F_i}}}{{\partial z}}\right)dz\end{align}

where i = $x,y,$ and $z$ which are the coordinates of the end point of the human arm in a 3D space. Notice that the human arm’s impedance can be expressed as a differential operator and is associated with infinitesimal force changes and infinitesimal displacements.

Equation (14) is correct only in the case of infinitesimal displacements. Therefore, a reliable method is measuring the force of many discontinuous points and then estimating the coefficients of the human arm’s impedance modeling by processing the obtained data.

The impedance modeling proposed in the previous section is represented as follows:

(15) \begin{align}F_h = M_h \left(\ddot{X}_c-\ddot{X}_d\right) + B_h \left(\dot{X}_c-\dot{X}_d\right) + K_h \left({X}_c-{X}_d\right)\end{align}

where ${X_d}$ and ${X_c}$ represent the desired and actual positions, respectively; ${\dot X_c}$ and ${\ddot X_c}$ are the actual velocity and acceleration of the arm, respectively; ${\dot X_d}$ and ${\ddot X_c}$ are the desired velocity and acceleration of the arm, respectively.

Assuming that the operator’s hand applies a force $d{F_h}$ to the joystick over a very short time, which causes the joystick to move a very short distance $d{X_h}$ , the impedance modeling can be defined as follows:

(16) \begin{align}F_h + dF_h = M_h (\ddot{X}_c + d\ddot{X}_h - \ddot{X}_d) + B_h(\dot{X}_c + d\dot{X}_h - \ddot{X}_d) + K_h({X}_c + d{X}_h - {X}_d),\end{align}

where $dX_h, d\dot{X}_h$ , and $d\ddot{X}_h$ are the distance, velocity, and acceleration of the arm’s motion over a short time, respectively. By subtracting Eq. (15) from Eq. (16), $d{F_h}$ can be expressed as follows:

(17) \begin{align}dF_h = M_h d\ddot{X}_h + B_h d\dot{X}_h + K_h d{X}_h.\end{align}

Assuming that ${\hat F_h}$ is the estimated force output of the human hand, ${\hat M_h},{\hat B_h}, $ and ${\hat K_h}$ are the estimated values of the arm’s mass, damping, and stiffness coefficients, respectively. The human arm’s force, ${\hat F_h}$ , can be represented as follows:

(18) \begin{align}\hat{F}_h = \hat{M}_h \ddot{X} + \hat{B}_h \dot{X} + \hat{K}_h {X}.\end{align}

By using the least squares method, square error S can be expressed as follows:

(19) \begin{align}S = {\sum\limits_{p = 1}^n {\left[ {{F_h}(p) - {{\hat F}_h}(p)} \right]} ^2}\end{align}

Notice that $\frac{{\partial S}}{{\partial {{\hat M}_h}}}$ , $\frac{{\partial S}}{{\partial {{\hat B}_h}}}$ , and $\frac{{\partial S}}{{\partial {{\hat K}_h}}}$ can be approximated as zeros [Reference Li and Song22]. Using the Eq. (18) and (19), impedance coefficients of human arm can be estimated as follows:

(20) \begin{align}\left[\begin{array}{c}\hat{M}_h\\[4pt]\hat{B}_h\\[4pt]\hat{K}_h\\\end{array}\right] = \left[\begin{array}{c@{\quad}c@{\quad}c}\sum_{p=1}^{n} \ddot{X}^2 (p) & \sum_{p=1}^{n} \ddot{X}(p)\dot{X}(p) & \sum_{p=1}^{n} \ddot{X}(p){X}(p)\\[7pt]\sum_{p=1}^{n} \dot{X}(p)\ddot{X}(p) & \sum_{p=1}^{n} \dot{X}^2(p) & \sum_{p=1}^{n} \dot{X}(p){X}(p)\\[7pt]\sum_{p=1}^{n} {X}(p)\ddot{X}(p) & \sum_{p=1}^{n} {X}(p)\ddot{X}(p) & \sum_{p=1}^{n} {X}^2(p)\end{array}\right]^{-1} \cdot \left[\begin{array}{c}\sum_{p=1}^{n} \ddot{X}(p)F_h(p)\\[4pt]\sum_{p=1}^{n} \dot{X}(p)F_h(p)\\[4pt]\sum_{p=1}^{n} {X}(p)F_h(p)\end{array}\right]\end{align}

where n is the number of sampling points and t is the sampling time.

To estimate the time-varying impedance coefficients of the human arm in a short time, the impedance can be represented as follows:

(21) \begin{align}{\hat Z_h}(t) = V{(t)^{ - 1}}W(t)\end{align}

where $\left( t \right) = \left[ {\begin{array}{c@{\quad}c@{\quad}c}{\mathop \sum _{p = t - n + 1}^n {{\ddot X}^2}\left( p \right)} & {\mathop \sum _{p = t - n + 1}^n \ddot X\left( p \right)\dot X\left( p \right)} & {\mathop \sum _{p = t - n + 1}^n \ddot X\left( p \right)X\left( p \right)}\\[5pt]{ \sum _{p = t - n + 1}^n \dot X\left( p \right)\ddot X\left( p \right)} & { \sum _{p = t - n + 1}^n {{\dot X}^2}\left( p \right)} & { \sum _{p = t - n + 1}^n \dot X\left( p \right)X\left( p \right)}\\[5pt]{ \sum _{p = t - n + 1}^n X\left( p \right)\ddot X\left( p \right)} & { \sum _{p = t - n + 1}^n X\left( p \right)\ddot X\left( p \right)} & { \sum _{p = t - n + 1}^n {X^2}\left( p \right)}\end{array}} \right]$ and $W\left( t \right) = \left[ {\begin{array}{c}{ \sum _{p = t - n + 1}^n \ddot X\left( p \right){F_h}\left( p \right)}\\[5pt]{ \sum _{p = t - n + 1}^n \dot X\left( p \right){F_h}\left( p \right)}\\[5pt]{ \sum _{p = t - n + 1}^n X\left( p \right){F_h}\left( p \right)}\end{array}} \right]$ .

Notice that n is the number of samples and is set to 150 experimentally considering the control cycle. The elements of $V(t)$ and $W(t)$ can be estimated by the simulated annealing (SA) method as follows [Reference Song, Pan, Xu and Li27]:

(22) \begin{align} {v_{p,q}}(t + 1) & = \sum\limits_{i = t - n + 2}^n {{X^{(3 - p)}}(i)} {X^{(3 - q)}}(i)\nonumber\\[5pt]& = {v_{p,q}}(t) + {X^{(3 - p)}}(t + 1){X^{(3 - q)}}(t + 1) - {X^{(3 - p)}}(t - n + 1){X^{(3 - q)}}(t - n + 1),\end{align}

where $p = 1,2,3$ and $q = 1,2,3$ .

(23) \begin{align}{w_p}(t + 1) & = \sum\limits_{i = t - n + 2}^n {{X^{(3 - p)}}(i)} {F_h}(i) = {w_p}(t) + {X^{(3 - p)}}(t + 1){F_h}(t + 1) - {X^{(3 - p)}}(t - n + 1){F_h}(t - n + 1).\nonumber\\\end{align}

where $p = 1,2,3$ .

The resulting matrix can be approximated to a diagonal matrix to observe the changes in the human arm’s impedance.

4.2. Impedance control of joystick

The estimated human arm’s impedance has been matched to that of the joystick during the control of the joystick. This joystick used in this research is a parallel structure consisting of three serial manipulators, as shown in Fig. 4.

Each serial manipulator consists of two active joints and one passive joint; the gray joints are the active joints and the white joint is the passive joint. Three serial manipulators with 3-DOF are connected to the upper plate by passive spherical joints, and the motions of the serial manipulators generate the motion of the upper plate relative to the lower plate.

Figure 5 shows related angles in kinematic analysis of one of the serial manipulators. Because the joystick consists of three serial manipulators in parallel, the position of each joint point on the upper plate is identical to the position of the end effector of each serial manipulator. Therefore, if the positions of the end effectors of the three series manipulators are obtained, the coordinates of the central position can be obtained by finding the average of the coordinates of the three points. Thereafter, the position and direction of movement of the upper plate can be derived.

Table I. Denavit–Hartenberg parameter of each serial manipulator.

Figure 4. Joystick in this research.

Figure 5. Link parameters of a serial manipulator.

Table I shows the link parameter table for the serial manipulator in Fig. 5. The $a$ value of the DH parameter is length of the common normal, $d$ is offset along previous $z$ to the common normal, $\alpha $ is angle about common normal, from ${z_{i - 1}}$ c to ${z_i}$ and $\theta $ is angle about previous $z$ , from ${x_{i - 1}}$ to ${x_i}$ . The position kinematic equations for the manipulator can be obtained as follows:

(24) \begin{align}\left[ \begin{array}{l}{{\hat X}_i}\\[4pt]{{\hat Y}_i}\\[4pt]{{\hat Z}_i}\end{array} \right] = \left[ \begin{array}{l}{l_1}\cos ({\theta _{2i}}) + {l_2}\cos ({\theta _{2i}} + {\theta _{3i}})\\[4pt] \sin ({\theta _{1i}})({l_1}\sin ({\theta _{2i}}) + {l_2}\sin ({\theta _{2i}} + {\theta _{3i}}))\\[4pt] {l_1}\cos ({\theta _{1i}})\sin ({\theta _{2i}}) + {l_2}\cos ({\theta _{1i}})\sin ({\theta _{2i}} + {\theta _{3i}})\end{array} \right]\end{align}

The position of the haptic joystick’s end effector is considered to be the center of the upper plate and is marked as X. The prediction of the target position is based on the force in the x-y-z directions sensed by the F/T sensor and the previously estimated human arm’s impedance coefficients. The differential equation is solved and the required position increments in the x-y-z directions are obtained. This joystick can be controlled by dynamic modeling as follows:

(25) \begin{align}{M_h}\left( {{{\ddot X}_c} - {{\ddot X}_d}} \right) + {B_h}\left( {{{\dot X}_c} - {{\dot X}_d}} \right) + {K_h}\left( {{X_c} - {X_d}} \right) = F.\end{align}

With the definition of $\Delta X = {X_c} - {X_d}$ , which represents the position increment caused by human’s motion intention, Eq. (25) can be represented as

(26) \begin{align}\Delta \ddot{X} + M_h^{-1}B_h \Delta \ddot{X} + M_h^{-1}K_h \Delta {X} = M_h^{-1} F.\end{align}

The analytical solution of differential Eq. (26) can be obtained ignoring the initial conditions as follows:

(27) \begin{align}\Delta X = {e^{ - {\lambda _1}t}} + {e^{ - {\lambda _2}t}} + K_h^{ - 1}F \end{align}

where ${\lambda _1}$ and ${\lambda _2}$ are the general solutions of the following characteristic equation:

(28) \begin{align}{\lambda ^2} + M_h^{ - 1}{B_h}\lambda + M_h^{ - 1}{K_h} = 0\end{align}

To manipulate this joystick to follow the human hand compliantly, a PD control with the gravity compensation algorithm has been adopted to drive all the motor parameters in the joystick. The designed controller can be represented as

(29) \begin{align}{F_c} = K_p \Delta X + K_D \Delta \dot{X} + G(\theta)\end{align}

where ${K_P}$ and ${K_D}$ are the controller gains, $\Delta X$ is the position increment, and $G(\theta )$ is the gravity compensation force.

To obtain the gravity term, it is crucial to determine the mapping relationship between the gravity force and joint configuration. Considering that all joints of the joystick are revolute joints, if link $p$ ’s center of mass is located at a point in link $p$ , gravity $G(\theta )$ can be expressed as

(30) \begin{align}G(\theta ) = g(\theta )A\end{align}

where $g(\theta )$ is the gravity acceleration on each link, which can be expressed as follows:

(31) \begin{align}g(\theta ) = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{{g_{11}}(\theta )} & {{g_{12}}(\theta )} & \ldots & {{g_{1(n - 1)}}(\theta )} & {{g_{1n}}(\theta )}\\[5pt]0 & {{g_{22}}(\theta )} & \ldots & {{g_{2(n - 1)}}(\theta )} & {{g_{2n}}(\theta )}\\[5pt]\ldots & \ldots & \ldots & \ldots & \ldots\\[5pt]0 & 0 & \ldots & {{g_{(n - 1)(n - 1)}}(\theta )} & {{g_{(n - 1)n}}(\theta )}\\[5pt]0 & 0 & \ldots & 0 & {{g_{nn}}(\theta )}\end{array}} \right]\end{align}

Each element in matrix $g(\theta )$ can be calculated by

(32) \begin{align}{g_{pq}}(\theta ) = g({{\partial R_0^p} / {\partial {\theta _p}}}){u_p}\end{align}

where $g = {[0,0,9.8]^T}m/{s^2},$ $R_0^p$ is the rotation matrix of coordinate system $p$ with respect to coordinate system $0$ , ${u_p}$ is the coordinate system $p$ ’s directional unit vector that points to link $p$ ’s mass center.

In Eq. (30), $A$ , a constant matrix, which is represented as follows:

(33) \begin{align}A = \left[ {\begin{array}{c}{{m_1}l_1^c + {l_1}\sum\limits_{i = 2}^n {{m_i}} }\\[4pt]{{m_2}l_2^c + {l_2}\sum\limits_{i = 3}^n {{m_i}} }\\[4pt]{{m_3}l_3^c + {l_3}\sum\limits_{i = 4}^n {{m_i}} }\\[4pt]\ldots\\[4pt]{{m_{n - 1}}l_{n - 1}^c + {l_{n - 1}}{m_i}}\\[4pt]{{m_n}l_n^c}\end{array}} \right]\end{align}

where $l_p^c$ is the distance between coordinate system $p$ ’s origin and link $p$ ’s mass center, ${m_i}$ is the mass of link i, and ${l_p}$ is the distance between coordinate system $p$ ’s origin and coordinate system $p + 1$ ’s origin.

5.1. Human arm’s impedance estimation experiment

The following experiment is conducted to obtain the coefficients of the human arm’s impedance modeling. A subject is a 25-year-old man with no disabilities. The subject sits on the chair and holds the handle and moves the joystick along the target trajectory with only gravity compensation. The target trajectory is a square with a side length of 0.3 m, centered on the initial position of the handle. The four vertices of the square are A, B, C, and D, respectively. The data collected are used to estimate the impedance of the human arm.

The data are sampled and stored for offline analysis. The data include the joystick joint angles and the forces in the x–y–z directions. The position of the joystick’s end effector is calculated by considering the joint angles in the kinematic equation of the joystick. Details regarding the kinematics of the joystick are described in Section 5. The impedance coefficients obtained by the least squares method are shown in Fig. 6.

Figure 6. Estimated human arm’s impedance coefficients.

These figures show that the impedance of the human arm changes with the movement of the arm. Therefore, it is necessary to perform a real time impedance control to accurately match the desired impedance to the human arm’s impedance during the joystick control.

Table II. Interaction force in the x, y and z directions before and after using the impedance matching algorithm.

5.2. Human motion prediction and reaction experiment

The estimated human arm’s impedance coefficient is applied to the joystick impedance control to confirm the efficiency. The operator holds the handle and moves the joystick along the target trajectory. During this process, the force applied by the operator to the joystick is included in the impedance modeling to determine the target position, which is the target position at which the human brain sends a motion command to the arm through the nerve. The target position is a motion vector indicating the position increments in the x, y, and z directions. The designed PD controller with gravity compensation drives the motor settings in the joystick to allow the end effector to follow the guidance or traction of the operator’s hand. The F/T sensor disposed at the end effector of the joystick senses the interaction force between the operator’s hand and the joystick. The position of the joystick’s end effector is transmitted to the mobile robot through the communication module to control the movement of the mobile robot. The position of the joystick’s end effector has been calculated using the encoder recording angle and the link length, which has been used to calculate the trajectory tracking error against the reference trajectory.

Figure 7 shows the interaction forces between the operator’s hand and the joystick before and after applying the impedance matching algorithm. It can be observed that the average values of the interaction force in the x, y, z directions become very small after applying the impedance matching algorithm. The specific values are shown as follows:

Table II shows the average interaction forces in x, y, z directions after applying the impedance matching algorithm, respectively. Notice that the interaction force becomes almost 7 times smaller by applying the impedance matching algorithm. The maximum forces in x, y, z directions after applying the impedance matching algorithm become 4 to 5 times smaller than before.

Figure 8 shows the desired trajectory calculated from the estimated human arm’s impedance modeling and the actual trajectory in the impedance matching algorithm. It can be quantified that the average values of the trajectory tracking error are shown in Table III. In the experiment, there are three paths: square reference path (represented by A, B, C and D points), desired trajectory calculated by impedance model and force and real trajectory actuated by motors in the joystick. Through Figure 8 and Table III, the error values between the desired trajectory and the reference path occurred as much as 0.015 m, 0.012 m, and 0.011 m.

Table III. Trajectory tracking accuracy.

Figure 7. Interaction forces before and after applying impedance matching algorithm.

Figure 8. Trajectory tracking performance.

By the impedance matching control, the interaction force between the human and the joystick becomes 5 times smaller with the tracking error of 0.02 m compared to the PD control without impedance matching. This proves the effectiveness of the human arm’s impedance matching algorithm proposed in this research.