2.1. Admittance system model

The collaborative motion of human–robot contact and human–human contact can be equivalent to the contact motion between two mass-damping-spring model systems. The equivalent model of the active force applicator (leader) is the impedance system, which outputs the force according to the input position. The equivalent model of the follower is the admittance system, which outputs the position according to the input force. Both systems show the dynamic relationship between position and interaction force, with equal parameters and the same formula [Reference Kang, Oh, Seo, Kim and Choi17, Reference Žlajpah and Petrič28]:

(1) \begin{equation} \boldsymbol{{C}}({\ddot{\boldsymbol{{X}}}_{d}}-{\ddot{\boldsymbol{{X}}}_{r}})+\boldsymbol{{C}}({\dot{\boldsymbol{{X}}}_{d}}{-}{\dot{\boldsymbol{{X}}}_{r}})+\boldsymbol{{K}}({{\boldsymbol{{X}}}_{d}}{-}{{\boldsymbol{{X}}}_{r}})=\boldsymbol{{F}} \end{equation}

where $\boldsymbol{{M}}\in{{R}^{3\times 3}}$ , ${\boldsymbol{{C}}}\in{{R}^{3\times 3}}$ , and ${\boldsymbol{{K}}}\in{{R}^{3\times 3}}$ are the virtual mass matrix, virtual damping matrix, and virtual matrix, respectively; ${\boldsymbol{{F}}}={{[{{f}_{x}},{{f}_{y}},{{f}_{z}}]}^{T}}$ is the cartesian space interaction force vector; ${{{\boldsymbol{{X}}}}_{d}}=[{{x}_{d}},{{y}_{d}},{{z}_{d}}]^T$ is the desired position vector; and the current position of robot is set as the reference position vector ${{\boldsymbol{{X}}}_{r}}=[{{x}_{r}},{{y}_{r}},{{z}_{r}}]^T$ .

For convenient applications in the control of industrial robots, Eq. (1) can be converted to its discrete format to solve the desired trajectory of the robot [Reference Lecours, Mayer-St-Onge and Gosselin29, Reference Duan, Gan, Chen and Dai30]:

(2) \begin{equation}{{{\ddot{\boldsymbol{{X}}}}}_{{d}}}(t)={{\boldsymbol{{M}}}^{{-1}}}({\boldsymbol{{F}}-\boldsymbol{{C}}\dot{\tilde{\boldsymbol{{X}}}}-\boldsymbol{{K}}\tilde{\boldsymbol{{X}}}})+{{{\ddot{\boldsymbol{{X}}}}}_{r}}(t) \end{equation}

(3) \begin{equation}{{{\dot{\boldsymbol{{X}}}}}_{{d}}}(t)={{{\dot{\boldsymbol{{X}}}}}_{{r}}}(t)+T{{{\ddot{\boldsymbol{{X}}}}}_{{d}}}(t)\ \end{equation}

(4) \begin{equation}{{\boldsymbol{{X}}}_{{d}}}(t){=}{{\boldsymbol{{X}}}_{{r}}}(t){+}T{{{\dot{\boldsymbol{{X}}}}}_{{d}}}(t)\ \end{equation}

where ${{\boldsymbol{{X}}}_{{d}}}(t)$ , ${{{\ddot{\boldsymbol{{X}}}}}_{{d}}}(t)$ are the desired position, desired velocity, and acceleration at moment t, respectively; ${\tilde{\boldsymbol{{X}}}}={{{\boldsymbol{{X}}}}_{{r}}}(t)-{{{\boldsymbol{{X}}}}_{{r}}}(t-1)$ ; ${\tilde{\boldsymbol{{X}}}}={{\boldsymbol{{X}}}_{{r}}}(t)-{{\boldsymbol{{X}}}_{{r}}}(t-1)$ ; and T is the system communication cycle between the robot controller and the servo driver.

In traditional robot admittance control, admittance parameters are fixed and pre-set according to the environment. In the human–robot collaborative motion, the human is the dynamic system, and fixed admittance parameters are insufficient to achieve natural and compliant human–robot collaborative motion. To achieve better human–robot collaborative performance, the admittance parameters of the robot should be changed according to the change of the collaborative behavior to ensure a high matching degree between the admittance parameters of the robot and the collaborative motion of the human. The best method to improve the matching degree between robot admittance parameters and human collaborative motion is to imitate the change rule of follower admittance parameters in the human–human collaborative motion [Reference Roveda, Vicentini and Tosatti9, Reference Ikeura, Moriguchi and Mizutani27]. Therefore, a human–human collaborative experiment is designed to study the variation of follower admittance parameters in the human–human collaborative motion.

2.2. Human–human collaborative motion experiment

The cartesian space trajectory of humans and robots can be synthesized by the motion trajectory of the X, Y, and Z axes. Therefore, the human–human collaborative motion experiments along the X, Y, and Z axes were designed, as shown in Fig. 1(a)–(c). Figure 1(d) is the experimental scene diagram. The average arm length of all participants is 57.2 cm. So, the maximum motion distance along the X-axis is 114.4 cm. The maximum motion distance along the Y-axis is 57.2 cm, and the maximum motion distance along the Z-axis is 114.4 cm. The combination with the common workspace of the human arm in human engineering, the starting and target points along the X-axis, Y-axis, and Z-axis, is set. The human–human collaborative motion experiment’s starting point and target point along the X, Y, and Z axes are in the workspace that can be reached by the human arm. The distance between adjacent points is marked in Fig. 1(a)–(c). The human–human collaborative motion experiments along the X-axis and Y-axis are shown in Fig. 1(a) and (b), respectively. The starting point of the human–human collaborative motion is the blue square in the figure, and the target is the red circle in the figure. The human–human collaborative motion experiment along the X-axis has 3 starting points and 12 target points. The human–human collaborative motion experiment along the Y-axis has 5 starting points and 10 target points. The human–human collaborative motion experiment along the Z-axis is shown in Fig. 1(c). The starting point of the human–human collaborative motion is the blue cuboid in the figure, with 15 starting points in total, and the target point is the red cube in the figure. Three target point regions (C1, C2, and C3) are set along the Z-axis, with 15 target points in each layer. The telescopic rod on the bracket in Fig. 1(c) is used to mark the target point position of each human–human collaborative motion.

Figure 1. The human–human collaborative motion experiment scheme.

In the experiment, the Onrobot force sensor (HEX-E, Onrobot, Denmark) collects the interaction force of the human–human collaborative motion. Both ends of the force sensor are connected with a handle, and participants A and B hold one handle. The follower (participant B) completes the human–human collaborative motion from the starting point to the target point under the leader’s guidance (participant A). The red arrow in the figure indicates the motion direction. In point-to-point motion, the leader tries to keep the guiding path straight. The follower does not need to know the exact motion path but only actively follows the motion of the leader by their sense of displacement, velocity, and force, leaving an impression of obedience on the leader. During the experiment, the followers are required to wear motion capture clothing (Perception Legacy, Noitom, China) to collect data such as motion position.

