Bayesian Estimation of Human Impedance and Motion Intention for Human-Robot Collaboration

—This paper proposes a Bayesian method to acquire the estimation of human impedance and motion intention in a human-robot collaborative task. Combining with prior knowledge of human stiffness, estimated stiffness obeying Gaussian distribution is obtained by Bayesian estimation and human motion intention can be also estimated. An adaptive impedance control strategy is employed to track a target impedance model and neural networks are used to compensate for uncertainties in robotic dynamics. Comparative simulation results are carried out to verify the effectiveness of estimation method and emphasize the advantages of the proposed control strategy. The experiment, performed on Baxter r robot platform, illustrate a good system performance.

Bayesian Estimation of Human Impedance and Motion Intention for Human-Robot Collaboration that cooperator will need to know motion intention of initiator and adapt to movement and interaction force of initiator. Obviously, collaborative robots, which are centered on human task requirement, have the ability to assist human partner and supersede cooperator's work in such kind of tasks [5]- [8].
Let us consider a classical physical human-robot interaction (pHRI) scenario as in Fig. 1. Abundant control strategies are developed for pHRI [9], [10] and various adaptive or learning control strategies also draw much attention from scholars [11]- [16]. Impedance control, first proposed by Hogan [17], is used to relate interactive force with deviations from desired states. Adaptive impedance control methods are proposed subsequently, e.g., [18]- [21]. Compared with hybrid force/position control, impedance control shows better robustness and does not need transitions between contact and noncontact situations. Although traditional impedance control has shown good performance in pHRI [22], it only enables human to change the robot's actual trajectory but not the robot's desired trajectory [23]. If the robot has knowledge of human motion intention [24], it can regard human motion intention as its own desired trajectory and human will cost less effort to accomplish the task. In [25], human motion intention has been estimated by online neural networks (NNs) based on available sensory information, an updating law is designed and the robot moves to time-varying human's intended position actively. In [26], an inversion-based approach is proposed to estimate the human intent by demonstration and it is used in input updating for improving trajectory tracking accuracy. The effectiveness of human-guided iterative learning control has been proven by human-in-loop trajectory tracking experiment. Mainprice et al. [27] proposed a method to predict the next movement of the human partner who is collaborating with robot by applying inverse optimal control and goal set iterative replanning. In [28], human motion intention is identified to enable the robot to follow human compliantly in fast point-to-point tasks.
When robot interacts with human in a constrained motion form, an estimation method of human impedance should be considered for improving the system stability during pHRI [29]- [32]. By tuning a target impedance based on human impedance estimation, variable target impedance parameters extend the robot learning skills beyond trajectory tracking, in which robot is gifted with submissive performance and more advanced skills that involve, among others, contacting with human partner. Some common contact impedance estimation methods are analyzed in [33], which include recursive least-squares method, model reference and indirect adaptive method, and signal processing method. Using information extracted from programming by human demonstration, Rozo et al. [34] proposed a method to estimate environmental stiffness which is obtained according to the covariance of the Gaussian mixture model. In [35], the tutor transfers a specific sawing skill to the robot successfully, by using electromyography (EMG) signals to estimate tutor stiffness in pHRI. In [36], desired impedance parameters are obtained based on the gradient-following and betterment methods. In [37], the optimal desired stiffness is designed by using human operator's EMG signals in an upper limb robotic exoskeleton application. In [38] and [39], in order to estimate human impedance characteristics, a small external perturbation to the human arm is required in the cooperative task.
The abundant control strategies of nonlinear systems are proposed in recent years [40]- [45]. The model-based control strategies have more precise tracking capacities than classical PID control. In addition, it can avoid spending time to find the proper PID values of gains. The researches on adaptive control also draw much attention [46]- [48]. However, uncertainties in model dynamics are ubiquitous [49], [50] and have attracted attention of researchers [51]- [55]. In [56], radial basis function NNs (RBFNNs) are used to handle uncertainties in robotic dynamics, and the backstepping method is used to design a stable controller. This RBFNN method has been used in the applications of robotic flexible joints [57], output and input constraints [58]- [63], and teleoperation [64]. In [65], NNs are employed to compensate for uncertainties in the presence of unknown dynamics of both the grasped object and dual robotic manipulators. A switching method is integrated into controller to achieve global stability. In [66], an adaptive robust control design is proposed for multiple mobile manipulators, a common object in contact with a rigid surface is grasped by multiple mobile manipulators and they show robustness not only to parametric uncertainties but also to external disturbances. Some observer-based adaptive control strategies are also proposed for solving unknown disturbance or unknown states [67]- [72].
The Bayesian estimated methods are widely utilized in dealing with uncertainties in robot motion planning [73] and robot visual tracking [74]. Some works have been done about tactile perception in recent years [75]. In this article, a Bayesian method is proposed for human impedance and motion intention estimation, and the neural impedance control strategy is used to achieve efficient human-robot cooperation.
The construction of this article is described as follows. In Section II, the dynamics of human and robot are presented and the task objective is introduced. In Section III, a Bayesian estimation method is employed in human stiffness estimation, and the human motion intention is estimated according to the dynamic relationship between human stiffness and motion intention. In Section IV, impedance control is analyzed, NNs are used to handle model uncertainties in control design, and stability analysis is proven by constructing Lyapunov function candidates. In Section V, comparative simulations are carried out to show the advancement of our proposed method. In Section VI, an experiment is designed to evaluate the Fig. 1. Scenario where human and robot collaborate to perform an object transporting task. Human is an initiator of the task, that is, human will lead the task and he/she knows the task target position, and robot will be obedient completely to help human to finish the task, that is, robot will be a cooperator.
performance of our controller design on Baxter robot platform. In Section VII, the conclusion is presented.

