Introduction

In recent years, robots are expected to reduce human workload and ease the threats of fatigue working. Since the accuracy and repeatability of the robotic tasks, most of operation modes are implemented in predefined tasks without human partners. Furthermore, robots are inadequate for the tasks that require quick judgment, flexibility, and adaptability under the unstructured environment; thus, humans need to complement the properties that robots lack, as well as these tasks are named physical human–robot interaction (pHRI) [1]. However, conventional robot is difficult to apply to these complex fields under the tasks of pHRI, such as space exploration [2], disaster relief [3] and medical assistance [4]. Modular robot manipulator (MRM) is a robot with multiple modules, which can be added or subtracted modules according to the requirements or variable tasks. In addition, these modules have standard electrical and mechanical interfaces, each of which integrates communication unit, drive unit, control unit and sensing unit [5]. According to the above advantages, MRMs are applicable to pHRI in complex fields because of the properties that are oriented to varied tasks and different human partners [6]. The core of pHRI is human motion intention identification, the robot can naturally respond to human intention to ensure the safety of human partners. Moreover, the robots should communicate the current state of human partner and predict human motion intention in pHRI [7].

In fact, the signals can be divided into biological signals and robot signals, which are utilized for intention prediction. The human body generates a variety of biological signals, such as surface electromyography [8], electroencephalogram [9] and electromyography (EMG) [10]. Recently, these signals have been widely utilized to predict the human motion intention in pHRI [11]. Many scholars have carried out a lot of researches on the human motion intention estimation using EMG signals. In the application of pHRI systems, EMG signals are mainly utilized to predict human motion intention. The typical continuous motion estimation research [12] proposed a control framework to obtain the robot motion intention. Human motion intention prediction approach using EMG consist in model-based approaches and model-free approaches. The biomechanical model method is meaningful to reveal the mechanism of movement from the microscopic perspective, namely the point of view of muscle tissue. In the 1940s, Hill first proposed the classical muscle mechanics model [13]. Subsequently, several studies concentrated on the human motion intention prediction based on Hill muscle model [14, 15]. Meanwhile, the method of model-free approach for intention recognition of continuous joint motion has been developed rapidly, such as nonlinear autoregressive model [16], and high-order polynomial modeling [17].

Although EMG signals are generally used to predict human motion intention, they are easily affected by electrode position, adipose tissue thickness, body temperature and perspiration. In addition, EMG signal has a large amount of information, which requires a complex pre-processing process to extract the signal of control input [18]. Therefore, many researchers proposed a number of human motion intention estimation-based pHRI control methods via the sensor information mounted on the robot, e.g., force sensors, joint torque sensors, position sensors, etc., are used in pHRI tasks. The goal of pHRI actually acts on the robot, and the interaction performance is reflected by these sensors that are easily obtained which mounted on the robot. Wakita et al. [19] designed an intelligent walking cane robot and proposed ”Intentional Direction” to describe the motion intention. Erden et al. [20] proposed the law of zero momentum conservation-based sensorless interactive control strategy, in which the robot received the estimated human intention and updated its reference position. Under human–robot cooperative tasks, the problem of continuous human motion estimation can be dealt with by modified least squares algorithm [21], growing hidden Markov models [22]. Rahul et al. [23] proposed an iterative learning control to predict the human motion intention based on the inversion method. Furthermore, many control methods are concentrated on the human motion intention recognition under the rapid point-to-point motion process, such as, admittance control [24], learning control [25]. Meanwhile, Li et al. [25] presented the neural network control theory to deal with the arbitrary continuous motion trajectory. On this basis, perceptual enhancement technique [26] and joint space observer-based admittance adaptive control scheme [27] are utilized to improve the tracking performance of the robot and the performance at the interaction point, respectively. Li et al. [28] proposed a Bayesian estimation method to predict the impedance and the human motion intention of pHRI tasks and an adaptive impedance control method is utilized to track the target impedance model. However, the interaction force information of the existing human motion intention estimation method is obtained from the six-dimensional force/torque sensor or joint torque sensor. Indeed, the installation of joint torque sensor and six-dimensional force sensor will increase the cost burden, and furthermore, the installation of these sensors will lead to more measurement errors and calibration errors, which will cause the decline of the overall accuracy of the MRM. Therefore, it is necessary to design a novel torque-sensorless human motion intention estimation method for pHRI tasks, which can eliminate the necessity of force/torque sensors.

To achieve natural physical interaction, the performance of the task is not limited to fast and accurate trajectory tracking performance, but also includes security issues during interaction. One way to ensure the safety is to maintain the interaction force within an acceptable scale. To create the ’safety’ interaction environment, the approaches are divided into passive method and active method, which are utilized to reduce the interaction force for pHRI tasks. The passive method can reduce the interaction force and increase the compliance of robot, which includes the mechanical design of the robot. A changeable stiffness mechanism structure is proposed, which installs on the robot joint [29]. The mechanism maintains high joint stiffness in classical work state of the robot, and reduces joint stiffness dramatically in the process of pHRI, thus as to realize the compliance property of the robot. In addition, the serial elastic actuation design is represented to improve the compliance property of the robot [30]. However, this kind of design will increase the complexity of the robot system and result a large mount of unimaginable difficulties. Another method to improve the safety of the interaction environment is active control mode, which includes switch control strategy based on residual momentum observer [31], collision detection-based compliance control method [32]. However, the above-mentioned works mainly focus on human-robot collision detection, which are difficult to use or limited in practice. In the context of pHRI, the appropriate control method is utilized to pledge the stability and security of the interaction process. For the past few years, many scholars have conducted research on the control problem under the background of pHRI, such as impedance control [33], force/position control [34], human–robot cooperative control [35, 36]. However, most of the existing control methods under pHRI have not considered the scheme of adopting human motion intention estimation, and these methods either assist human motion with a given curvature, or estimate the expected trajectory of the robot with accurate dynamic model information. Indeed, since humans always handle the movement of end-effectors, and the robots maintain the trajectory tracking, so that the human is not really free to manipulate. Estimating the motion intention of the human using the limb model can reflect the dominant operation of the human during pHRI. Therefore, it is an urgent requirement to investigate an approximate interaction control method that realizing completely free operation under human motion intention estimation with high-performance position tracking in pHRI.

Fig. 1
figure 1

Harmonic drive compliance model

In this paper, decentralized robust interaction control method of pHRI-oriented MRM system is presented based on harmonic drive compliance model. First, dynamic model of MRMs is developed based on the harmonic drive compliance model, and the state space of is given. Second, a novel torque-sensorless human motion intention is estimated based on the harmonic drive compliance model by only using the information of local dynamic position measurements. In addition, decentralized robust interaction control strategy is derived by using the robot dynamics model, which maintains interaction contact force in a safe range on the basis of high-precision trajectory tracking. Based on the Lyapunov theory, the uniformly ultimately boundedness (UUB) of the trajectory position error and the contact torque error are proved. Finally, pHRI experiments confirm the effectiveness and superiority of the proposed method.

The main contributions of the paper can be summarized as follows:

  1. 1.

    To the best of our knowledge, it is the first time to propose a novel torque-sensorless human motion intention estimation method based on the harmonic drive compliance model to solve the pHRI problem of MRM, which will recognize the human motion intention by only using the information of local dynamic position measurements of each joint module.

  2. 2.

    This paper developed a novel decentralized robust interaction control method based on the human motion intention estimation, which maintained the interaction contact force within the scope of the security and realized the high performance of tracking trajectory by considering the model uncertainties of MRM. The interaction control method is proposed based on the safety contact force, which can establish the ‘safety’ interaction environment. In addition, the uncertainty decomposition-based robust control can compensate the dynamic uncertainties of MRM. The superiority of the proposed method is confirmed using the established experimental platform.

The construction of this paper is described as follows. In Sect. “Dynamic Model Formulation”, joint torque estimation method based on the harmonic drive compliance model is firstly presented, and the dynamics model of the modular robot manipulators under the physical human–robot interaction environment is introduced. In Sect. “Human Motion Intention Estimation-based Decentralized Robust Interaction Control”, the human motion intention is identified via harmonic drive compliance model when human expect to tow the robot. Meanwhile, the decentralized robust interaction controllers of each joint module are presented. In Sect. “Experiments”, the effectiveness and advantages of the proposed method are verified by experiments. Finally, Sect. “Conclusion” draws the conclusion.

Dynamic model formulation

Harmonic drive compliance model

MRM’s modulars integrate a rotary joint with a harmonic drive (HD). It is presented by [37] that HD is a mechanical structure composed of a wave generator, a flexspline and a circular spline, which has the ability of kinematic and torque transmission. The ideal kinematic and torque transport relationship of HD is formulated as:

$$\begin{aligned} {\theta _{wi}} = {\gamma _i}{\theta _{fi}}, \end{aligned}$$
(1)
$$\begin{aligned} {\tau _{wi}} = {\tau _{fi}}/{\gamma _i}, \end{aligned}$$
(2)

where i represents the ith module, \({\gamma _i}\) denotes reduction ratio, \({\theta _{wi}}\),\({\tau _{wi}}\) reflects reflect the position and torque of wave generator, \({\theta _{fi}}\) and \({\tau _{fi}}\) represent the same meaning of flexspline.

Since the circular spline is usually fixed, so that \({\theta _{ci}} = 0\), and the circular spline torque is ignored. [38] obviously presents that the kinematic and torque transport relationship are nonlinear, the harmonic drive compliance model is proposed by considering motor friction, compliance of HD and kinematic error as shown in Fig. 1.

By considering motor friction in the transport process, we can define the torque of wave generator in (3),

$$\begin{aligned} {\tau _{wi}} = {\tau _i} - {\tau _{mri}} = \left( {{\tau _{fi}} - {\tau _{fri}}} \right) /{\gamma _i}, \end{aligned}$$
(3)

where \({\tau _i}\) is the motor torque. Meanwhile, the torsional angle of flexspline and wave generator are calculated as follows:

$$\begin{aligned} \Delta {\theta _{wi}} = {\theta _{woi}} - {\theta _{wli}},\Delta {\theta _{fi}} = {\theta _{foi}} - {\theta _{fli}}, \end{aligned}$$
(4)

where \(\theta _{woi}\),\(\theta _{wli}\) reflect the angle of wave generator at center side and motor side, \(\theta _{foi}\),\(\theta _{fli}\) reflect the angle of the flexspline at gear-part and load-part, \({\theta _{foi}}\) and \({\theta _{wli}}\) are the information which are measured by an absolute encoder mounted at link side and an incremental encoder mounted at motor side. The torsional angle of HD is formulated by utilizing these position angle as:

$$\begin{aligned} \Delta {\theta _i} = {\theta _{foi}} - \left( {{\theta _{wli}}/{\gamma _i}} \right) . \end{aligned}$$
(5)

By rewriting (5), one obtains,

$$\begin{aligned} \begin{array}{l} \Delta {\theta _i} = \left( {\frac{{{\theta _{woi}}}}{{{\gamma _i}}} - \frac{{{\theta _{wli}}}}{{{\gamma _i}}}} \right) + {\theta _{foi}} - {\theta _{fli}} + \left( {{\theta _{fli}} - \frac{{{\theta _{woi}}}}{{{\gamma _i}}}} \right) \\ {} \mathrm{{ }} = \Delta {\theta _{fi}} + \left( {\Delta {\theta _{wi}}/{\gamma _i}} \right) + \Delta {\theta _{upi}}\\ {}\mathrm{{ }} \approx \Delta {\theta _{fi}} + \left( {\Delta {\theta _{wi}}/{\gamma _i}} \right) , \end{array} \end{aligned}$$
(6)

where \(\Delta {\theta _{upi}}\) is the kinematic error. The kinematic error is infinite which approximately equals to zero.

To best of our knowledge, HD has the hysteresis and stiffness properties. In addition, these properties will lead local elastic coefficient aggrandizes with the increasing of the flexspline torque, and the elastic coefficient of flexible \(k_{fi}\) is formulated as follows:

$$\begin{aligned} {k_{fi}}= & {} d{\tau _{fi}}/d\Delta {\theta _{fi}}, \end{aligned}$$
(7)
$$\begin{aligned} {k_{fi}}= & {} {k_{fi0}}\left( {{{\left( {{c_{fi}}{\tau _{fi}}} \right) }^2} + 1} \right) , \end{aligned}$$
(8)

where \({k_{fi0}},{c_{fi}}\) are constants. If \({k_{fi0}} \ne 0\), flexspline torsion can be defined as:

$$\begin{aligned} \Delta {\theta _{fi}} = \int _0^{{\tau _{fi}}} {\frac{{d{\tau _{fi}}}}{{{k_{fi}}}}} = \frac{{\arctan \left( {{c_{fi}}{\tau _{fi}}} \right) }}{{{c_{fi}}{k_{fi0}}}}. \end{aligned}$$
(9)

Based on the purpose of replicating the hysteresis shape of this stiffness property, the local elastic coefficient of the wave generator is formulated as:

$$\begin{aligned} {k_{wi}} = {k_{wi0}}{e^{{c_{wi}}\left| {{\tau _{wi}}} \right| }}, \end{aligned}$$
(10)

where \({k_{wi0}},{c_{wi}}\) reflect constants. If \({k_{wi0}} \ne 0\), local elastic coefficient of wave generator is formulated as:

$$\begin{aligned} \Delta {\theta _{wi}} = \int _0^{{\tau _{wi}}} {\frac{{d{\tau _{wi}}}}{{{c_{wi}}}}} = \frac{{{\mathop {\mathrm{sgn}}} \left( {{\tau _{wi}}} \right) }}{{{c_{fi}}{k_{wi0}}}}\left( {1 - {e^{{c_{wi}}\left| {{\tau _{wi}}} \right| }}} \right) . \end{aligned}$$
(11)

The total torsion angle of HD is formulated by substituting (9), (11) into (6),

$$\begin{aligned} \begin{array}{lllll} \Delta {\theta _i} &{}= \Delta {\theta _{wi}}/{\gamma _i} + \Delta {\theta _{fi}} = \frac{{{\mathop {\mathrm{sgn}}} \left( {{\tau _{wi}}} \right) }}{{{c_{fi}}{k_{wi0}}}}\left( {1 - {e^{{c_{wi}}\left| {{\tau _{wi}}} \right| }}} \right) \\ &{}\quad + \frac{{\arctan \left( {{c_{fi}}{\tau _{fi}}} \right) }}{{{c_{fi}}{k_{fi0}}}}. \end{array} \end{aligned}$$
(12)

Deforming (12), the joint torque is estimated as (13) which equals to flexspline output torque of HD.

$$\begin{aligned}&{\tau _{fi}} = \tan \left( {{c_{fi}}{k_{fi0}}\left( {\Delta {\theta _i} - \frac{{{\mathop {\mathrm{sgn}}} \left( {{\tau _{wi}}} \right) }}{{{\gamma _i}{c_{fi}}{k_{wi0}}}}\left( {1 - {e^{{c_{wi}}\left| {{\tau _{wi}}} \right| }}} \right) } \right) } \right) \nonumber \\&/{c_{fi}} = {\tau _{fci}}, \end{aligned}$$
(13)

where \({\tau _{fci}}\) is determined part of estimated joint torque, and the dynamic model of MRM is developed based on harmonic drive compliance model in the next part.

Dynamic model formulation

An MRM system is developed under pHRI, which consists a human limb and a robot manipulator, and the motion of the MRM is constrained by human.

According to the robot kinematics, one obtains,

$$\begin{aligned} z\left( t \right)= & {} \xi \left( q \right) , \end{aligned}$$
(14)
$$\begin{aligned} \dot{z}\left( t \right)= & {} J\left( q \right) \dot{q}, \end{aligned}$$
(15)
$$\begin{aligned} \ddot{z}\left( t \right)= & {} \dot{J}\left( q \right) \dot{q} + J\left( q \right) \ddot{q}, \end{aligned}$$
(16)

where \(z\left( t \right) \in {R^n}, q \in {R^n}\) represent Cartesian space position as well as joint space position, \(\xi \left( \cdot \right) \) reflects the mapping matrix of joint space to Cartesian space and \(J\left( q \right) \in {R^{n \times n}}\) represents Jacobian matrix.

According to the structure of the robot module, because of the angular position of flexspline is the displacement of the linkside of the robot as shown in Fig. 1, we can get the relation \({\theta _i} = {\theta _{foi}}\), where \({\theta _i}\) is the ith joint position. Meanwhile, combine joint torque estimation which is given in (13) and HD-based dynamic model of MRM, the dynamic model of ith joint subsystem under physical human–robot interaction is formulated as follows:

$$\begin{aligned} {I_{mi}}{\gamma _i}{{\ddot{q}}_i}+ & {} {f_i}\left( {{q_i},{{\dot{q}}_i}} \right) + {I_{mi}}\sum \limits _{j = 1}^{i - 1} {o_{mi}^T{o_j}{{\ddot{q}}_j}} + \frac{{{\tau _{fci}}}}{{{\gamma _i}}}\nonumber \\&\quad +{I_{mi}}\sum \limits _{j = 2}^{i - 1} {\sum \limits _{k = 1}^{j - 1} {o_{mi}^T\left( {{o_k} \times {o_j}} \right) {{\dot{q}}_k}{{\dot{q}}_j}} } = {\tau _i} + {\tau _{iext}},\nonumber \\ \end{aligned}$$
(17)