A total of five groups of participants (two in each group) participated in the human–human collaborative motion experiment along the X, Y, and Z axes. The participants (all males, with mean $age=27\pm 3$ years, ${{height}}=175\pm 6$ mm, and $weight=71.4\pm 6.8$ kg) are right handed and have no physical problems. The same point-to-point motion experiment was repeated 20 times. Therefore, 1200 $(5\times 20\times 12=1200)$ , 1000 $(5\times 20\times 10=1000)$ , and 4500 $(5\times 20\times 15\times 3=4500)$ instances of motion data were collected along the X, Y, and Z axes, respectively. Both the force sensor and the motion-capture clothing collect data at 100 Hz. Before the experiment, each group of participants performed a preliminary experiment to familiarize themselves with the experiment’s content.

3.1. Calculation method of admittance parameters of human arm

Studies show that the admittance parameters of the human arm do not change significantly in a brief time (10 ms) during normal motion, that is, the admittance parameter values at adjacent moments are similar [Reference Li and Ge31]. Therefore, this paper calculates admittance parameters based on the interaction force and position data of the follower at three adjacent moments in the human–human collaborative motion. According to Eq. (1), admittance models at three adjacent moments can be obtained:

(5) \begin{equation} \left [ \begin{array}{c@{\quad}c@{\quad}c}{{{{\ddot{\boldsymbol{{X}}}}}}_{e}}(t-1) &{{{{\dot{\boldsymbol{{X}}}}}}_{e}}(t-1) &{{\boldsymbol{{X}}}_{e}}(t-1) \\ {{{{\ddot{\boldsymbol{{X}}}}}}_{e}}(t) &{{{{\dot{\boldsymbol{{X}}}}}}_{e}}(t) &{{{\boldsymbol{{X}}}}_{e}}(t) \\ {{{{\ddot{\boldsymbol{{X}}}}}}_{e}}(t+1) &{{{{\dot{\boldsymbol{{X}}}}}}_{e}}(t+1) &{{\boldsymbol{{X}}}_{e}}(t+1) \\ \end{array} \right ]\left [ \begin{matrix} \boldsymbol{{M}}(t) \\ {\boldsymbol{{C}}}(t) \\ {\boldsymbol{{K}}}(t) \\ \end{matrix} \right ]=\left [ \begin{matrix} \boldsymbol{{F}}(t-1) \\ \boldsymbol{{F}}(t) \\ \boldsymbol{{F}}(t+1) \\ \end{matrix} \right ] \end{equation}

where ${{\boldsymbol{{X}}}_{e}}(t)={{\boldsymbol{{X}}}_{{d}}}(t)-{{\boldsymbol{{X}}}_{{r}}}(t)$ , ${{{\boldsymbol{{X}}}}_{{d}}}(t)$ , and ${{\boldsymbol{{X}}}_{{r}}}(t)$ are the desired position and actual position at moment t, respectively; $\boldsymbol{{M}}(t)$ , ${\boldsymbol{{C}}}(t)$ , and ${\boldsymbol{{K}}}(t)$ are virtual mass matrix, virtual damping matrix, and virtual stiffness matrix at moment t, respectively, and $\boldsymbol{{F}}(t)$ is the interaction force at moment t.

According to Eq. (5), the admittance parameters ${\boldsymbol{{M}}}(t)$ , ${\boldsymbol{{C}}}(t)$ , and ${\boldsymbol{{K}}}(t)$ , at moment t can be obtained:

(6) \begin{equation} \left [ \begin{matrix} \boldsymbol{{M}}(t) \\ {\boldsymbol{{C}}}(t) \\ {\boldsymbol{{K}}}(t) \\ \end{matrix} \right ]={{\left [ \begin{array}{c@{\quad}c@{\quad}c}{{{{\ddot{\boldsymbol{{X}}}}}}_{e}}(t-1) &{{{{\dot{\boldsymbol{{X}}}}}}_{e}}(t-1) &{{\boldsymbol{{X}}}_{e}}(t-1) \\ {{{{\ddot{\boldsymbol{{X}}}}}}_{e}}(t) &{{{{\dot{\boldsymbol{{X}}}}}}_{e}}(t) &{{\boldsymbol{{X}}}_{e}}(t) \\ {{{{\ddot{\boldsymbol{{X}}}}}}_{e}}(t+1) &{{{{\dot{\boldsymbol{{X}}}}}}_{e}}(t+1) &{{\boldsymbol{{X}}}_{e}}(t+1) \\ \end{array} \right ]}^{-1}}\left [ \begin{matrix} \boldsymbol{{F}}(t-1) \\ \boldsymbol{{F}}(t) \\ \boldsymbol{{F}}(t+1) \\ \end{matrix} \right ] \end{equation}

By substituting the experimental data of the human–human collaborative motion from Section 2.2 into Eq. (6), the admittance parameters of the follower in the process of the human–human collaborative motion can be obtained.

3.2. Fitting of admittance parameters of human arm

The admittance parameters of the human arm were divided into the static and dynamic admittance parameters, which reflect the static and dynamic relationship between the end position and the interaction force of the human arm when the arm is in the motion preparation phase and the motion phase, respectively [Reference Perreault, Kirsch and Crago32]. In the human–human collaborative motion, the follower undergoes a transition from the state of preparation to the state of motion, so the admittance parameters of the follower change from static to dynamic admittance parameters.

In the first three frames of the human–human collaborative motion, the position and velocity change are small, so the human arm is in the motion preparation phase during this period. By substituting the data of the first three frames of the human–human collaborative motion into Eq. (6), the static admittance parameters of the human arm in the preparatory phase of motion can be obtained. The relationship between the static admittance parameters and the interaction forces was obtained from all experimental data, as shown in Fig. 2.

Figure 2. The relationship between the static admittance parameters and the interaction forces of the human–human collaborative motion. ${f}_{sx}$ , ${f}_{sy}$ , and ${f}_{sz}$ are the interaction forces of the X, Y, and Z axes, respectively; ${m}_{sx}$ , ${m}_{sy}$ , and ${m}_{sz}$ are the static virtual masses of the X, Y, and Z axes, respectively; ${c}_{sx}$ , ${c}_{sy}$ , and ${c}_{sz}$ are the static virtual damping of the X, Y, and Z axes, respectively; and ${k}_{sx}$ , ${k}_{sy}$ , and ${k}_{sz}$ are the static virtual stiffnesses of the X, Y, and Z axes, respectively.

By substituting the experimental data of the human–human collaborative motion into Eq. (6), the dynamic admittance parameters of the human arm in the motion phase can be obtained. Through the analysis of the experimental data, it is found that the change rate of admittance parameter changes with the change rate of interaction force, showing a certain regularity. Therefore, the relationship between the admittance parameter’s change rate and the interactive force’s change rate of the X-axis, Y-axis, and Z-axis was statistically analyzed. The relationship between the dynamic admittance parameter’s change rate and the change rate of interaction force was obtained from all experimental data, as shown in Fig. 3.