II. PROBLEM FORMULATION
In this article, we consider an object transporting task as shown in Fig. 1. In this task, human will lead by applying an interaction force to the object and robot will cooperate with human to lift the object on the other side.

1) Robot's Dynamic Model:
We consider the robot as an m-DOF rigid manipulator, so the robotic dynamics in the joint space can be described as follows: where q,q,q ∈ R m are the joint angle, velocity, and acceleration vectors, respectively. M(q) ∈ R m×m is the symmetric and positive-definite inertia matrix, C(q,q)q ∈ R m is Coriolis and centripetal vector, G(q) ∈ R m denotes the gravity vector, τ ∈ R m denotes the control input vector, f r ∈ R h is the vector of the interaction force between the robot and the transferred object, and J(q) ∈ R h×m is the Jacobian matrix, where h denotes the dimension in the Cartesian space. The forward kinematics of the robot is given by x = (q), differentiating x with respect to time we obtainẋ = J(q)q. Based on the inverse kinematics,q andq in the joint space can be described aṡ where J −1 (q) denotes the inverse of J(q), and x,ẋ,ẍ ∈ R h denote the position, velocity, and acceleration vectors in the Cartesian space, respectively. By substituting (2) into (1), we obtain the robot's dynamic model in the Cartesian space as follows: where the inertia matrix M r (x) ∈ R h×h , the Coriolis and centripetal force vector C r (x,ẋ)ẋ ∈ R h , the gravitational force vector G r (x) ∈ R h , and the control force vector u ∈ R h in the Cartesian space in (3) can be calculated as 2) Human's Dynamic Model: In pHRI, the dynamic model of human in the Cartesian space in h dimension can be simply described as a spring model where K h ∈ R h×h denotes the human's stiffness matrix, x h ∈ R h denotes the human's target position vector in h dimension, that is, human motion intention, x ∈ R h denotes the actual position, and f h ∈ R h denotes the interaction force vector between human and transferred object.