where \({I_{mi}}\) reflects moment of inertia, \({\gamma _i}\) represents HD reduction ratio \(\left( {{\gamma _i} \ge 1} \right) \), \(q_i\) reflects the ith joint position, \({f_i}\left( {{q_i},{{\dot{q}}_i}} \right) \) reflect ith joint friction part \({I_{mi}}\sum \nolimits _{j = 1}^{i - 1} {o_{mi}^T{o_j}{{\ddot{q}}_j}} + {I_{mi}}\sum \nolimits _{j = 2}^{i - 1} {\sum \nolimits _{k = 1}^{j - 1} {o_{mi}^T\left( {{o_k} \times {o_j}} \right) {{\dot{q}}_k}{{\dot{q}}_j}} }\) reflect couplings terms, \({o_{mi}}\) denote the unit vectors of ith rotor axis, \({o_{i}}\) represents the unity vectors of ith joint axis, \({\tau _{fci}}\) is joint torque estimation which is given in (13), \({\tau _i}\) represents the control input, \({\tau _{iext}} = {\tau _{fci}} - {\tau _{dbi}},{\tau _{ext}} = {[{\tau _{1ext}}, \cdots {\tau _{iext}}, \cdots {\tau _{next}}]^T} \in {R^n}\) where \({\tau _{iext}}\) denotes the interaction torque, \({\tau _{db}}\) is the torque of the offline database (see in Remark 1). The interaction force \({f _{ext}} = {[{f _{ext1}}, \cdots {f _{exti}}, \cdots {f _{ext}}]^T} \in {R^n}\) on the end-effector of MRM is obtained by transforming the joint interaction torque \({\tau _{ext}}\) as follows:

$$\begin{aligned} {f_{ext}} = {J^{{-T}}}\left( q \right) {\tau _{ext}}, \end{aligned}$$
(18)

where \({f _{ext}} = {[{f _{ext1}}, \cdots {f _{exti}}, \cdots {f _{extn}}]^T} \in {R^n}\),\(f _{exti}\) is ith scalar of the vector \(f _{ext}\), which is the counterpart of ith interaction torque \({\tau _{ext}}\). The interaction contact force \(f _{exti}\) is utilized to estimate human motion intention on the third section, which is formulated based on harmonic drive compliance model.

Remark 1

An offline database is developed which reflects the relation between the joint torque \({\tau _{db}}\mathrm{{ = [}}{\tau _{db1,}} \cdots {\tau _{dbi,}} \cdots {\tau _{dbn}}{\mathrm{{]}}^{\mathrm{T}}} \in {R^n}\) and the corresponding joint position q in conventional position control without human. Furthermore, the harmonic drive compliance model-based joint torque estimation \({\tau _{fci}}\) is constituted by \({\tau _{db}}\) and \({\tau _{ext}}\) in pHRI. Therefore, the joint interaction torque \({\tau _{ext}}\) can be obtained by the joint torque estimation method and the offline database.

According to [39], we can obtain linearization joint friction \({f_i}\left( {{q_i},{{\dot{q}}_i}} \right) \) in (17),

$$\begin{aligned} {f_i}({q_i},{\dot{q}_i})\approx & {} {{\hat{b}}_i}{\dot{q}_i} + \left( {{{{\hat{f}}}_{ci}} + {{{\hat{f}}}_{si}}{e^{\left( { - {{{\hat{f}}}_{\tau i}}\dot{q}_i^2} \right) }}} \right) {\mathop {\mathrm{sgn}}} \left( {{{\dot{q}}_i}} \right) \nonumber \\&+ {f_{qi}}\left( {{q_i},{{\dot{q}}_i}} \right) + Y\left( {{{\dot{q}}_i}} \right) {{\tilde{F}}_i}, \end{aligned}$$
(19)

where \({f_{ci}}\),\({f_{si}}\),\({f_{\tau i}}\),\({b_i}\) denote the linearization joint friction-related parameter, which are unknown and nonconstant. The nominal values of \({b_i},{f_{ci}},{f_{si}}\) and \({f_{\tau i}}\) are reflected as \({{\hat{b}}_i},{{\hat{f}}_{ci}},{{\hat{f}}_{si}}\) and \({{\hat{f}}_{\tau i}}\) respectively, \({f_{qi}}({q_i},{\dot{q}_i})\) denotes the friction modeling errors, \({{\tilde{F}}_i} = \left[ - {{\hat{b}}_{fi}} + {b_{fi}}, - {{\hat{f}}_{ci}} + {f_{ci}}, - {{\hat{f}}_{si}}\right. \)\(\left. + {f_{si}}, - {{\hat{f}}_{\tau i}} + {f_{\tau i}}\right] ^{\mathrm{T}}\) represents friction parametric uncertainty, \(Y\left( {{{\dot{q}}_i}} \right) \) formulates as follows:

$$\begin{aligned} Y\left( {{{\dot{q}}_i}} \right)= & {} [{\dot{q}_i},{\mathop {\mathrm{sgn}}} \left( {{{\dot{q}}_i}{} } \right) ,{} {e^{\left( { - {{{\hat{f}}}_{\tau i}}\dot{q}_i^2} \right) }}{\mathop {\mathrm{sgn}}} \left( {{{\dot{q}}_i}} \right) , \nonumber \\&- {{\hat{f}}_{si}}\dot{q}_i^2{e^{\left( { - {{{\hat{f}}}_{\tau i}}\dot{q}_i^2} \right) }}{} {\mathop {\mathrm{sgn}}} \left( {{{\dot{q}}_i}{} } \right) ].{} {} \end{aligned}$$
(20)

Rewrite the terms \({I_{mi}}\sum \nolimits _{j = 1}^{i - 1} {o_{mi}^{\mathrm{T}}{o_j}{{\ddot{q}}_j}} \) and \({I_{mi}}\sum \nolimits _{j = 2}^{i - 1} \sum \nolimits _{k = 1}^{j - 1} o_{mi}^{\mathrm{T}}({o_k} \times {o_j}){{\dot{q}}_k}{{\dot{q}}_j} \) for facilitating as follows:

$$\begin{aligned}&\begin{array}{l} {I_{mi}}\sum \limits _{j = 1}^{i - 1} {o_{mi}^{\mathrm{T}}{o_{qj}}{{\ddot{q}}_{qj}}} = {I_{mi}}\sum \limits _{j = 1}^{i - 1} {D_j^i{{\ddot{q}}_j} } \\ \quad =\sum \limits _{j = 1}^{i - 1} {[{I_{mi}}{\hat{D}}_j^i{} {} {} {} {} {} {} {I_{mi}}]{{[{{\ddot{q}}_j}{} {} {} {} {} {} {} {\tilde{D}}_j^i{{\ddot{q}}_j}]}^{\mathrm{T}}}} = \sum \limits _{j = 1}^{i - 1} {I_j^iU_j^i}, \end{array} \end{aligned}$$
(21)
$$\begin{aligned}&\begin{array}{l} {I_{mi}}\sum \limits _{j = 2}^{i{} - 1} {\sum \limits _{k = 1}^{j{} - 1} {o_{mi}^{\mathrm{T}}} } ({o_{qk}} \times {o_{qj}}){{\dot{q}}_k}{{\dot{q}}_j} = {I_{mi}}\sum \limits _{j = 2}^{i{} - 1} {\sum \limits _{k = 1}^{j{} - 1} {\Theta _{kj}^i{{\dot{q}}_k}{{\dot{q}}_j}} } \\ \quad = \sum \limits _{j = 2}^{i{} - 1} {\sum \limits _{k = 1}^{j{} - 1} {[{I_{mi}}{\hat{\Theta }} _{kj}^i{} {} {} {} {} {} {I_{mi}}]{{[{{\dot{q}}_k}{{\dot{q}}_j}{} {} {} {} {} {} {\tilde{\Theta }} _{kj}^i{{\dot{q}}_k}{{\dot{q}}_j}]}^{\mathrm{T}}}} } \\ \quad = \sum \limits _{j = 2}^{i{} - 1} {\sum \limits _{k = 1}^{j{} - 1} {J_{kj}^iV_{kj}^i} }, \end{array} \end{aligned}$$
(22)
$$\begin{aligned}&\begin{array}{l} {U_{zi}} = \sum \limits _{j = 1}^{i - 1} {I_j^iU_j^i} \\ {V_{zi}} = \sum \limits _{j = 2}^{i - 1} {\sum \limits _{k = 1}^{j - 1} {J_{kj}^iV_{kj}^i} } . \end{array} \end{aligned}$$
(23)