Figure 3. The relationship between the change rate of dynamic admittance parameter and the change rate of interaction force of the human–human collaborative motion. ${\dot{f}}_{x}$ , ${\dot{f}}_{y}$ , and ${\dot{f}}_{z}$ are the interaction forces’ change rates of the X, Y, and Z axes, respectively; ${\dot{m}}_{x}$ , ${\dot{m}}_{y}$ , and ${\dot{m}}_{z}$ are the dynamic virtual masses’ change rates of the X, Y, and Z axes, respectively; ${\dot{c}}_{x}$ , ${\dot{c}}_{y}$ , and ${\dot{c}}_{z}$ are the dynamic virtual damping’ change rates of the X, Y, and Z axes, respectively; and ${\dot{k}}_{x}$ , ${\dot{k}}_{y}$ , and ${\dot{k}}_{z}$ are the dynamic virtual stiffnesses’ change rates of the X, Y, and Z axes, respectively.

Figures 2 and 3 show that the static and dynamic admittance parameters show certain regularity. To apply the variation rule of admittance parameters to the robot admittance parameter regulator, the relationship between the static admittance parameters and interaction force, as well as the relationship between the change rate of dynamic admittance parameter and the change rate of interaction force, was fitted by the data fitting method. According to the experiment in Section 2, the human–human collaborative motion experiments along the X-axis, Y-axis, and Z-axis were performed 1200, 1000, and 4500 times, respectively. The average duration of each collaborative motion is 160 frames (1 s is equal to 100 frames), so there are enough data for function fitting. According to observations derived from the scatter diagram in Figs. 2 and 3, the data fitting function can choose the polynomial function. According to the observation of the shape of the curve of the polynomial function and the observation of the shape of the scatter plots in Figs. 2 and 3, it is known that the fitting times (the order of the polynomial) of fitting functions do not exceed six times.

The polynomial function H of degree h is

(7) \begin{equation} H(x,p)=\sum \limits _{j=0}^{h}{{{p}_{j}}{{x}^{j}}} \end{equation}

where x is the input and p is the parameter of the polynomial function with h + 1 parameters and $h\le 6$ .

To determine the complexity of the model, that is, to determine the number of polynomials and prevent over-fitting, the regularization method is used to solve the polynomial parameters. The regularized cost function $L(p)$ is

(8) \begin{equation} L(p)=\frac{1}{N}\sum \limits _{i=1}^{N}{{{(H({{x}_{i}},p)-{{y}_{i}})}^{2}}+\frac{\lambda }{2}}{{\left \| p \right \|}_{2}} \end{equation}

where the first term is the average loss function, the second term is the regularization term, $\lambda \ge 0$ is the coefficient that adjusts the relationship between the two, N is the number of fitting data, ${x}_{i}$ is the observed value of the input x, that is, the observed value of the interaction force or the observed value of the interaction force’s change rate, ${y}_{i}$ is the observed value of the input y, that is, the observed value of the admittance parameter or the observed value of the admittance parameter’s change rate, and ${ \| p \|}_{2}$ represents the L2 norm of a polynomial function parameter vector.

By substituting the X-axis data in Fig. 2 into Eq. (8), the fitting function of each static admittance parameter and interaction force when the human–human collaboration moves along the X-axes can be obtained by minimizing the regularized cost function:

(9) \begin{equation} \left [ \begin{matrix}{{m}_{sx}} \\ {{c}_{sx}} \\ {{k}_{sx}} \\ \end{matrix} \right ]=\left [ \begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{p}_{m\text{sx}1}} &{{p}_{m\text{sx}2}} &{{p}_{m\text{sx}3}} &{{p}_{m\text{sx}4}} \\ {{p}_{csx1}} &{{p}_{csx2}} &{{p}_{csx3}} &{{p}_{csx4}} \\ {{p}_{ksx1}} &{{p}_{ksx2}} &{{p}_{ksx3}} &{{p}_{ksx4}} \\ \end{array} \right ]\left [ \begin{matrix}{{f}_{sx}} \\ f_{sx}^{2} \\ f_{sx}^{3} \\ f_{sx}^{4} \\ \end{matrix} \right ]+\left [ \begin{matrix}{{p}_{msx0}} \\ {{p}_{csx0}} \\ {{p}_{ksx0}} \\ \end{matrix} \right ] \end{equation}

where ${m}_{sx}$ , ${c}_{sx}$ , and ${k}_{sx}$ represent the virtual mass, virtual damping, and virtual stiffness of the static admittance parameters, respectively, and ${f}_{sx}$ is the interaction force of X-axis.

Table I shows the parameters p of the fitting function of each static admittance parameter and interaction force when the human–human collaborative moves along the X-axes. The fitting function between the static admittance parameters of the Y-axis and Z-axis and the interaction forces of the Y-axis and Z-axis is obtained by similar methods, which will not be described here.

Table I. Parameters of the fitting function of X-axes static virtual mass, static virtual damping, and static virtual stiffness.

By substituting the X-axis data in Fig. 3 into Eq. (8), the fitting function of the dynamic admittance parameters’ change rate and the interaction forces’ change rate when human–human collaboration moves along the X-axes can be obtained by minimizing the regularized cost function:

(10) \begin{equation} \left [ \begin{matrix}{{{\dot{m}}}_{x}} \\ {{{\dot{c}}}_{x}} \\ {{{\dot{k}}}_{x}} \\ \end{matrix} \right ]=\left [ \begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{p}_{m\text{x}1}} &{{p}_{m\text{x}2}} &{{p}_{m\text{x}3}} &{{p}_{m\text{x}4}} \\ {{p}_{cx1}} &{{p}_{cx2}} &{{p}_{cx3}} &{{p}_{cx4}} \\ {{p}_{kx1}} &{{p}_{kx2}} &{{p}_{kx3}} &{{p}_{kx4}} \\ \end{array} \right ]\left [ \begin{matrix}{{{\dot{f}}}_{x}} \\ \dot{f}_{x}^{2} \\ \dot{f}_{x}^{3} \\ \dot{f}_{x}^{4} \\ \end{matrix} \right ]+\left [ \begin{matrix}{{p}_{mx0}} \\ {{p}_{cx0}} \\ {{p}_{kx0}} \\ \end{matrix} \right ] \end{equation}

where ${\dot{m}}_{x}$ , ${\dot{c}}_{x}$ , and ${\dot{k}}_{x}$ represent the virtual mass change rate, virtual damping change rate, and virtual stiffness change rate of the dynamic admittance parameters, respectively, and ${\dot{f}}_{x}$ is the change rate of admittance parameter.

Table II shows the parameter p of the fitting function of the change rate of dynamic admittance parameter and the change rate of the interaction force when the human–human collaborative moves along the X-axes. The fitting function between the dynamic admittance parameters change rate of the Y-axis and Z-axis, and the interaction forces change rate of the Y-axis and Z-axis is obtained by similar methods, which will not be described here.