B. Task Objective
In this task, the most important problems are how to acquire human stiffness and how to obtain human motion intention in (5). If robot knows human motion intention and human stiffness, it will be convenient to design impedance controller for efficient human-robot interaction. In our task, we want to make human and robot act with a same behavior for performing tasks successfully. If they have different behaviors during a cooperative task, the task will be inefficient or unsuccessful. The same behavior means that the robot and human have a same initial position and a same moving target position, and the human's stiffness matrix K h should be same as the robot's stiffness matrix K d . Therefore, we first design a target impedance model for the robot, which is described as follows: where d is the desired inertia matrix, D d is the desired damper matrix, K d is the desired stiffness matrix, and x d denotes the robot's desired target position. Considering a slow speed human-robot interactive process, (6) can be simplified as becauseẋ andẍ are close to zero. The simplified target impedance model (7) shows the dynamic relationship between displacement and interaction force clearly. As it can be seen from (5), in this cooperative object transporting task, we should design the robot desired target position x d as human motion intention x h and design K d as human stiffness K h . However, human motion intention x h and human stiffness K h are unknown to robot. Therefore, we need to propose an estimation method to obtain an estimate of human motion intentionx h and an estimate of human stiffnessK h . We can write the estimate of (5) in one dimension as follows: whereK h1 ,x h1 , andf h1 denote the estimates of human stiffness parameter, human target position, and interaction force between human and object in one dimension, respectively. In this article, we regard the transporting object as a mass point of which the tiny mass and volume can be ignored. Therefore, the interaction force between human and transferred object f h is the same as the interaction force between robot and transferred object f r . This leads to a scenario where human and robot has a direct physical contact. When we measure f r by the force sensor mounted on the end-effector of robot, f h can be obtained.
The control architecture is shown in Fig. 2. In the following two sections, we first explain how to estimate human's target position and stiffness, and then design a controller to achieve desired robot's impedance.

III. HUMAN STIFFNESS LEARNING AND MOTION INTENTION ESTIMATION
The Bayesian parameter estimation method is an important method to estimate unknown parameters. We use this method to get the estimation of K h1 and x h1 .
First, we establish a quadratic cost function to evaluate the estimation accuracy as follows: ] can be regarded as K h according to (5), so we can use (9) to evaluate the estimation accuracy of K h . f h1 andẋ 1 can be measured by force and velocity sensors, respectively.
Remark 2: Similar idea has been used in [76] for estimating human stiffness in real time.
We assume that ] obeys the following distribution: where N( * ) denotes the Gaussian distribution function, μ denotes the mathematical expectation, and σ 2 denotes the variance of random variable set κ 1 . Regarding that the actual human stiffness parameter K h1 can be deemed as μ, we can estimate K h1 according to the Bayesian parameter estimation method if σ 2 is known to the control designer. We rewrite the cost function (9) as follows: whereμ is the estimate of μ. We can obtain the predictor probability distribution of stiffness parameter p(μ) as follows: where μ 0 and σ 2 0 denote the predictor expectation and variance of μ, and their values can be found based on the literature about human stiffness measurement [77]. We can obtain the updater probability distribution p(μ|κ) as follows: where p(κ|μ) denotes the joint probability distribution, and it can be calculated as where (12) and (14) into (13), we can obtain the updater probability distribution p(μ|κ) as follows: where α is introduced to absorb the irrelevant terms about μ. Considering that p(κ|μ) and p(μ) follow the Gaussian distribution, we can rewrite (15) as where α 1 and α 2 are the parameters used to absorb the irrelevant items of μ. Note that p(μ|κ) follows the Gaussian distribution, so we can conclude that Because the coefficient in exponential term in (17) equals its counterpart in (16), we can obtain We can conclude that If we use the quadratic cost function like (9), the Bayesian parameter estimationμ can be described as the conditional expectation when κ is given and μ can be estimated as follows: Thus, the Bayesian estimation of μ can be rewritten aŝ From (22), we can conclude that the estimate of human stiffness parameterK h1 remains in the interval from (μ −σ ) to (μ +σ ), that is Then, we can obtain the corresponding human motion intention estimatex h1 as follows: SinceK h1 obeys the Gaussian distribution, the corresponding human motion intention estimatex h1 also obeys the Gaussian distribution, that isx where μ x and σ x are the expectation and the variance ofx h1 , respectively. They can be described as where x 1 denotes the position in one dimension.
Along with increasing n,σ converges to a small value, and μ converges toμ n . μ x converges to (f h1 /μ n ) + x 1 and σ x converges to zero. Using this method, we can estimate K h1 and x h1 in one dimension. In a similar way, human's stiffness matrixK h and motion intention vectorx h can be obtained by Bayesian parameter estimation.