We obtain relationships between \(D_j^i = o_{mi}^{\mathrm{T}}{o_{qj}}\) and \({\hat{D}}_j^i = D_j^i - {\tilde{D}}_j^i\) in (21) and (22), where \(D_j^i\) is dot product of \({o_{mi}}\), \({o_{qj}}\), and \({\tilde{D}}_j^i\) denotes the alignment error. We obtain \(\Theta _{kj}^i = o_{mi}^{\mathrm{T}}({o_{qk}} \times {o_{qj}})\), \({\hat{\Theta }} _{kj}^i = \Theta _{kj}^i - {\tilde{\Theta }} _{kj}^i\), where \(\Theta _{kj}^i\) is the dot product of \({o_{mi}}\), \({o_{qk}} \times {o_{qj}}\), \({\tilde{\Theta }} _{kj}^i\) reflects the same meaning as \({\tilde{D}}_j^i\). In addition, \(U_j^i\) and \(V_{kj}^i\) are considered as model uncertainties, and the previous model uncertainties have upper bounds by [40].

State space formulation of MRM system

The rewriting dynamic model is obtained by substituting (19), (21), (22) and (23) into (17),

$$\begin{aligned} {\ddot{q}_i} = - {B_i}\left[ \begin{array}{l} {{{\hat{b}}}_{fi}}{{\dot{q}}_i} + ({{{\hat{f}}}_{ci}} + {{{\hat{f}}}_{si}}{e^{( - {{{\hat{f}}}_{\tau i}}\dot{q}_i^2)}})sgn({{\dot{q}}_i}) + {f_{qi}}({q_i},{{\dot{q}}_i})\\ + Y({{\dot{q}}_i}){{{\tilde{F}}}_i} + \frac{{{\tau _{fci}}}}{{{\gamma _i}}} + {U_{zi}} + {V_{zi}} - {\tau _i} - {\tau _{iext}} \end{array} \right] , \nonumber \\ \end{aligned}$$
(24)

where \({B_i} = {({I_{mi}}{\gamma _i})^{ - 1}}\). Meanwhile, redefine the state vector \({x_i} = {[{x_{i1}}{} {} {} {} {} {} {} {} {x_{i2}}]^{\mathrm{T}}} = {[{q_i}{} {} {} {} {} {} {} {} {\dot{q}_i}]^{\mathrm{T}}} \in {R^{2 \times 1}}\), and control input of the system is expressed as \({u_i} = {\tau _i} \in {R^{1 \times 1}}\). Then, the state space is formulated as:

$$\begin{aligned} {S_i}\left\{ \begin{array}{l} {{\dot{x}}_{i1}} = {x_{i2}}\\ {{\dot{x}}_{i2}} = {\phi _i}\left( {{x_i}} \right) + {\zeta _i}\left( x \right) + {B_i}{u_i} + \\ {y_i} = {x_{i1}} \end{array} \right. {B_i}{\tau _{iext}}, \end{aligned}$$
(25)

in (25), the precise modeling part is formulated as \({\phi _i}\left( {{x_i}} \right) = - {B_i}({{\hat{b}}_{fi}}{x_{i2}} + ({{\hat{f}}_{ci}} + {{\hat{f}}_{si}}{e^{( - {{{\hat{f}}}_{\tau i}}x_{i2}^2)}})sgn({x_{i2}}) + \frac{{{\tau _{fci}}}}{{{\gamma _i}}})\). \({\zeta _i}\left( x \right) = - {B_i}{f_{qi}}\left( {{x_{i1}},{x_{i2}}} \right) - {B_i}( {( {{x_{i2}}}){{{\tilde{F}}}_i}} ) - {B_i}( {{U_{zi}} + {V_{zi}}} )\) reflects the previous model uncertainty of the subsystem. The control objective of proposed control method is designed to keep interaction contact force within a safe range based on high performance trajectory tracking, which is described in the next section.

Human motion intention estimation-based decentralized robust interaction control

Harmonic drive compliance model-based human motion intention estimation

Different from the conventional position control that the desired trajectories are artificially given and predefined, the motion trajectory of MRM is unknown in pHRI, the ideal performance of pHRI should be determined by human motion intention. In this section, human motion intention is estimated by radial basis function neural network (RBFNN) via the human limb model. Meanwhile, the human limb can be expressed as an impedance model, which includes mass, damper and spring terms. As discussed and verified in [41], since the mass component is small, one can define the simplified human limb model as:

$$\begin{aligned} - {C_H}\dot{z} + {G_H}({z_{Hd}} - z) = {f_{ext}}, \end{aligned}$$
(26)

in (26) \({C_H}\),\({G_H}\) denote damper as well as spring matrices, \({z_{Hd}}\) denotes the estimated human motion intention at Cartesian space.

\({C_H}\), \({G_H}\) are unknown functions relate to z, \(\dot{z}\),i.e., rewriting \({C_H}\), \({G_H}\) as \({C_H}(z,\dot{z})\), \({G_H}(z)\) by referencing in [42]. Then, the motion intention \({z_{Hd}}\) can be rewritten as a nonlinear function related to the interaction contact force \({f_{ext}}\), actual velocity and position \(\dot{z}\), z i.e. as:

$$\begin{aligned} {z_{Hd}} = F({f_{ext}},z,\dot{z}), \end{aligned}$$
(27)

where \(F( \cdot )\) is unknown function, interaction contact force is formulated based on harmonic drive compliance model in (18). Meanwhile, one can present the following assumption for (27):

Assumption 1

The information determined the human motion intention in pHRI, which includes interaction contact force \(f_{ext}\), actual position z, and velocity \(\dot{z}\). Meanwhile, (27) is a continuous function.

\(F( \cdot )\) is a nonlinear continuous function caused by time-varying and uncertain properties of human limb. In addition, it is difficult to estimate \({z_{Hd}}\) in pHRI which caused by the human limb parameters are variable. Because of the approximation ability of RBFNN, the human motion intention can be estimated using RBFNN, and structure of the previous neural network is given by,

$$\begin{aligned} \begin{array}{l} \varphi (W,l) = {W^{\mathrm{T}}}H(l),W,S(l) \in {R^P},\\ H(l) = {[{h_1}(l),{h_2}(l),...,{h_p}(l)]^{\mathrm{T}}},\\ {h_k}(l) = \exp \left[ \frac{{ - {{(l - {\mu _k})}^{\mathrm{T}}}(l - {\mu _k})}}{{\eta _k^2}}\right] ,\\ k = 1,2,...,p, \end{array} \end{aligned}$$
(28)

in (28) \(\varphi (W,l)\) denotes the function relate to l, \(l \in {\Omega _l} \in {R^m}\) is neural network’s input, p reflects the nodes number of the RBFNN, \({\mu _k} = {[{\mu _{k,1}},{\mu _{k,2}},...,{\mu _{k,m}}]^{\mathrm{T}}}\) represents center field, \({\eta _k}\) reflects width of Gaussian function, and W denotes weight vector.

The actual value and estimation value of the motion intention are given by,

$$\begin{aligned} \begin{array}{l} {z_{Hd}} = W_{xi}^{\mathrm{T}}{H_i}({l_i}) + {\varepsilon _i},\\ {{{\hat{z}}}_{Hd}} = {\hat{W}}_{xi}^{\mathrm{T}}{H_i}({l_i}),i = 1,2,...,n, \end{array} \end{aligned}$$
(29)

in (29), \({l_i} = {[f_{exti}^{\mathrm{T}},z_i^{\mathrm{T}},\dot{z}_i^{\mathrm{T}}]^{\mathrm{T}}}\) reflects the same meaning as l in (28), \({\varepsilon _i}\) represents the estimation error, \({{\hat{W}}_{xi}}\) reflects the weight estimation, and \({H_i}\) denotes an activation function.

The control objective of pHRI is to ensure that robot track the human motion intention estimation actively and maintain the interact contact force \({f_{ext}}\) to minimum. By employing the gradient descent algorithm, one can formulate the updating law of the weight vector \({{\hat{W}}_{xi}}\) according to the following cost function:

$$\begin{aligned} {E_i} = \frac{1}{2}f_{exti}^2, \end{aligned}$$
(30)

where \({f_{exti}}\) reflects the ith joint interact force which is mapping by ith joint interact torque.

Then, one can define the weight update law \({{\dot{{\hat{W}}}_{xi}}}\) as follows:

$$\begin{aligned} \begin{array}{l} {{\dot{{\hat{W}}}_{xi}}}\left( t \right) = - {\beta _i}\frac{{\partial {E_i}}}{{\partial {{{\hat{W}}}_{xi}}}} = - {\beta _i}\frac{{\partial {E_i}}}{{\partial {f_{exti}}}}\frac{{\partial {f_{exti}}}}{{\partial {z_{Hd,i}}}}\frac{{\partial {z_{Hd,i}}}}{{\partial {{{\hat{W}}}_{xi}}}}\\ \mathrm{{{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{} {}{}{}{}{}{}{}{} = }} - {\beta _i}{f_{exti}}\frac{{\partial {f_{exti}}}}{{\partial {z_{Hd,i}}}}\frac{{\partial {z_{Hd,i}}}}{{\partial {{{\hat{W}}}_{xi}}}}, \end{array} \end{aligned}$$
(31)

where \({\beta _i}\) represents a positive scalar.

Equivalently, according to (26), \(\frac{{\partial {f_{exti}}}}{{\partial {z_{Hd,i}}}}\) is formulated as:

$$\begin{aligned} \frac{{\partial {f_{exti}}}}{{\partial {z_{Hd,i}}}} = {G_{H,i}}. \end{aligned}$$
(32)

At the same time, \(\frac{{\partial {z_{Hd,i}}}}{{\partial {{{\hat{W}}}_{xi}}}}\) can be obtained by (29) as:

$$\begin{aligned} \frac{{\partial {z_{Hd,i}}}}{{\partial {{{\hat{W}}}_{xi}}}} = {H_i}\left( {{l_i}} \right) . \end{aligned}$$
(33)

The weight update rate can be obtained by substituting (32),(33) into (31) as:

$$\begin{aligned} {\dot{{\hat{W}}}_{xi}}\left( t \right) = - {\alpha _i}{f_{exti}}{H_i}\left( {{l_i}} \right) , \end{aligned}$$
(34)

where \({\alpha _i} = {\beta _i}{G_{H,i}}\).

In addition, \({{\hat{W}}_{xi}}\left( t \right) \) is satisfied with (35),

$$\begin{aligned} {\dot{{\hat{W}}}_{xi}}\left( t \right) = {{\hat{W}}_{xi}}\left( 0 \right) - {\alpha _i}\int _0^t {\left( {{f_{exti}}\left( \omega \right) {H_i}\left( {{l_i}\left( \omega \right) } \right) } \right) } d\omega . \end{aligned}$$
(35)

Therefore, the human motion intention estimation in Cartesian coordinates \({{\hat{z}}_{Hd}}\) is received according to (29) by the previous equations.

Decentralized robust controller implementation

To the best of our knowledge, the position relationship between Cartesian coordinates and joint coordinates is given by (14), one can obtain human motion intention at joint coordinates \({\hat{q}} = {[{{\hat{q}}_1}, \cdots {{\hat{q}}_i}, \cdots {{\hat{q}}_n}]^T} \in {R^n}\) by substituting \({{\hat{z}}_{Hd}}\) into (36),

$$\begin{aligned} {\hat{q}}\left( t \right) = {\zeta ^{ - 1}}\left( {{{{\hat{z}}}_{Hd}}} \right) . \end{aligned}$$
(36)

One can design the decentralized interaction robust control strategy based on the torque-sensorless human motion intention estimation via harmonic drive compliance model, which utilized the dynamic model formulation of MRM in (17) and the state space formulation in (25).

It is obviously guaranteed that robot not only tracks \({{\hat{q}}_i}\), but also ensures the comfortable and safety in pHRI. Then, we can define the position error \({e_{qi}}\) and contact torque tracking error \({e_{\tau i}}\) as follows:

$$\begin{aligned} {e_{qi}} = {q_i} - {{\hat{q}}_i},{e_{\tau i}} = {J^{\mathrm{T}}}\left( q \right) \left( {{f_{exti}} - {f_{di}}} \right) , \end{aligned}$$
(37)

where \({f _{di}}\) reflects the desired constrained contact force of the end effector. we can obtain the control variable related to the position error \({e_{qi}}\) and contact torque tracking error \({e_{\tau i}}\) in (37) as follows:

$$\begin{aligned} {s_i} = {\dot{e}_{qi}} + {k_{pi}}{e_{qi}} + {k_{\tau i}}\int _0^t {{e_{\tau i}}d\tau }, \end{aligned}$$
(38)

where \({k_{pi}}\) and \({k_{\tau i}}\) are positive constants.

Remark 2

\({f_{di}}\) is the desired constrained contact force of end effector which is obtained from interact force \({f _{exti}}\) by saturation function method, and by considering the safety of human partner in pHRI, the threshold range of the saturation function is artificiality given for 7.5 Nm, the formulation of the desired constrained contact torque is given as \({f _{di}} = \left\{ \begin{array}{l} 7.5,\left| {{f _{exti}}} \right| > 7.5\\ {f _{exti}},\left| {{f _{exti}}} \right| \le 7.5 \end{array} \right. .\)

First, one can define the decentralized control law \({u_{i1}}\) which is exploited for compensating the precise modeling terms as:

$$\begin{aligned} {u_{i1}} = {{\hat{b}}_i}{\dot{q}_i} + \left( {{\hat{f}}_{ci}} + {{\hat{f}}_{si}}{e^{\left( - {{{\hat{f}}}_{\tau i}}\dot{q}_i^2\right) }}\right) sgn({\dot{q}_i}) + \frac{{{\tau _{fci}}}}{{{\gamma _i}}}. \end{aligned}$$
(39)

According to [40], the friction modeling uncertainty \({f_{qi}}\left( {{q_i},{{\dot{q}}_i}} \right) \) is considered as an unknown constant, which will be compensated by the integral compensator. Meanwhile, the parameter uncertainties are non-constant, which will be divided into an uncertain constant vector \({{\tilde{F}}_{ic}}\) and an uncertain variable \({{\tilde{F}}_{iv}}\).

$$\begin{aligned}&{{\tilde{F}}_i} = {{\tilde{F}}_{ic}} + {{\tilde{F}}_{iv}}, \end{aligned}$$
(40)
$$\begin{aligned}&\left| {{F_{ivn}}} \right| < \rho _n^i,n = 1,2,3,4. \end{aligned}$$
(41)

in which \(\rho _n^i = \left[ {{\rho _{i1}},{\rho _{i2}},{\rho _{i3}},{\rho _{i4}}} \right] \) is a determined constant vector.

Based on the idea of the decomposed control, the constant part \({{\tilde{F}}_{ic}}\) will be compensated by adaptive control and the variable part \({{\tilde{F}}_{iv}}\) will be dealed with by a robust item, the controller \(u_{if}\) is designed to deal with friction uncertainty as:

$$\begin{aligned} {u_{if}} = u_u^i + Y\left( {{{\dot{q}}_i}} \right) \left( {u_{pc}^i + u_{pv}^i} \right) , \end{aligned}$$
(42)

where \({f_{qi}}\left( {{q_i},{{\dot{q}}_i}} \right) \) is dealed with by \(u_u^i\), \(u_{pc}^i\) and \(u_{pv}^i\) are utilized to deal with \({{\tilde{F}}_{ic}}\) and \({{\tilde{F}}_{iv}}\), respectively. Since the compensation for friction is identical for each joint, and the specific formulation of \(u_{pc}^i,u_{pv}^i\) and \(u_u^i\) are given by:

$$\begin{aligned} \begin{array}{l} u_u^i = \left\{ \begin{array}{l} - {\rho _{fpi}}\frac{{{s_i}}}{{\left| {{s_i}} \right| }}{} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} \left| {{s_i}} \right|> {\varepsilon ^i}\\ - {\rho _{fpi}}\frac{{{s_i}}}{{{\varepsilon ^i}}}{} {} {} {} {} {} {} {} {} {} {} {} {}{}{} {}{}{} {}{} {} {} \left| {{s_i}} \right| \le {\varepsilon ^i} \end{array} \right. ,\\ u_{pc}^i = - k\displaystyle \int \limits _0^t {Y{{\left( {{{\dot{q}}_i}} \right) }^{\mathrm{T}}}} {s_i}d\tau ,\\ u_{pvn}^i = \left\{ \begin{array}{l} - \rho _n^i\frac{{\zeta _n^i}}{{\left| {\zeta _n^i} \right| }}{} {} {} {} {} {} {} {} {} {} {} {} {} {} {} \left| {\zeta _n^i} \right| > \varepsilon _{pn}^i\\ - \rho _n^i\frac{{\zeta _n^i}}{{\varepsilon _{pn}^i}}{} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} \left| {\zeta _n^i} \right| \le \varepsilon _{pn}^i \end{array} \right. ,\;n = 1,2,3,4, \end{array} \end{aligned}$$
(43)

where \(\zeta _n^i = Y{\left( {{{\dot{q}}_i}} \right) ^T}{s_i}\), \({\varepsilon ^i}\), \(\varepsilon _{pn}^i\) are positive control parameters. We can decompose the \({U_{zi}}\) and \({V_{zi}}\) as (44), and (45) illustrates the up-bounds of the variable terms \({U_{ziv}}\) and \({V_{ziv}}\).

$$\begin{aligned}&{U_{zi}} = {U_{zic}} + {U_{ziv}},\nonumber \\&{V_{zi}} = {V_{zic}} + {V_{ziv}}, \end{aligned}$$
(44)
$$\begin{aligned}&\left| {{U_{ziv}}} \right| \le {\rho _{Uiv}},\left| {{V_{ziv}}} \right| \le {\rho _{Viv}}. \end{aligned}$$
(45)
Fig. 2
figure 2

The whole control diagram

Adopt the same ideas as in \({{\tilde{F}}_i}\), the uncertainty constant parts \({U_{zic}}\) and \({V_{zic}}\) will be compensated by the adaptive compensators \(u_{uicn}^i\) and \(u_{vicn}^i\), respectively, the robust controllers \(u_{uivn}^i\), \(u_{vivn}^i\)are designed to compensate the uncertainty variable part \({U_{ziv}}\), \({V_{ziv}}\), respectively.

$$\begin{aligned} \begin{array}{l} {u_{uicn}} = - {k_{icn}}\displaystyle \int \limits _0^t {{k_{icn}}{s_i}d\tau },\\ {u_{vicn}} = - {k_{ivn}}\displaystyle \int \limits _0^t {{k_{ivn}}{s_i}d\tau }, \end{array} \end{aligned}$$
(46)
$$\begin{aligned} \begin{array}{l} {u_{uivn}} = \left\{ \begin{array}{l} - {\rho _{Uiv}}\frac{{{\upsilon _{ivn}}}}{{\left| {{\upsilon _{ivn}}} \right| }}\mathrm{{ } {} {} {} }\left| {{\upsilon _{ivn}}} \right|> \varepsilon _{Un}^i\\ - {\rho _{Uiv}}\frac{{{\upsilon _{ivn}}}}{{\varepsilon _{Un}^i}}\mathrm{{ } {} {} {} {} {} {} {} {} }\left| {{\upsilon _{ivn}}} \right| \le \varepsilon _{Un}^i \end{array} \right. ,\\ {u_{vivn}} = \left\{ \begin{array}{l} - {\rho _{Viv}}\frac{{{\chi _{viv}}}}{{\left| {{\chi _{viv}}} \right| }}\mathrm{{ } {} {} {} }\left| {{\chi _{viv}}} \right| > \varepsilon _{vn}^i\\ - {\rho _{Viv}}\frac{{{\chi _{viv}}}}{{\varepsilon _{vn}^i}}\mathrm{{ } {} {} {} {} {} {} {} {} }\left| {{\chi _{viv}}} \right| \le \varepsilon _{vn}^i\mathrm{{ }} \end{array} \right. , \end{array} \end{aligned}$$
(47)

where \({\upsilon _{ivn}} = {k_{icn}}{s_i}\), \({\chi _{viv}} = {k_{ivn}}{s_i}\) and \(\varepsilon _{pin}^i\), \(\varepsilon _{pin}^i\) are positive control parameters.

We can obtain the control law \({u_{i2}}\) by combining with (42), (43), (44), (46) and (47), which compensates the precise model uncertainties.

$$\begin{aligned} {u_{i2}}&= u_u^i + Y\left( {{{\dot{q}}_i}} \right) \left( {u_{pc}^i + u_{pv}^i} \right) + {u_{uicn}}\nonumber \\&\quad + {u_{uivn}} + {u_{vicn}} + {u_{vivn}}. \end{aligned}$$
(48)

Decentralized interaction controller implementation

Since desired position trajectory of MRM is predefined in traditional position control task, the desired position trajectory is unknown in pHRI, which is dominated by human in this paper. The planned trajectory of robot is the human motion intention generation when human interact with robot. However, it is circumscribed that traditional position control method can only provide the trajectory tracking, human motion interaction control method can guarantee the accuracy of tracking both position and the convergence of interaction force to a safe range.

According to (38), the time derivative of \({s_i}\) is formulated as follows:

$$\begin{aligned} {\dot{s}_i} = {\ddot{e}_{qi}} + {k_{pi}}{\dot{e}_{qi}} + {k_{\tau i}}{e_{\tau i}}. \end{aligned}$$
(49)

Substituting (49) into (25), we can obtain as:

$$\begin{aligned} {\dot{s}_i} = {\phi _i}\left( {{x_i}} \right) + {\zeta _i}\left( x \right) + {B_i}{u_i} + {B_i}{\tau _{iext}} + {v_i}, \end{aligned}$$
(50)

where \({v_i} = - {\ddot{{\hat{q}}}_i} + {k_{pi}}{\dot{e}_{qi}} + {k_{\tau i}}{e_{\tau i}}\).

Since the control objective of the decentralized human–robot interaction control method is not only satisfied the accuracy of position tracking, but also guarantee the interaction force convergence to an acceptable range. Meanwhile, one can design the controller as:

$$\begin{aligned} {u_{i3}} = {I_{mi}}{\gamma _i}\left( {{{\ddot{{\hat{q}}}}_i} - {k_{pi}}{{\dot{e}}_{qi}} - {k_{\tau i}}{e_{\tau i}}} \right) - {\tau _{iext}}. \end{aligned}$$
(51)

As the human motion intention estimation is obtained in the previous section, the decentralized interaction control law is designed based on the safety contact force to guarantee that the convergence of interaction force within the safety threshold, so that the robot can show the performance of compliance. The flow diagram of proposed control diagram in this paper is show in Fig. 2.

By combining (39) (48) and (51), one can design the decentralized interaction robust control law as follows:

$$\begin{aligned} \begin{array}{l} {u_i} = {u_{i1}} + {u_{i2}} + {u_{i3}} = {{{\hat{b}}}_i}{{\dot{q}}_i} + ({{{\hat{f}}}_{ci}} + {{{\hat{f}}}_{si}}{e^{( - {{{\hat{f}}}_{\tau i}}\dot{q}_i^2)}})sgn({{\dot{q}}_i})\\ {}\mathrm{{ }} + u_u^i + Y\left( {{{\dot{q}}_i}} \right) \left( {u_{pc}^i + u_{pv}^i} \right) + {u_{uicn}} + {u_{uivn}} + {u_{vicn}}\\ {}\mathrm{{ }} + {u_{vivn}} + {I_{mi}}{\gamma _i}\left( {{{\ddot{{\hat{q}}}}_i} - {k_{pi}}{{\dot{e}}_{qi}} - {k_{\tau i}}{e_{\tau i}}} \right) - {\tau _{iext}} + \frac{{{\tau _{fci}}}}{{{\gamma _i}}}. \end{array}\nonumber \\ \end{aligned}$$
(52)

Stability proof

Theorem 1

Consider an MRM under pHRI, with the subsystem dynamic model formulated in (17), the model uncertainties existed in (19), (21) and (22), if the human motion intention is estimated by using (29),then the weight error of the human motion intention estimation, the contact torque tracking error and the position tracking error of each robotic joint are UUB under the proposed decentralized interaction robust control law (52).

Proof

Select Lyapunov function

$$\begin{aligned} {V_i} = \frac{1}{2}{M_i}s_i^2 + \frac{1}{2}k\psi _i^{\mathrm{T}}{\psi _i} + \eta _i^{\mathrm{T}}{\eta _i} + \Lambda _i^{\mathrm{T}}{\Lambda _i} + \frac{1}{2}{k_2}{\tilde{W}}_{xi}^{\mathrm{T}}{{\tilde{W}}_{xi}}, \nonumber \\ \end{aligned}$$
(53)

where \({\psi _i} = \frac{1}{{{k_1}}}{{\tilde{F}}_{ci}} + \int \nolimits _0^t {Y{{\left( {{{\dot{q}}_i}} \right) }^{\mathrm{T}}}{s_i}d\tau },{\eta _i} = \frac{1}{{{k_{icn}}}}{U_{zic}} + \int \nolimits _0^t {{k_{icn}}{s_i}d\tau },\) and \({\Lambda _i} = \frac{1}{{{k_{ivn}}}}{V_{_{zic}}} + \int \nolimits _0^t {{k_{ivn}}{s_i}d\tau }.\) \({{\dot{\psi }} _i} = Y{\left( {{{\dot{q}}_i}} \right) ^{\mathrm{T}}}{s_i},{{\dot{\eta }} _i} = {k_{icn}}{s_i},{{\dot{\Lambda }} _i} = {k_{ivn}}{s_i}\). \({M_i} = {I_m}_i{\gamma _i},k,{k_2},{{\tilde{F}}_{ci}},{U_{zic}},\) and \({V_{zic}}\) are constants. Differentiating (53) yields,

$$\begin{aligned} \begin{array}{l} {{\dot{V}}_i} = {k_1}\left( {\frac{1}{{{k_1}}}{{{\tilde{F}}}_{ci}} + \displaystyle \int \limits _0^t {Y{{\left( {{{\dot{q}}_i}} \right) }^{\mathrm{T}}}{s_i}d\tau } } \right) \left( {Y{{\left( {{{\dot{q}}_i}} \right) }^{\mathrm{T}}}{s_i}} \right) + {k_2}{\tilde{W}}_{xi}^{\mathrm{T}}{{\dot{{\tilde{W}}}}_{xi}}\\ {}\mathrm{{ }} \qquad + \left( {\frac{1}{{{k_{icn}}}}{U_{zic}} + \displaystyle \int \limits _0^t {{k_{icn}}{s_i}d\tau } } \right) \\ \left( {{k_{icn}}{s_i}} \right) + \left( {\frac{1}{{{k_{ivn}}}}{V_{zic}} + \displaystyle \int \limits _0^t {{k_{ivn}}{s_i}d\tau } } \right) \left( {{k_{ivn}}{s_i}} \right) \\ {}\mathrm{{ }} + {s_i}\!\!\left( \begin{array}{l} \left( {{u_{ui}} - {f_{qi}}\left( {{q_i},{{\dot{q}}_i}} \right) } \right) + Y\left( {{{\dot{q}}_i}} \right) \left( {u_{pc}^i - {{{\tilde{F}}}_{ci}}} \right) \\ + Y\left( {{{\dot{q}}_i}} \right) \left( {u_{pv}^i - {{{\tilde{F}}}_{vi}}} \right) \\ \quad + \left( {{u_{uicn}} + {u_{uivn}}} \right) + \left( {{u_{vicn}} + {u_{vivn}}} \right) \\ - \left( {{U_{zic}} + {U_{ziv}}} \right) - \left( {{V_{zic}} + {V_{ziv}}} \right) \end{array} \right) \\ {}\mathrm{{ }} = {s_i}\left( {{u_{ui}} - {f_{qi}}\left( {{q_i},{{\dot{q}}_i}} \right) } \right) + {s_i}Y\left( {{{\dot{q}}_i}} \right) \left( {{u_{pvi}} - {{{\tilde{F}}}_{vi}}} \right) \\ \qquad + {s_i}\left( {{u_{uivn}} - {U_{ziv}}} \right) \\ {}\mathrm{{ }} \quad + {s_i}\left( {{u_{vivn}} - {V_{ziv}}} \right) + {k_2}{\tilde{W}}_{xi}^{\mathrm{T}}{{\dot{{\tilde{W}}}}_{xi}}. \end{array}\nonumber \\ \end{aligned}$$
(54)

According to (43), (46) and (47), we obtain the following relations:

$$\begin{aligned}&\left\{ \begin{array}{l} {s_i}Y{\left( {{{\dot{q}}_i}} \right) ^{\mathrm{T}}}\left( { - {{{\tilde{F}}}_{vi}} + {u_{pvi}}} \right) < 0, if\left| {{s_i}Y{{\left( {{{\dot{q}}_i}} \right) }^{\mathrm{T}}}} \right| > {\varepsilon _{pni}},\\ {s_i}Y{\left( {{{\dot{q}}_i}} \right) ^{\mathrm{T}}}\left( { - {{{\tilde{F}}}_{vi}} + {u_{pvi}}} \right) \le {s_i}Y{({{\dot{q}}_i})^{\mathrm{T}}}\sum \limits _{n = 1}^4 {\left( {\rho _n^i\frac{{{s_i}Y{{({{\dot{q}}_i})}^{\mathrm{T}}}}}{{\left| {{s_i}Y{{({{\dot{q}}_i})}^{\mathrm{T}}}} \right| }} - \rho _n^i\frac{{{s_i}Y{{({{\dot{q}}_i})}^{\mathrm{T}}}}}{{{\varepsilon _{pni}}}}} \right) }, \\ if\left| {{s_i}Y{{({{\dot{q}}_i})}^{\mathrm{T}}}} \right| \le {\varepsilon _{pni}}, \end{array} \right. \nonumber \\ \end{aligned}$$
(55)
$$\begin{aligned}&\left\{ \begin{array}{l} {s_i}\left( {{u_{uivn}} - {U_{ziv}}} \right) = {s_i}\left( { - {\rho _{Uiv}}\frac{{{\upsilon _{ivn}}}}{{\left| {{\upsilon _{ivn}}} \right| }}{} - {U_{ziv}}} \right)< 0, if\left| {{k_{icn}}{s_i}} \right| > \varepsilon _{Un}^i,\\ {s_i}\left( {{u_{uivn}} - {U_{ziv}}} \right) < {s_i}\left( {{\rho _{Uiv}}\frac{{{\upsilon _{ivn}}}}{{\left| {{\upsilon _{ivn}}} \right| }} - {\rho _{Uiv}}\frac{{{\upsilon _{ivn}}}}{{\varepsilon _{Un}^i}}} \right) , if\left| {{k_{icn}}{s_i}} \right| \le \varepsilon _{Un}^i, \end{array} \right. \nonumber \\ \end{aligned}$$
(56)
$$\begin{aligned}&\left\{ \begin{array}{l} {s_i}\left( {{u_{vivn}} {-} {V_{ziv}}} \right) {=} {s_i}\left( { - {\rho _{Viv}}\frac{{{\chi _{viv}}}}{{\left| {{\chi _{viv}}} \right| }} {-} V_{kjc}^i} \right) < 0, if\left| {{k_{ivn}}{s_i}} \right| > \varepsilon _{vn}^i,\\ {s_i}\left( {{u_{vivn}} - {V_{ziv}}} \right) {=} {s_i}\left( {{\rho _{Viv}}\frac{{{\chi _{viv}}}}{{\left| {{\chi _{viv}}} \right| }} {-} {\rho _{Viv}}\frac{{{\chi _{viv}}}}{{\varepsilon _{pin}^i}}} \right) ,if\left| {{k_{ivn}}{s_i}} \right| \le \varepsilon _{vn}^i. \end{array} \right. \nonumber \\ \end{aligned}$$
(57)

\(\square \)

Fig. 3
figure 3

Experiment platform

Fig. 4
figure 4

pHRI experiment based on human motion intention estimation via harmonic drive compliance model

For the NN items, one obtain:

$$\begin{aligned} \begin{array}{l} {k_2}{\tilde{W}}_{xi}^{\mathrm{T}}{{\dot{{\tilde{W}}}}_{xi}} = {k_2}\left( {{e_{xiH}}{\tilde{W}}_{xi}^{\mathrm{T}}{S_i}\left( {{l_i}} \right) - {{\left\| {{\tilde{W}}_{xi}^{\mathrm{T}}{S_i}\left( {{l_i}} \right) } \right\| }^2}} \right) \\ {}\mathrm{{ }} \le {k_2}\left( {\frac{1}{2}e_{xiH}^2 - \frac{1}{2}{{\left\| {{\tilde{W}}_{xi}^{\mathrm{T}}{S_i}\left( {{l_i}} \right) } \right\| }^2}} \right) , \end{array} \end{aligned}$$
(58)

if \({{\tilde{W}}_{xi}}\) follows the constraints which is given by (59), (58) is negative:

$$\begin{aligned} {\Omega _x} = \left\{ {{{{\tilde{W}}}_{xi}}:\left\| {{{{\tilde{W}}}_{xi}}} \right\| \le \left\| {\frac{{{e_{xM}}}}{{{S_{xM}}}}} \right\| } \right\} . \end{aligned}$$
(59)

At \(\left| {\zeta _n^i} \right| \le \frac{{\varepsilon _{pn}^i}}{2}\), \(\left| {\upsilon _{jn}^i} \right| \le \frac{{\varepsilon _{Un}^i}}{2}\) and \(\left| {{\chi _{viv}}} \right| \le \frac{{\varepsilon _{vn}^i}}{2}\), the last terms of (55), (56), and (57) become the maximum value, and we have

$$\begin{aligned} \dot{V}&\le \sum \limits _{n = 1}^4 {\left( {\frac{{\rho _n^i{\varepsilon _{pni}}}}{4}} \right) } + \frac{{{\rho _{fpi}}{\varepsilon _i}}}{4} + \frac{{{\rho _{Uiv}}\varepsilon _{Un}^i}}{4} + \frac{{{\rho _{Viv}}\varepsilon _{vn}^i}}{4}\nonumber \\&\quad + \frac{1}{2}e_{xiH}^2 - \frac{1}{2}{\left\| {{\tilde{W}}_{xi}^{\mathrm{T}}{S_i}\left( {{l_i}} \right) } \right\| ^2}. \end{aligned}$$
(60)

According to (60), Lyapunov function can only be found if satisfied the following equation.

$$\begin{aligned} \left| {{s_i}} \right| \ge \sqrt{\left( {\sum \limits _{n = 1}^4 {\left( {\rho _n^i{\varepsilon _{pni}}} \right) + {\rho _{Viv}}\varepsilon _{vn}^i} + {\rho _{Uiv}}\varepsilon _{Un}^i + {\rho _{fpi}}{\varepsilon _i}} \right) } . \nonumber \\ \end{aligned}$$
(61)

Definition:

$$\begin{aligned} s_i^2 \le \left( {\sum \limits _{n = 1}^4 {\left( {\rho _n^i{\varepsilon _{pni}}} \right) } + {\rho _{Uiv}}\varepsilon _{Un}^i + {\rho _{Viv}}\varepsilon _{vn}^i + {\rho _{fpi}}{\varepsilon _i}} \right) . \nonumber \\ \end{aligned}$$
(62)

Then, on the surface \({s_{ir}}\) and \(\partial {s_{ir}}\) above, we can get:

$$\begin{aligned} \dot{V} \le - \left( {\left[ {\sum \limits _{n = 1}^4 {\left( {\rho _n^i{\varepsilon _{pni}}} \right) + {\rho _{Uiv}}\varepsilon _{Un}^i} + {\rho _{Viv}}\varepsilon _{vn}^i + {\rho _{fpi}}{\varepsilon _i}} \right] /4} \right) . \nonumber \\ \end{aligned}$$
(63)

When the root locus intersects the surface \(\partial {s_{ir}}\) at \({T_{ir}}\), we have:

$$\begin{aligned}&{V_i}\left( {{r_{ir}}\left( {{T_{ir}}} \right) } \right) - {V_i}\left( {{r_{ir}}\left( 0 \right) } \right) \nonumber \\&\quad \le - \left[ {\left( \begin{array}{l} \sum \limits _{n = 1}^4 {(\rho _n^i{\varepsilon _{pni}})} + {\rho _{Uiv}}\varepsilon _{Un}^i\\ + {\rho _{Viv}}\varepsilon _{vn}^i + {\rho _{fpi}}{\varepsilon _i} \end{array} \right) /4} \right] {T_{ir}}. \end{aligned}$$
(64)

According to (38), the previous \({s_i}\) upper bound that can be represented by \({e_{qi}},{\dot{e}_{qi}}\) and \({e_{\tau i}}\) This completes the proof.

Table 1 Parameters setting
Table 2 Performance comparisons in each control method

Experiments

In this section, the proposed control scheme is verified by the established experimental platform, and the experiment platform is illustrated in Fig. 3. We developed a modular robot with two degrees of freedom under pHRI, the joint module of the robot integrates DC motor, incremental encoder to measure motor displacement, absolute encoder to measure connecting rod displacement and HD device. The experimental platform is equipped with computer host, monitor, QPIDE data acquisition card and power amplifier device, and experiments are carried out through Simulink software. In these experiments, the harmonic drive compliance model is utilized to estimate the human motion intention when human interact with MRM, as shown in Fig. 4. Table 1 illustrates the parameters of the MRM system, which includes model parameters, controller parameters and upper bound parameters of the model uncertainty.

In this part, the experiment results illustrate to discuss the performance of harmonic drive compliance model-based torque estimation and human motion intention estimation, position tracking, tracking error, control torque and human–robot interaction contact force, respectively. Three types of control methods are utilized to confirm the effectiveness of the proposed control method which includes the human motion intention estimation-based existed decentralized robust control method via six-dimensional force sensor, e.g., [43] and the existed decentralized integral sliding mode control based on human motion intention estimation via six-dimensional force sensor, e.g., [44] as well as the proposed human motion intention estimation-based decentralized robust interaction control method via harmonic drive compliance model. All these control methods are used for the same type of pHRI task, i.e., handshaking. Table 2 lists the performance comparisons for each control method.

Fig. 5
figure 5

Joint torque and torque estimation error curves under the proposed method

Fig. 6
figure 6

Trajectory curves of human–robot interaction point in three-dimensional space under the proposed method

(1) Harmonic drive compliance model-based torque estimation and human motion intention estimation

Figure 5 shows joint torque curves under handshaking task that using the proposed harmonic drive compliance model-based joint torque estimation method. Meanwhile, joint torque estimation error is represented in the same figure, which can analyze the performance of joint torque estimation clearly. As we can see from the experimental results, it is concluded that estimated joint torque is highly consistent with the torque sensing measurements. One observes that the joint torque curves exist jump phenomenon at the inflection point of the motor inversion. Fortunately, the proposed harmonic drive compliance-based joint torque estimation method, which can effectively capture the instantaneous variation trend of the joint torque. Figure 6 reflects the human motion intention estimation curve under hand shaking pHRI task via the proposed harmonic drive compliance model. Since the experimental platform shown in Fig. 3, which is a 2-DOF MRM, and the joint axis is parallel assembly, the curve of human motion intention estimation is represented as a plane track in three-dimensional space. We can see in Fig. 6 that human motion intention can be well estimated at the inflection point, which indicates that the robot can ‘actively’ complete the task of following the human hand. Therefore, the trajectory of the robotic end-effector can track human motion intention estimation precisely via harmonic drive compliance model.

(2) Position tracking performance

Figures 7, 8, 9, 10, 11, 12 illustrate the position tracking and tracking error curves in joint space under handshaking via the two existing control methods as well as the proposed control method, respectively. From the detailed drawing of the trajectory tracking figures and the position tracking error figures, we can get the conclusion that the tracking performance of the proposed control method is better than other existing control method. Since the performance of the six-dimensional force sensor, tracking errors of these two existing methods possess jump phenomenon when the direction of the human motion intention changes. The proposed human motion intention estimation method is based on the harmonic drive compliance model, which possesses more accurate identification performance. Therefore, the availability and superiority of the proposed harmonic drive compliance model-based human motion intention estimation approach are verified.

Fig. 7
figure 7

Joint position tracking curves under existing decentralized robust control method

Fig. 8
figure 8

Joint position error curves under existing decentralized robust control method

Fig. 9
figure 9

Joint position tracking curves under existing decentralized integral sliding mode control method

Fig. 10
figure 10

Joint position error curves under existing decentralized integral sliding mode control method

Fig. 11
figure 11

Joint position tracking curves under proposed decentralized robust interaction control method

Fig. 12
figure 12

Joint position error under proposed decentralized robust interaction control method

Fig. 13
figure 13

Control torque curves under existing decentralized robust control method

Fig. 14
figure 14

Control torque curves under existing integral sliding mode control method

Fig. 15
figure 15

Control torque curves under proposed decentralized robust interaction control method

Fig. 16
figure 16

Human–robot interaction contact force curves under existing decentralized robust control method

Fig. 17
figure 17

Human–robot interaction contact force curves under existing integral sliding mode control method

Fig. 18
figure 18

Human–robot interaction contact force curves under proposed decentralized robust interaction control method

(3) Control torque

Figures 13, 14 and 15 represent the control torque curves under handshaking via the existing decentralized robust control method, decentralized integral sliding mode control method and the proposed control method, respectively. As we can see from these figures, when the trajectory changes fiercely, the control torque will increase instantly, which will affect the durability of the DC motor. As shown in Fig. 13, since the strong robustness of the controller, the control torque curve possesses serious chattering effect. The control torque curve is smooth because of the property of the integral sliding mode that can attenuate the chattering effect of the controller, which is illustrated in Fig. 14. As we can see in Fig. 15, using the proposed control method, the instantaneous growth of control torque is kept within the safe range. Meanwhile, the chattering effect of the proposed controller is decreased. Therefore, the proposed controller is superior.

(4) Human–robot interaction contact force

Figures 16, 17 and 18 show the human–robot interaction contact force curves under handshaking via the existing decentralized robust control method, the existing decentralized integral sliding model control method as well as the proposed control method. Since the experimental platform is assembled in parallel by a 2-DOF robot, since the contact force of the end-effector along the Z-axis has no component, the trajectory of the end effector is two dimensional. The assignment of human–robot interaction contact force of the proposed decentralized robust interaction control method not only converges within the acceptable threshold range, but also has no obvious chattering effect. Since these existing two control methods only control the trajectory tracking, the human–robot interaction contact force curve of the existing decentralized robust control method is far beyond the safe range. In addition, the proposed decentralized robust interaction control method is based on the harmonic drive compliance model, so that the human–robot interaction contact force is relatively smooth and without strong chattering phenomenon.

Conclusion

A novel torque-sensorless human motion intention estimation-based decentralized robust interaction control method is proposed for controlling HD-based MRM using only position measurements on the motor and link sides of each module in pHRI, which is utilized to solve the problem of unknown human motion intention. The human motion intention estimation of Cartesian space is reflected by interaction force which applied harmonic drive compliance model, and the human motion intention estimation of each subsystem is obtained by the mapping matrix from Cartesian space to joint space. The estimated motion intention of joint space has been utilized to design the decentralized robust interaction controller of MRM, as a result, MRM can actively follow the motion of human, and the ’safety’ interaction environment is developed. The uniformly ultimately boundedness (UUB) of the tracking error is proved by the Lyapunov theory. The advantage of the proposed method is verified by pHRI experiment.