Table II. Parameters of the fitting function of X-axes dynamic virtual mass change rates, dynamic virtual damping change rates, and dynamic virtual stiffness change rates.

To assess the fitting effect of the fitting function, the goodness and root mean square error (RMSE) of fit of these polynomial fitting functions were calculated. The expressions of the goodness $R^{2}$ and RMSE of fit are as follows:

(11) \begin{equation}{{R}^{2}}=1-\frac{\sum \limits _{i=1}^{N}{{{({{{\hat{y}}}_{i}}-{{y}_{i}})}^{2}}}}{\sum \limits _{i=1}^{N}{{{({{{\bar{y}}}_{i}}-{{y}_{i}})}^{2}}}} \end{equation}

(12) \begin{equation} \text{RMSE}=\sqrt{\frac{1}{N}\sum \limits _{i=1}^{N}{{{({{{\hat{y}}}_{i}}-{{y}_{i}})}^{2}}}} \end{equation}

where ${y}_{i}$ represents the actual value of the data, ${\hat{y}}_{i}$ represents the fitted value of the data, ${\bar{y}}_{i}$ represents the mean value of the actual value of the data, and N represents the number of data. A value of ${R}^{2}$ closer to 1 indicates a better fit, and a smaller RMSE indicates a smaller deviation between the fitted value and the actual value.

The results of the goodness and RMSE of fit calculated for each fitting function are shown in Tables III and IV. Tables III and IV show that the goodness of fit of the polynomial fitting function is all greater than 0.9, and the maximum RMSE is 4.138, indicating that the fitting curve has a good fitting effect and the degree of data dispersion is small. Therefore, the fitting curve can better represent the regularity of data.

Table III. Goodness and RMSE of fit of the fitting functions for X, Y, and Z axes static admittance parameters.

Table IV. Goodness and RMSE of fit of the fitting function of change rates of X, Y, and Z axes dynamic admittance parameters.

3.3. Determination of admittance parameters of the human arm

According to the definition of the static and dynamic admittance parameters, the dynamic admittance parameter function can be obtained by adding the static admittance parameter function and change rate function of the dynamic admittance parameters:

(13) \begin{equation} \boldsymbol{{M}}(t)={{\boldsymbol{{M}}}_{{s}}}+\sum \limits _{1}^{t}{{\dot{\boldsymbol{{M}}}}(t)} \end{equation}

(14) \begin{equation}{\boldsymbol{{C}}}(t)={{{\boldsymbol{{C}}}}_{{s}}}+\sum \limits _{1}^{t}{{\dot{\boldsymbol{{C}}}}(t)} \end{equation}

(15) \begin{equation}{\boldsymbol{{K}}}(t)={{{\boldsymbol{{K}}}}_{{s}}}+\sum \limits _{1}^{t}{{\dot{\boldsymbol{{K}}}}(t)} \end{equation}

where $\boldsymbol{{M}}(t)=\text{diag}[{{m}_{x}}(t),{{m}_{y}}(t),{{m}_{z}}(t)]$ , ${{\boldsymbol{{M}}}_{{s}}}=\text{diag}[{{m}_{sx}},{{m}_{sy}},{{m}_{sz}}]$ , ${\boldsymbol{{C}}}(t)=\text{diag}[{{c}_{x}}(t),{{c}_{y}}(t),{{c}_{z}}(t)]$ ; ${{{\boldsymbol{{C}}}}_{{s}}}=\text{diag}[{{c}_{sx}},{{c}_{sy}},{{c}_{sz}}]$ , ${\boldsymbol{{K}}}(t)=\text{diag}[{{k}_{x}}(t),{{k}_{y}}(t),{{k}_{z}}(t)]$ , ${{{\boldsymbol{{K}}}}_{{s}}}= \text{diag}[{{k}_{sx}},{{k}_{sy}},{{k}_{s\text{z}}}]$ , ${\boldsymbol{{M}}}_{{s}}$ , ${{\boldsymbol{{C}}}}_{{s}}$ , and ${{\boldsymbol{{K}}}}_{{s}}$ are functions of interacting forces, which can be obtained from Eq. (9), ${\dot{\boldsymbol{{M}}}}(t)$ , ${\dot{\boldsymbol{{C}}}}(t)$ , and functions of the change rate of interaction force, which can be obtained from Eq. (10), and t is the motion time, in frames (1 s is equal to 100 frames). In general, the three-axis admittance parameters of the human arm play a leading role in the motion of the human arm [Reference Bascetta and Ferretti33]. Therefore, these admittance parameter matrices $\boldsymbol{{M}}, \boldsymbol{{C}}$ , and $\boldsymbol{{K}}$ are set as diagonal matrices [Reference Duchaine and Gosselin34].

According to the expression of admittance parameters of the human arm, the human arm can adjust the admittance parameters automatically and in real-time according to the change rate of interaction force, thus ensuring compliance in the human–human collaborative motion.

4.1. Robot control architecture

Based on the analysis of the various rules of admittance parameters of the human–human collaborative motion, the human–robot admittance control architecture is proposed, as shown in Fig. 4. In the human–robot collaborative motion, the interaction force is measured by a force sensor, and the differentiator obtains the change rate of interaction force. The human-like admittance parameter regulator takes the interaction force and the change rate of interaction force as an input and outputs admittance parameters matching the current human motion. The force measured by the force sensor and the admittance parameter output by the human-like admittance parameter regulator serves as the input of the human-like admittance controller. The human-like variable admittance controller can obtain the desired position, desired velocity, and desired acceleration according to Eqs. (2), (3), and (4), respectively. The robot’s current position, current velocity, and current acceleration are input into the robot together with the desired position, desired velocity, and desired acceleration of the leader. The robot moves according to the input information and returns the robot’s end position, velocity, and acceleration information in real time.

Figure 4. Block diagram of human-like variable admittance control in the human–robot collaborative motion system.

4.2. Controllability and stability analysis

To facilitate the controllability and stability analysis of the system, Eq. (1) is written as a state-space expression:

(16) \begin{equation} \left. \begin{aligned} &{{{{\dot{\boldsymbol{{X}}}}}}_{e}}=\boldsymbol{{A}}{{\boldsymbol{{X}}}_{e}}+{\boldsymbol{{B}}\boldsymbol{{F}}} \\ &{\boldsymbol{{Y}}}={\boldsymbol{{D}}}{{\boldsymbol{{X}}}_{e}} \\ \end{aligned} \right \} \end{equation}