IV. CONTROL DESIGN
Afterx h andK h are obtained, we set x d asx h , and set K d asK h to achieve the task objective. We set D d as where denotes a proper coefficient between 0 and 1. We set inertia matrix d close to the robot's inertia matrix M r . According to (6), we construct the error signal as where e = x d −x, and if we want to achieve the relationship in (6), we should make converge to zero. To facilitate analysis, we define another impedance error ω as where We choose two positive-definite matrices A and B as And we defineḟ According to (29) and (30), we rewrite (28) as Similar to [78], we define an auxiliary variable z as so we can rewrite (31) as When z converges to zero, we can conclude thatż → 0 if its limit exists. We define a virtual state variable vector x r aṡ so z can be rewritten as In the following text, we employ z to design an impedance controller and analyze control stability. Consider the following Lyapunov function candidate as: Differentiating V 1 with respect to time, we havė is skew-symmetric. Thus, we can rewritė V 1 asV Considering (35), we rewrite (3) as soV 1 can be written aṡ and the model-based impedance controller u can be designed as where K g is a positive-definite matrix, when u is designed as (41), we can obtaiṅ To address uncertainties in the robot's dynamic model, that is, M r (x), C r (x,ẋ), and G r (x) are unknown in practical situations, an adaptive impedance control is designed in this part. The adaptive law is designed aṡ whereŴ i is the weight estimate of NN, i = T i is a positive gain matrix, and δ i is a small positive constant which is used to improve the system robustness.
is the input of NN.Ŵ T S(Z) is used to estimate W * T S(Z) as follows: where W * i is the actual weight of NN, S(Z) denotes the basis function, and the estimation error (Z) stays in bounds over the compact set , ∀Z ∈ , || (Z)|| <¯ , with¯ as a positive constant.
The NN impedance controller can be designed as We consider another Lyapunov function V 2 as We define the weight errorW i =Ŵ i − W * i . Differentiating V 2 with respect to time, we havė Substituting (45) into (47), we can obtaiṅ Substituting (43) into (48), we havė We can obtaiṅ where where ς denotes the eigenvalue of a matrix, and¯ denotes the bound of . For ensuring ρ > 0, we should make ς min (K g − (1/2)I) > 0, ς max ( −1 i ) > 0. Theorem 1: For each compact set 0 , the initial conditions z 0 andŴ 0 are in bounds, the controller (45) guarantees that the closed-loop error signal z remains in the compact set z , and the weight errorW remains in the compact set W , that is where D = 2(V 2 (0) + C)/ρ with positive constants C and ρ is given in (51).

V. SIMULATION
In this section, we consider a scenario where a human partner is holding hand grasp on robotic end-effector with a force sensor. A two-link revolute joint robot shown in Fig. 3 is considered and an interaction force is applied to the end-effector by the human partner.
If M(q), C(q, q), G(q), and J are obtained, we can calculate robot's dynamic parameter matrices in the Cartesian space M r (x), C r (x,ẋ), and G r (x) in (4).
We consider that a human partner applies interaction force to the hand grasp on the end-effector from the initial position