where ${\boldsymbol{{A}}}= \left[ \begin{matrix}{\textbf{0}} &{\boldsymbol{{I}}} \\{-\boldsymbol{{K}}}{{{\boldsymbol{{M}}}}^{-1}} & -{\boldsymbol{{C}}}{{\boldsymbol{{M}}}^{-1}} -{\boldsymbol{{C}}}{{\boldsymbol{{M}}}^{-1}} \\ \end{matrix} \right]$ , ${{\boldsymbol{{X}}}_{e}}={{ [ \begin{matrix}{{\boldsymbol{{X}}}_{e1}} &\quad {{\boldsymbol{{X}}}_{e2}} \\ \end{matrix} ]}^{T}}$ , ${{\boldsymbol{{X}}}_{e1}}={{\boldsymbol{{X}}}_{{d}}}{-}{{\boldsymbol{{X}}}_{{r}}}$ , ${{{\boldsymbol{{X}}}}_{e2}}={{{\dot{\boldsymbol{{X}}}}}_{e1}}$ , ${\boldsymbol{{B}}}={{ [ \begin{matrix}{\textbf{0}} &\quad {{\boldsymbol{{M}}}^{-1}} \\ \end{matrix} ]}^{T}}$ , ${\boldsymbol{{D}}}={{ [ \begin{matrix}{\boldsymbol{{I}}} &\quad {\textbf{0}} \\ \end{matrix} ]}^{T}}$ . For simplicity, $\boldsymbol{{M}}, \boldsymbol{{C}}, \boldsymbol{{K}}$ , and $\boldsymbol{{F}}$ are substituted for $\boldsymbol{{M}}({t}), \boldsymbol{{C}}({t}), \boldsymbol{{K}}({t})$ , and $\boldsymbol{{F}}({t})$ . The meanings of $\boldsymbol{{M}}, \boldsymbol{{C}}, \boldsymbol{{K}}$ , and $\boldsymbol{{F}}$ are the same as those of $\boldsymbol{{M}}, \boldsymbol{{C}}, \boldsymbol{{K}}$ , and $\boldsymbol{{F}}$ in Eq. (1).

The controllability analysis. By the Popov–Belevitch–Hautus test [Reference Hespanha35], we have

(17) \begin{equation}{\boldsymbol{{Q}}=}\left [ \begin{matrix}{\boldsymbol{{B}}} &\quad \boldsymbol{{AB}} \\ \end{matrix} \right ]=\left [ \begin{array}{c@{\quad}c}{0} &{{\boldsymbol{{M}}}^{{-1}}} \\ {{\boldsymbol{{M}}}^{{-1}}} &{-\boldsymbol{{C}}}{{\boldsymbol{{M}}}^{{-1}}}{{\boldsymbol{{M}}}^{{-1}}} \\ \end{array} \right ] \end{equation}

where $\boldsymbol{{Q}}$ is the controllability matrix.

The necessary and sufficient condition for the system to be controllable is that $\boldsymbol{{Q}}$ is full rank. Therefore, for $\boldsymbol{{Q}}$ to have full rank, $\boldsymbol{{M}}$ should satisfy:

(18) \begin{equation} \text{rank}(\boldsymbol{{M}})=3 \end{equation}

The stability analysis. Theorem 1: The condition of the system (14) BIBO (bounded input bounded output) stability is the existence of a diagonal matrix $W\in{{R}^{6\times 6}}$ (a positive definite matrix) and a real constant $\gamma \in R$ such that the following matrix inequality:

(19) \begin{equation}{{\boldsymbol{{A}}}^{T}}\boldsymbol{{W}}+\boldsymbol{{WA}}+\frac{1}{\gamma }\boldsymbol{{WB}}{{{\boldsymbol{{B}}}}^{T}}\boldsymbol{{W}}\lt{0} \end{equation}

Let the Lyapunov function be:

(20) \begin{equation}{\boldsymbol{{V}}(}t{)=\boldsymbol{{X}}}_{e}^{T}\boldsymbol{{W}}{{\boldsymbol{{X}}}_{e}} \end{equation}

The derivative with respect to $\boldsymbol{{V}}({t})$ :

(21) \begin{align} {\dot{\boldsymbol{{V}}}}\text{(}t\text{)}&={\dot{\boldsymbol{{X}}}}_{e}^{T}\boldsymbol{{W}}{{\boldsymbol{{X}}}_{e}}+\boldsymbol{{X}}_{e}^{T}{\boldsymbol{{W}}\dot{\boldsymbol{{X}}}}_{e}^{T} \nonumber \\[3pt] &={{\text{(}\boldsymbol{{A}}{{\boldsymbol{{X}}}_{e}}+\boldsymbol{{BF}}\text{)}}^{T}}\boldsymbol{{W}}{{\boldsymbol{{X}}}_{e}}+{\boldsymbol{{X}}}_{e}^{T}\boldsymbol{{W}}\text{(}\boldsymbol{{A}}{{\boldsymbol{{X}}}_{e}}+\boldsymbol{{BF}}\text{)} \nonumber \\[3pt] &={\boldsymbol{{X}}}_{e}^{T}\text{(}{{\boldsymbol{{A}}}^{T}}{\boldsymbol{{W}}+\boldsymbol{{WA}}}\text{)}{{\boldsymbol{{X}}}_{e}}+{{\boldsymbol{{F}}}^{T}}{{{\boldsymbol{{B}}}}^{T}}\boldsymbol{{W}}{{\boldsymbol{{X}}}_{e}}{+\boldsymbol{{X}}}_{e}^{T}\boldsymbol{{WBF}} \nonumber \\[3pt] &={\boldsymbol{{X}}}_{e}^{T}\text{(}{{\boldsymbol{{A}}}^{T}}{\boldsymbol{{W}}+\boldsymbol{{WA}}}\text{)}{{\boldsymbol{{X}}}_{e}}+{2\boldsymbol{{X}}}_{e}^{T}\boldsymbol{{WBF}} \end{align}

According to Young inequality, there exists a positive constant $\gamma$ , such that

(22) \begin{align} {\dot{\boldsymbol{{V}}}}\text{(}t\text{)}&\le \boldsymbol{{X}}_{e}^{T}\text{(}{{\boldsymbol{{A}}}^{T}}{\boldsymbol{{W}}+\boldsymbol{{WA}}}\text{)}{{\boldsymbol{{X}}}_{e}}+\frac{1}{\gamma }\boldsymbol{{X}}_{e}^{T}\boldsymbol{{WB}}{{{\boldsymbol{{B}}}}^{T}}\boldsymbol{{W}}{{\boldsymbol{{X}}}_{e}}+\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}} \nonumber\\[3pt] &={\boldsymbol{{X}}}_{e}^{T}\text{(}{{\boldsymbol{{A}}}^{T}}{\boldsymbol{{W}}+\boldsymbol{{WA}}}+\frac{1}{\gamma }\boldsymbol{{WB}}{{{\boldsymbol{{B}}}}^{T}}\boldsymbol{{W}}\text{)}\boldsymbol{{X}}_{e}+\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}} \end{align}

Using above inequality conditions, we may find a small enough positive constant $\mu \gt 0$ , such that

(23) \begin{align} {\dot{\boldsymbol{{V}}}}\text{(}t\text{)}&\le -\mu \boldsymbol{{X}}_{e}^{T}\boldsymbol{{W}}{{\boldsymbol{{X}}}_{e}}+\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}} \nonumber\\[3pt] & \le -\mu \boldsymbol{{V}}(t)+\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}} \end{align}

According to the comparison principle, we obtain

(24) \begin{align} \boldsymbol{{V}}\text{(}t\text{)}& \le{{e}^{-\mu t}}\boldsymbol{{V}}({{t}_{0}})+{{e}^{-ut}}\int _{0}^{t}{\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}}{{e}^{ut}}dt} \nonumber \\[3pt] & ={{e}^{-\mu t}}\boldsymbol{{V}}({{t}_{0}})+{{e}^{-ut}}(\left. \gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}}\frac{1}{\mu }{{e}^{ut}} \right |_{0}^{t}) \nonumber \\[3pt] & ={{e}^{-\mu t}}(\boldsymbol{{V}}({{t}_{0}})-\frac{1}{\mu }\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}})+\frac{1}{\mu }\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}} \nonumber \\[3pt] & \le{{e}^{-\mu t}}({{\lambda }_{\max }}({W})\left \|{{{X}}_{e0}} \right \|-\frac{1}{\mu }\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}})+\frac{1}{\mu }\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}} \end{align}

From the definition of $\boldsymbol{{V}}\text{(}t\text{)}$ , we have

(25) \begin{equation}{{\lambda }_{\min }}(\boldsymbol{{W}})\left \|{{\boldsymbol{{X}}}_{e}} \right \|\le \boldsymbol{{V}}(t)\le{{\lambda }_{\max }}(\boldsymbol{{W}})\left \|{{\boldsymbol{{X}}}_{e}} \right \| \end{equation}

where ${{\lambda }_{\max }}({\boldsymbol{{W}}})$ and ${{\lambda }_{\min }}({\boldsymbol{{W}}})$ are the largest and smallest eigenvalues of $\boldsymbol{{W}}$ , respectively.

(26) \begin{equation} \left \|{\boldsymbol{{X}}}_{e} \right \|\le{e}^{-\mu t}\left(\frac{{\lambda }_{\max }({W})}{{\lambda }_{\min }({W})}\left \|{{{{\boldsymbol{{X}}}}}_{e0}}\|\right .-\frac{\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}} }{\mu{{\lambda }_{\min }}({W})}\right)+\frac{\gamma{{\boldsymbol{{F}}}^{T}}\boldsymbol{{F}}}{\mu{{\lambda }_{\min }}(\boldsymbol{{W}})} \end{equation}

Since the interaction force exerted by humans on the robot cannot be increased infinitely, $\boldsymbol{{F}}$ (t) is bounded. The boundness of $ \|{{\boldsymbol{{X}}}_{e}} \|$ is obtained from the boundedness of $\boldsymbol{{F}}$ (t). Therefore, the stability condition of the system (16) is that there exists positive definite $\boldsymbol{{W}}$ and normal $\gamma$ for given $\boldsymbol{{M}}$ (t), $\boldsymbol{{C}}$ (t), and $\boldsymbol{{K}}({t})$ at any time such that Eq. (19) holds. Therefore, the condition that the system (16) is controllable and stable at every moment is that Eqs. (18) and (19) are both valid.

4.3. Human-like variable admittance parameter regulator

From the analysis in Section 3, it can be seen that the admittance parameters of followers in the human–human collaborative motion change with interaction force and the change rate of interaction force. To make robots collaborate with humans naturally and smoothly in human–robot collaborative motion, the admittance parameter variation rule of the follower is applied to the admittance parameter regulator of the robot. The flow chart of the automatic adjustment of the human-like variable admittance parameters of the robot is shown in Fig. 5. The steps of the human-like variable admittance parameters adjustment are as follows:

• Step 1: Before the human–robot collaborative motion, the robot’s admittance parameters are static when the robot is in standby state.

• Step 2: When the end motion velocity of the robot is $\exists i\in [1,2,3],{{{\dot{\boldsymbol{{X}}}}}_{r}}(i,1)\lt 1\,\text{mm/s}$ , this indicates that the robot is in the motion preparation phase. According to the interaction force, the static admittance parameters were obtained from Eq. (9), and the static admittance parameters were updated according to the real-time change of the interaction force.

• Step 3: When $\exists i\in [1,2,3],{{{\dot{\boldsymbol{{X}}}}}_{r}}(i,1)\ge 1\,\text{mm/s}$ , this indicates that the robot will be in a state of rapid motion. The virtual mass, virtual damping, and virtual stiffness of the robot were updated according to Eqs. (13), (14), and (15), respectively. Otherwise, return to step 2.

• Step 4: To determine whether the updated admittance parameters meet the system stability conditions: $\text{rank}(\boldsymbol{{M}})=3$ and ${{\boldsymbol{{A}}}^{T}}{\boldsymbol{{W}}+\boldsymbol{{WA}}}+\frac{1}{\gamma }\boldsymbol{{WB}}{{{\boldsymbol{{B}}}}^{T}}\boldsymbol{{W}}\lt{0}$ .

• Step 5: When the updated admittance parameters meet the stability conditions, the updated admittance parameters are output to the human-like variable admittance controller. Otherwise, the admittance parameter value of the last time (t − 1) is output to the human-like variable admittance controller.

• Step 6: When $\exists i\in [1,2,3],{{{\dot{\boldsymbol{{X}}}}}_{r}}(i,1)\lt 1\,\text{mm/s}$ indicates that the human–robot collaborative motion is about to end. Otherwise, the admittance parameters are updated through steps 3−5.

Figure 5. The flow chart of automatic adjustment of human-like variable admittance parameter regulator of the robot.

5.1. Experimental setup

A human–robot collaborative motion experimental system was designed to verify the effect of the HVAPR, as shown in Fig. 6. The human–robot collaborative motion experimental system was mainly composed of a epigynous machine (REN9000-28ICO, Lenovo, China), a 7-DOF (degree of freedom) manipulator (robot LBR IIWA, KUKA, Germany), and a 6-dimensional force sensor (HEX-E, Onrobot, Denmark) installed on the robot’s end-effector. The epigynous machine was mainly used to process the data from the robot and the force sensor and output the robot motion command information. MATLAB software was used to write data processing, control, communication programs, etc. The robot was mainly used to execute motion commands from the epigynous machine and collaborate with humans to complete tasks. Java language was used to write robot motion programs. Force sensor was mainly used to measure the interaction force exerted by the human on the robot. The communication mode between the epigynous machine and the force sensor was serial port communication. The communication mode between the epigynous machine and the robot was TCP $\backslash$ IP communication. After 10,000 times of data transfer time statistics, the average data transfer time between the epigynous machine and the robot and the force sensor is 11.20 ms.

Figure 6. Human–robot collaborative point-to-point motion experiment.

The participant worked collaboratively with the robot by using a handle attached to a force sensor. The force sensor collected the interaction force between the human and robot and transmitted the collected force data to the robot. In Fig. 6, building blocks A and B are the starting point and target point of human–robot collaborative point-to-point motion, and the blue cube is the position of the building block. The starting point and target point of the human–robot collaborative motion can be changed by adjusting the position and height of building blocks. The participant (leader) guided the robot (follower) from the starting position along the desired path (the red dotted line in Fig. 6) to the vicinity of the target point through interactive forces.