A. Estimation of Human Stiffness and Motion Intention
We suppose that the human's real dynamic model in the X-direction can be described as f h1 = −300(x(1) − 0.75), where the actual human stiffness in the X-direction K h1 = 300 Nm, and human motion intention in the X-direction x d1 = 0.75 m. We use the Bayesian method to estimate human stiffness K h1 in the X-direction, and the same method is used for estimating K h2 in the Y-direction. First, we set a predictor probability distribution of human stiffness parameter p(μ 1 ) as follows: In the random variable set κ 1 that obeys the distribution κ 1 ∼ N(μ, 10 2 ), using the proposed method in Section III, different predictor probability distributions of human stiffness parameter K h1 can be estimated as shown in Fig. 4(a). From this figure, we can conclude that K h1 can be estimated with different mathematical expectations or different variances of p(μ 1 ). In Fig. 4(b), we can see that by setting different stiffness parameters 200 N/m, 300 N/m, and 400 N/m, respectively, K h1 can be successfully estimated by our proposed method in the same predictor probability distribution of human stiffness parameter. In Fig. 5, human motion intention estimationx d1 and variance ofx d1 have been obtained by dynamic relationship between x d1 and K h1 . Different human motion intentions have been set as 0.45 m, 0.75 m, and 0.90m when p(μ 1 ) ∼  N(200, 10 2 ). We can see that with the proposed method, different human motion intention can be estimated.

B. Impedance Control With Neural Networks
As discussed in Section II, we set the target impedance model as a simplified spring model   Fig. 6(b) shows the velocity and velocity error in the X-direction between x(2) andẋ d1 . Note that when there exists no interaction force, the position error and velocity error will converge to zero according to the dynamical relationship in (6). Fig. 7(a) and (b) shows the position and position error, and the velocity and velocity error in the Y-direction, respectively. Fig. 8(a) shows the tracking performance of velocity x (2) in the X-direction, and Fig. 8(b) shows the tracking performance of auxiliary variable z 1 in the X-direction. We can conclude that under the proposed method, the error signal converges to zero. Fig. 9 shows the interaction force between human and object f r1 in the X-direction.

VI. EXPERIMENT
In this section, we consider a scenario where an interaction force is applied to the arm of a robotic manipulator by a human partner. We use the S 0 shoulder joint on the right arm of dualarm humanoid robot Baxter in our experiment. A human robot   interactive experiment is developed to prove the validity of our proposed control method.

A. Experimental Settings
Baxter robot has torque sensors in every joint of both two arms. Angle, angle velocity, and torque information can be read from its dedicated controller. Seen from Fig. 10, there are two computers (A and B) for controlling robot and calculation in this experiment. Computer A is used to calculate the NN compensation by MATLAB Simulink and transform the compensation value to the computer B by UDP. Computer B is used to receive the robot state signals and generate control signal to control the robot by Baxter Robot Operating System SDK (RSDK) in Ubuntu 14.04 LTS. We rewrite the target impedance model in the joint space as τ f r = K S 0 (x − x d ), and we consider the human impedance model in the joint space as τ f h = K h (x − x h ). K S 0 and K h denote S 0 and human joint stiffness parameter, respectively, and x h denotes the human target angle, x denotes the current angle, and τ f h denotes the interaction torque.

B. Case 1: No Estimation
In this part, we consider a scenario that a human partner operates S 0 shoulder joint of Baxter robot's right arm to the human target angle. We design robot target impedance stiffness parameter K S 0 as 3 Nm/rad, but different human target angles x h : 0.2 and 0.8 rad. Fixed desired angle x d of the robot is considered in this experiment and the robot initial position is set as 0 rad. An interaction force is applied to the robot arm from 3 to 13 s. Seen from Fig. 11(a), the robot moves from 0 to 0.5 rad driven by human partner and back to 0 rad under the impedance control method. Fig. 12 shows that when K S 0 is fixed, the interaction torques have the proportional relationships with the error between current angle x and desired angle x d . Larger error between current position and x d will generate greater interaction torque.