To assess the performance of the HVAPR, the HVAPR proposed in this paper and VAPR were applied to the human–robot collaborative motion and the performance differences between the two admittance parameter regulators were compared. The VAPR generally adjusts the admittance parameters according to human arm EMG, motion velocity, or interaction force. This paper chose to adjust the admittance parameters according to the human–robot collaborative motion velocity, and the equation is [Reference Kang, Oh, Seo, Kim and Choi17]

(27) \begin{equation} \left [ \begin{matrix}{{c}_{x}} \\ {{c}_{y}} \\ {{c}_{z}} \\ \end{matrix} \right ]=\left [ \begin{matrix}{{{F}_{s}}}/{\left |{\dot{x}} \right |}\; \\ {{{F}_{s}}}/{\left |{\dot{y}} \right |}\; \\ {{{F}_{s}}}/{\left |{\dot{z}} \right |}\; \\ \end{matrix} \right ] \end{equation}

where ${c}_{\text{x}}$ , ${c}_{\text{y}}$ , and ${c}_{\text{z}}$ are the virtual damping of the X-axis, Y-axis, and Z-axis, respectively; $\dot{x}$ , $\dot{y}$ , and $\dot{z}$ are the human–-robot collaborative actual motion velocities of the X-axis, Y-axis, and Z-axis, respectively. According to Eq. (27), when the velocity is constant, the larger the ${F}_{s}$ value, the greater the damping and the greater the interaction force, so the human–robot collaborative motion is laborious. The smaller the ${F}_{s}$ value, the smaller the damping and the interaction force, but there are some problems, such as the absence of large oscillation when the human–robot collaborative motion stops. Therefore, it is necessary to perform the human–robot collaborative motion experiments under the condition of different values of ${F}_{s}$ and select the appropriate ${F}_{s}$ values according to the feeling of the participants (the motion naturalness, smoothness, stability, overall performance, etc.) and the quantitative evaluation of performance (smoothness index and average energy). According to the experimental results, ${F}_{s}=3N$ .

5.2. Experiment task

To compare the performance differences of the two admittance parameter regulators, the human–robot collaborative performance based on HVAPR and VAPR was assessed from the perspective of humans and robots, respectively. This paper used a questionnaire to investigate participants’ feelings regarding the human–robot collaborative motion to assess and compare the human–robot collaborative performance qualitatively from the humans’ perspective. The smoothness of the robot’s trajectory, average motion energy, and interaction force was analyzed, and the human–robot collaborative performance was quantitatively assessed and compared from the robot’s perspective.

In the human–robot collaborative motion experiment, after completing a point-to-point motion, each participant answered the questions in Table V according to the motion’s naturalness, smoothness, stability, overall performance, etc. [Reference Roveda, Testa, Shahid, Braghin and Piga12]. Questions 1–4 and 6 adopt 5-point Likert scale [Reference Schmidtler, Bengler, Dimeas and Campeau-Lecours36].

Table V. Questionnaire question table.

The possible responses of the participants to the different questions are shown in Table VI. The semantic scale used for answering the six above questions was treated as a continuous scale since each interval of the scale was of equal proportion. Therefore, the data collected from the first, second, third, fourth, and sixth questions were treated as continuous measures [Reference Sheikholeslami, Moon and Croft37].

Table VI. Table of responses of participants to different questions.

Answers to the first question reflect how satisfied participants are with each human–robot collaborative motion. Questions 2–4 are a validation of the different dimensions of the answer to question 1. Answers to questions 1–4 showed a high level of internal consistency (each participant’s question answered with a Cronbach’s $\alpha$ greater than 0.9).

Answers to the fifth question reflect whether participants can correctly distinguish between the two admittance parameter regulators based on their feelings concerning the human–robot collaborative motion. The high recognition ratio of admittance parameter regulators reflects that the two admittance parameter regulators give different feelings to participants. Answers to the sixth question serve as a confidence measure (recognition confidence) of responses to the fifth question, that is, the participants’ confidence in recognizing the robot admittance parameter regulator. The higher the recognition confidence value is, the more confident the participants are in the admittance parameter regulator recognized by themselves. The answers to questions 5 and 6 show the difference in the influence of these two admittance parameter regulators on the human–robot collaborative motion.

Before the human–robot collaborative motion experiment, each participant underwent a preliminary experiment to assess the difference between the two admittance parameter regulators. In the human–robot collaborative motion experiment, the robot randomly selected one of the two admittance parameter regulators to adjust the admittance parameters. The participants did not know which admittance parameter regulator the robot adopted in each human–robot collaborative motion. To prevent participants from excessively adapting to the robot motion state [Reference Burdet, Osu, Franklin, Milner and Kawato38], participants’ correct assessment of human–robot collaborative performance was affected. After completing a point-to-point motion task, the robot’s current position was used as the starting point of the next point-to-point motion task. The new target point was determined by randomly changing the placement position and height of the building blocks, as shown in Fig. 6.

After each point-to-point motion experiment, each participant answered all the questions within 10 s. Six participants participated in the human–robot collaborative motion experiment. The participants (all males, with mean $age=27.5\pm 2.5$ years, ${{height}}=174\pm 6$ mm, and $weight=68.2\pm 4.8$ kg) are right handed and have no physical problems. Each participant performed 600 times of the human–robot collaborative point-to-point motion experiment, and the program was used to record the data of the robot’s motion and interaction force. After each participant completed all the experiments, the overall satisfaction of the human–robot collaborative motion experiment was assessed. Two volunteers counted all participants’ answers to the questionnaire to ensure the accuracy of the data. The results showed that the questionnaire data of six participants had a high level of internal consistency (Cronbach’s $\alpha =0.915$ ).

This paper selected the smoothness J of the motion trajectory and the average energy ${P}_{\text{ave}}$ as the motion performance assessment indexes to quantitatively assess the human–robot collaborative performance. The smaller the J value, the better the motion stability. The lower the ${P}_{\text{ave}}$ value, the lower the energy consumption, and the easier the motion. The smoothness index J and the average energy ${P}_{\text{ave}}$ are, respectively, defined as [Reference Kang, Oh, Seo, Kim and Choi17, Reference Hamad, Aydin and Basdogan24]:

(28) \begin{equation} J=\frac{1}{{{t}_{f}}-{{t}_{0}}}\int _{{{t}_{0}}}^{{{t}_{f}}}{\dddot{x}_{r}^{2}(t)dt} \end{equation}

(29) \begin{equation}{{P}_{\text{ave}}}=\frac{1}{{{t}_{f}}-{{t}_{0}}}\int _{{{t}_{0}}}^{{{t}_{f}}}{\left |{{f}_{x}}(t){{{\dot{x}}}_{r}}(t) \right |dt} \end{equation}

where ${t}_{0}$ and ${t}_{f}$ are the start time and end time, respectively.

5.3. Analysis of questionnaire survey results

The analysis of the questionnaire data shows that the answers to questions 1–4 have a high level of internal consistency. Questions 2–4 are responses to different dimensions of question 1. Therefore, only the results of the answer to question 1 were shown. The satisfaction assessment results of all the participants on the two admittance parameter regulators were obtained by statistical analysis of the answers to question 1, as shown in Table VII. As shown in Table VII, the participants’ satisfaction with HVAPR is mainly concentrated in the “general,” “satisfied,” and “very satisfied” categories, while satisfaction with VAPR is mainly concentrated in the “unsatisfied,” “general,” and “satisfied” categories. Overall, participants are more satisfied with HVAPR than with VAPR. HVAPR, which adjusts the admittance parameters based on the change rate of interaction force, enables the robot to provide admittance parameters with a good degree of matching to the current motion of the participant, thus resulting in better satisfaction in the human–robot collaborative motion. Since HVAPR is extracted from the variation rules of the admittance parameters in the human–human collaborative motion, the admittance parameters obtained cannot fully match the motion of each participant, resulting in “unsatisfied” and “general” assessments.

Table VII. Satisfaction assessment results of the two admittance parameter regulators.

The admittance parameters obtained by VAPR were poorly matched with participants’ current motion state, resulting in poor overall satisfaction with the human–robot collaborative motion. According to the participants’ assessment of the human–robot collaborative motion, when the admittance parameter regulator was VAPR, if participants actively adjusted their motion state to adapt to the admittance parameter regulation method, then the satisfaction regarding the human–robot collaborative motion was better. Therefore, there are “satisfied” and “very satisfied” assessments in the human–robot collaborative satisfaction survey. However, requiring participants to actively adjust their motion to adapt to VAPR is more difficult.

The recognition ratio of participants regarding the two admittance parameter regulators was obtained by statistical analysis of the answers to question 5, as shown in Table VIII. The ratio of the different recognition confidence of the two admittance parameter regulators was obtained by statistical analysis of the answers to question 6, as shown in Table IX. Table VIII shows that the recognition ratio of each participant for the two admittance parameter regulators was between 70% and 90%. The recognition confidence in recognizing the two admittance parameter regulators was mainly concentrated in the “certain” and “very certain” categories, as shown in Table IX. This indicates that the two admittance parameter regulators gave participants different feelings, so the participants could confidently recognize which admittance parameter regulator the robot used. Although the recognition ratio of the admittance parameter regulator is relatively high, it is not close to 100%. The “uncertain” and “general” categories’ ratio in the recognition confidence survey is about 15%. This indicates that the two admittance parameter regulators gave participants similar feelings in some cases, so the robot admittance parameter regulator could not be correctly recognized.

Table VIII. Recognition ratio of admittance parameter regulators for all participants.

Table IX. Recognition confidence ratio of admittance parameter regulators.

5.4. Analysis and assessment of human–robot collaborative motion performance

To evaluate the performance of the proposed HVAPR, the motion information of the X-axis of the robot is taken as an example. The change curves of the robot’s end path ( ${x}_{r}$ ), the interaction force ( ${f}_{x}$ ), and the actual velocity ( ${\dot{x}}_{r}$ ) are shown in Fig. 7. For clarity and brevity, we show only representative experimental results in Fig. 7. The oscillation ratio of the path near the target point and the RMSE between the actual velocity and the minimum jerk model velocity (it is ideal velocity model of human arm terminal motion) were calculated, and the results are shown in Table X.

Table X. Oscillation ratio, overshoot amount, and RMSE of human–robot collaborative motion.

When compared with the VAPR, the ${f}_{x}$ change curve of interaction force in Fig. 7 demonstrates that the human–robot collaborative interaction force based on the HVAPR is smaller and changes stably. According to the change curve of robot motion trajectory in Fig. 7(a), the robot motion trajectory appears to obviously oscillate and overshoot when it is near the target point. The oscillation ratio of the robot trajectory (the ratio of the experimental number of oscillations on the trajectory near the target point to the total number of experiments) is 70.19%, and the average overshoot amount [the distance between the maximum deviation value and the target point (Figs. 7a)] is 0.0589 m, as shown in Table X. This led participants to view the robot as dangerous and thus to provide a “less satisfied” rating. This is due to the poor degree of matching between the admittance parameters regulated by the VAPR and the participant’s (leader) motion state. As shown in Fig. 7(a), the robot reached the target point in 4.5 s, but the motion velocity remained high. According to Eq. (27), the damping was small at this time, and the robot could not stop moving in time, resulting in the robot exceeding the target point (with a large overshoot). To return to the vicinity of the target point, at 5 s, the participants applied a large reverse interaction force to reverse the motion of the robot, resulting in an oscillating trajectory of the robot. Although the human–robot collaborative motion trajectory based on a human-like admittance parameter regulator also appears to oscillate, the oscillation ratio is only 31.33%, and the average overshoot is 0.0275 m, as shown in Table X. This effectively improved the participants’ satisfaction with the human–robot collaborative motion.

It can also be seen from Fig. 7 that the human–robot collaborative motion velocity curves based on these two admittance parameter regulators show “bell-shaped” distribution rule, which conforms to the characteristics of the minimum jerk model [Reference Jing, Chengyun and Biyun39]. Therefore, the minimum jerk model was selected to represent the ideal motion velocity curve. As shown in Table X, the human–robot collaborative motion velocity of VAPR is significantly different from the minimum jerk model velocity, and the RMSE is 0.1100. The human–robot collaborative motion velocity of HVAPR has a small difference from the minimum jerk model velocity and has a smaller RMSE (0.0582).

In conclusion, the velocity curve of the human–robot collaborative motion based on the HVAPR was closer to that of the minimum jerk model than that of the VAPR. The HVAPR had a lower trajectory oscillation ratio (31.33%), smaller average overshoot (0.0275 m), and smaller interaction force.

The motion performance of human–robot collaboration based on these two admittance parameter regulators was quantitatively assessed using the aforementioned two motion performance assessment indexes (Eqs. 28 and 29). The mean value and standard deviation of n (n = 3600) times the human–robot collaborative motion assessment indexes were calculated, and the results are shown in Table XI.

Table XI. Average energy and smoothness of the human–robot collaborative motion.

Figure 7. Motion and interaction forces at the end of robots.

As seen in Table XI, HVAPR and VAPR have significant differences in average energy ( $F(1,10)=52.761\gt 4.965,p\lt 0.05$ ) and motion smoothness ( $F(1,10)=22.039\gt 4.965,p\lt 0.05$ ). HVAPR is significantly superior to VAPR in average energy and smoothness, indicating that HVAPR can make the robot move more smoothly and consume less energy in the human–robot collaborative motion. Meanwhile, HVAPR has some advantages in terms of standard deviation (0.67 $E^{-5}$ ) of motion trajectory smoothness and standard deviation (0.86) of average energy, indicating better repeatability of HVAPR. Combined with the results of the previous questionnaire, it can be seen that in the human–robot collaborative motion, HVAPR can not only give participants a better collaborative experience but also make the robot have better motion performance.