C. Case 2: Motion Intention Estimation
Motion intention estimationx d is involved in this part. In Fig. 11(b), the robot moves from initial angle 0 rad to target angle of 0.5 rad driven by human partner, an interaction torque is applied to the robot arm from 3 to 8 s, and the robot will remain the current angle after 8 s when motion intention estimation based on the Bayesian estimated method is involved. Fig. 13 shows the relationships between the interaction torques and x − x d when human moves robot to 0.2 and 0.8 rad. As can be seen from Fig. 14(b) and (d), we can conclude that the interaction torque under our proposed method is smaller than the torque under impedance control when motion intention estimation is not involved. And the robot will remain in the position when interaction torque disappears as can be seen Fig. 14(c). In this part, we also utilize the NN method to estimate human motion intention for comparison with our proposed method. Indicated from Fig. 14(e), the convergence of the NN estimation method is slower than our proposed Bayesian estimation method. NNs rely on online sensor information which will bring heavy computational burden to influence convergence.

D. Case 3: Impedance Estimation
In this part, the target angle impedance stiffness value is set as 3 Nm/rad and 15 Nm/rad, respectively. The experimental process is same as the process in case 1. Indicated from Fig. 15, we can see the proportional relationships withK h and interaction torque, that is, larger stiffness will generate greater  interaction torque at the same angle displacement. We also consider the human stiffness estimation based on the Bayesian method in Fig. 16, from which we can conclude that the joint stiffness can be estimated by our proposed method. In Fig. 17, we set the predictor probability distribution of joint stiffness parameter as p(μ) ∼ N(1, 0.1 2 ), N(5, 0.1 2 ), respectively. Joint stiffness can be estimated successfully considering different probability distributions of human stiffness parameter.

E. Case 4: Simultaneous Estimations
In this part, we use the Bayesian method to estimate joint stiffness and human target angle simultaneously, where the predictor probability distribution of stiffness parameter p(μ) is set as p(μ) ∼ N(1, 0.1 2 ). The experimental process is divided into two phases. In the first phase, it is the same as the process (S 0 joint) in case 2, where the human partner moves the robot to the target position 1. In the second phase, we utilize the joint S 1 to lift the robot to the target position 2. The mean and standard deviation of the above measures are computed using 50 data points (10 human subjects × 5 repetitions). Each of ten human subjects (P1, P2, . . . , and P10) repeats the task for five times (T1, T2, T3, T4, and T5). Indicated from Fig. 18, we  can see that both human motion intention and joint stiffness can be estimated successfully, which show the robustness of the proposed method. We provide the statistical analysis of estimated stiffness of one human subject for five repetitions when interacting with the robot's S 0 and S 1 joints. Indicated  from Fig. 18(b) and (d), we can see that all estimated stiffness parameters converge to a constant value. Table I shows that the convergence values are "8.78 ± 0.13 Nm/rad" and "9.33 ± 0.11 Nm/rad" in S 0 and S 1 joints, respectively. In Table II, we can find that the stiffness of ten human subjects can be estimated successfully, and all estimated parameters converge to constant values in reasonable times. And the experimental results in a 7-degree-of-freedom are shown in Fig. 19, the proposed controller and Bayesian estimation method are utilized in this task. The experimental results on a Baxter robot platform illustrate good performance.

VII. CONCLUSION
In this article, a Bayesian method has been proposed to estimate human impedance and motion intention in a human-robot collaborative task. Estimated stiffness obeying the Gaussian distribution has been obtained by Bayesian estimation combining with the prior knowledge of human stiffness. According to the dynamic relationship, human motion intention can be also estimated. NNs have been used to compensate for uncertainties in robotic dynamics and an adaptive impedance control strategy has been employed to track a target impedance model. Comparative simulation and experimental results have been carried out to verify the advantages of the proposed control strategy and the effectiveness of estimation method.