Estimation of Involuntary Components of Human Arm Impedance in Multi-Joint Movements via Feedback Jerk Isolation

Stable and efficient coordination in physical human-robot interaction requires consideration of human feedback behavior. In unpredictable tasks, where voluntary cognitive feedback is too slow to guarantee desired task execution, the human must rely on involuntary intrinsic and reflexive feedback. The combined effects of these two feedback mechanisms and the inertial characteristics can be summarized in the involuntary impedance components. In this work, we present a method for the estimation of the involuntary impedance components of the human arm in multi-joint movements. We apply force perturbations to evoke feedback jerks that can be isolated using a high pass filter and limit the duration of the estimation interval to guarantee exclusion of voluntary cognitive feedback. Dynamic regressor representation of the rigid body dynamics of the arm and first order Taylor series expansion of the feedback behavior yield a model that is linear in the involuntary impedance components. The constant values of the inertial parameters are estimated in a static posture maintenance task and subsequently inserted to estimate the remaining components in a dynamic movement task. The method is validated with simulated data of a neuromechanical model of the human arm and its performance is compared to established methods from the literature. The results of the validation demonstrate superior estimation performance for moderate movement velocities, and less influence of the variability of the movements. The applicability to real data and the plausibility of the limited estimation interval are successfully demonstrated in an experiment with human participants.

1. Introduction

Technological advancements in robotics are enabling robotic assistance through physical interaction of human and robot for applications in medical, domestic, and industrial domains. In physical human-robot interaction (pHRI), control strategies are designed to provide efficient and intuitive interaction, during which instabilities must be avoided to guarantee safety and comfort of the human. Fulfillment of these requirements requires consideration of human feedback behavior in the control design process.

During execution of a desired motor task, joint torques produced by the neuromuscular system are composed of a feedforward and a feedback component (Tee et al., 2004). Due to a priori calculation, the feedforward component cannot account for unpredictable dynamics, which can result from incomplete or incorrect internal models (Tomi et al., 2008), inherent neural noise (Faisal et al., 2008), and unexpected external perturbations (Gomi and Kawato, 1997). The resulting deviations are compensated by restoring torques that contain effects of muscle intrinsic viscoelastic properties, reflexive feedback, and cognitive feedback (Gomi and Osu, 1998). Apart from the effects of the muscle intrinsic viscoelastic properties, these torques are produced at different task-dependent delays. As cognitive feedback possesses the longest delays [in the order of 100 ms (Franklin and Wolpert, 2011)], it may be too slow to guarantee desired task execution in unpredictable tasks (Mehta and Schaal, 2002). Thus, in such situations, the human must rely on intrinsic and reflexive feedback. In this work, all cognitive feedback at supraspinal level is referred to as voluntary feedback and the combined effects of all feedback mechanisms with shorter delays than voluntary feedback are grouped into involuntary feedback (see Figure 1).

FIGURE 1

Figure 1. Classification of intrinsic, reflexive, and cognitive feedback due to unpredictable dynamics into involuntary and voluntary feedback. All cognitive feedback at supraspinal level is referred to as voluntary feedback and the combined effects of all feedback mechanisms with shorter delays than voluntary feedback are grouped into involuntary feedback.

The restoring torques generated by involuntary feedback can be modeled by a linear system composed of the involuntary impedance components, damping and stiffness (Bennett et al., 1992). Due to length- and velocity-tension relationships of muscle fibers and tendons, both components depend on joint angles and angular velocities. Additionally, the central nervous system (CNS) is able to modulate both components through co-contraction of antagonistic pairs of muscles at the respective joints (Mussa-Ivaldi et al., 1985), e.g., to achieve accuracy requirements (Lametti et al., 2007) or to compensate environmental instabilities (Burdet et al., 2001). Analysis of such modulation strategies in experiments that emulate realistic pHRI may provide valuable information for the control design process, e.g., for stability assessment or control behavior adaptation. However, acquisition of this information requires an estimation method that is compatible with the short estimation interval, in which only involuntary feedback is active. It must enable precise execution of perturbations and accurate isolation of resulting feedback behavior, as estimation errors in the unperturbed states directly influence the estimation accuracy of the involuntary impedance components (Burdet et al., 2013). To the best of the author's knowledge, existing methods are either not able to guarantee exclusion of voluntary feedback (e.g., Gomi and Kawato, 1997; Erden and Billard, 2015) or are limited to application to free, unfettered movements (Piovesan et al., 2013).

In this work, we present a method for the estimation of the involuntary impedance components of the human arm in multi-joint movements. We apply force perturbations in order to evoke deviations during two-dimensional point to point arm movements. These perturbations are designed such that the dominant frequencies of the jerks of the evoked feedback behavior lie above those of the unperturbed movements. Thus, the feedback behavior can be isolated by a high pass filter. The duration of the estimation interval is limited to guarantee the exclusion of voluntary feedback. Dynamic regressor representation of the rigid body dynamics and first order Taylor series expansion of the feedback behavior yield a model that is linear in the involuntary impedance components in joint space. The constant values of the inertial parameters are estimated in a static posture maintenance task and subsequently used for the estimation of the remaining involuntary impedance components in a dynamic movement task. Both the feedback jerk isolation and the involuntary impedance estimation are validated with simulated data of a neuromechanical model of the human arm. We compare the validation results to those obtained by application of the methods presented in Gomi and Kawato (1997) and Erden and Billard (2015). In the validation of the feedback jerk isolation, we additionally analyze the effects of different movement velocities as well as different frequencies and amplitudes of neural noise (see Figure 2). Finally, we perform an extensive evaluation of the applicability of the method to real data based on an experiment with human participants. It includes an analysis of the effects of different durations of the estimation interval.

FIGURE 2

Figure 2. Schematic diagram of the major elements of the method and the validation with simulated data. The colored boxes, which are defined in the data acquisition part of the diagram, indicate which sets of simulated data are applied and analyzed in which phases of the validation.

1.1. Related Work

Estimation of multi-joint arm impedance is generally either performed by application of a perturbation paradigm, in which perturbations are used to evoke deviations from unperturbed states to measure variational dynamics, or by application of electromyographic (EMG) sensing modalities. In the latter category of approaches, multiple studies combine measurements of muscle activity with parametric muscle models, e.g., linear models of muscle stiffness (Osu and Gomi, 1999), quadratic models of muscle tension (Shin et al., 2009), or calibrated models of musculotendon unit forces (Buzzi et al., 2017). In Kim et al. (2009), an artificial neural network produces a mapping between EMG data and stiffness estimates, that are obtained from measured joint torques and an empirically determined linear model. In Lakatos et al. (2013) and Ajoudani et al. (2015), similar mappings are produced by pairing stiffness estimates obtained by application of the perturbation paradigm with EMG data. All of these EMG approaches have in common that they exclusively estimate stiffness.

Estimation of multi-joint arm impedance by application of the perturbation paradigm can generally be divided into two main categories: estimation in static posture maintenance tasks and in dynamic movement tasks. In both categories, perturbations are used to evoke deviations from unperturbed states and impedance is estimated based on the variational dynamics. Impedance estimation in static posture maintenance tasks is significantly less complex, as it does not require estimation of unperturbed states of an underlying movement. Furthermore, it allows for application of position perturbations, which are executed by moving the hand along a perturbation trajectory and thus enable a priori definition of the variational kinematics. Some studies exploit this advantage by using perturbation trajectories with plateau phases, during which the deviation remains constant (Mussa-Ivaldi et al., 1985; Dolan et al., 1993; Gomi and Osu, 1998; Darainy et al., 2004; Masia and Squeri, 2015). During the plateau phase, variational velocities and accelerations are zero. Thus, the respective variational forces can be used to exclusively estimate stiffness. Some studies use position perturbations defined by third (Artemiadis et al., 2010; Patel et al., 2014) or fifth degree polynomials (Lakatos et al., 2011, 2013), which are specifically designed to improve the conditioning of the estimation. Other studies use stochastic position (Perreault et al., 2004; Ajoudani et al., 2015) or force perturbations (Palazzolo et al., 2007) in combination with significantly longer estimation intervals to estimate frequency domain impedance transfer functions.

Impedance estimation during dynamic movement tasks is significantly more challenging than during static posture maintenance tasks, as it requires precise estimation of the unperturbed states of the underlying movement. In Wang et al. (2001), mechanical impedance is solely estimated in terms of feedback forces. In Hondori and Tech (2011), it is represented by the ratio of frequency domain forces and velocities. In Burdet et al. (2000), position perturbations with plateau phases are used to exclusively estimate stiffness in a 60 ms interval that starts 120 ms after perturbation onset. In Darainy et al. (2007) and Wong et al. (2009b), similar methods are used to exclusively estimate stiffness in a 50 ms interval that starts 250 ms after perturbation onset. In Piovesan et al. (2013), force perturbations and time-frequency analysis are used to estimate impedance 135 ms after perturbation onset. Due to reliance on vibrational energy, this method is only applicable to free, unfettered movements. In Gomi and Kawato (1997), force perturbations are used to estimate impedance by least squares estimation in a 280 ms interval. Erden and Billard (2015) use a similar method to estimate impedance in a 250 ms interval.

To the best of the author's knowledge, among existing methods in the literature, only those presented in Gomi and Kawato (1997) and Erden and Billard (2015) are able to estimate stiffness, damping, and inertial characteristics during multi-joint arm movements that emulate realistic pHRI. In Gomi and Kawato (1997), participants perform transversal and longitudinal point to point movements. The movements are guided by a target position that moves along a reference trajectory and is displayed on a computer monitor. First order Taylor series expansion is used to obtain a linearized model of the rigid body dynamics and the feedback behavior in joint space. Inertial parameters are estimated a priori and used for the estimation of damping and stiffness. Variational dynamics are obtained by calculating the difference to the averaged unbiased dynamics of all perturbed movements, which are calculated by subtracting the offsets at perturbation onset. In Erden and Billard (2015), participants perform a manual tungsten inert gas welding task in cooperation with a robot. The movements are guided by a straight reference trajectory and the participants are instructed “to do their best to achieve a good welding.” Variational behavior is described by a linear model of diagonal inertia, damping, and stiffness matrices in Cartesian space, which are estimated simultaneously. Variational positions and forces are obtained by subtraction of the offsets at perturbation onset and variational velocities and accelerations are obtained by differentiation. Both Gomi and Kawato (1997) and Erden and Billard (2015) are unable to guarantee exclusion of voluntary feedback, as both methods require estimation intervals that are significantly longer than the delay of voluntary feedback. In Gomi and Kawato (1997), this limitation is compensated by application of the do-not-intervene-voluntarily paradigm, which means that participants are instructed to not intervene voluntarily in response to the perturbation. However, this approach poses two problems. Firstly, despite these instructions, voluntary intervention is nonetheless possible and exclusion of voluntary feedback cannot be guaranteed. Secondly, instructing the participants to not intervene voluntarily substantially limits the number of possible, realistic pHRI application scenarios.

A preliminary work has been presented in Börner et al. (2018). In this preliminary work, the rigid body dynamics were linearized in joint space and the involuntary impedance components were estimated in Cartesian space. Due to the linearization in joint space, assumptions concerning negligibility of inertia and Coriolis terms were necessary to obtain a linear second order system in Cartesian space. The elements of the Cartesian endpoint inertia were calculated with anthropometric data and subsequently inserted to estimate the remaining involuntary impedance components in a dynamic movement task. Furthermore, the application of a preliminary perturbation design required assumptions regarding correlations of the beginning of the estimation interval and noticeable deviations due to the perturbation.

The remainder of this article is structured as follows: section 2 contains the materials and methods, which include the derivation of the involuntary impedance model as well as the presentation of the feedback jerk isolation, the perturbation design, and the involuntary impedance estimation. Furthermore, the details of the simulation and the experiment are presented. Section 3 contains the results, which include the validation and the comparison with existing methods based on simulated data as well as the evaluation with experimental data. Section 4 contains the discussion and conclusions are summarized in section 5.

2. Materials and Methods

We begin this section by formulating the problem considered in this work. Subsequently, we present the feedback jerk isolation and the involuntary impedance estimation. Finally, we introduce the details of the simulation and the experiment.

2.1. Problem Formulation

In this section, we first derive a general human motor behavior model, which then serves as the basis for the derivation of a model of the involuntary impedance components.

2.1.1. Human Motor Behavior Model

The human arm is modeled as a multi-joint two-link system. In order to reduce complexity and neglect gravity, movements are constrained to the horizontal plane. The rigid body dynamics is given by

$\begin{array}{ll}{\tau }_{\text{int}}+{\tau }_{\text{ext}}=M_q(q)\ddot{q}+C(q,\dot{q})\dot{q},& (1)\end{array}$

where q is the 2 degree of freedom (DoF) arm configuration in joint space, M_q is the symmetric and positive-definite inertia matrix (Hogan, 1985), C is the Coriolis/centrifugal matrix, τ_int are the internally generated joint torques, and τ_ext are the external torques (Gomi and Kawato, 1997). The arm configuration $q={[q_1,q_2]}^{\text{T}}$ is defined by the angles of the shoulder q₁ and the elbow q₂.

The internal torques τ_int are produced by the muscle tensions m that act on the musculoskeletal system:

$\begin{array}{ll}{\tau }_{\text{int}}=J_m^{\text{T}}(\lambda )m(\lambda ,\dot{\lambda },a),& (2)\end{array}$

where λ are the muscle lengths, $\dot{\lambda }$ are the respective derivatives, a are the muscle activations, and J_m is the muscle Jacobian that contains the muscle moment arms that are necessary for the transformation to joint space (Shin et al., 2009). During execution of a motor task, the muscle activations a consist of a feedforward term a_FF, a feedback term a_FB, and a neural noise term a_N (Franklin et al., 2008):

$\begin{array}{ll}a=a_{\text{FF}}(\theta )+a_{\text{FB}}+a_{\text{N}}(a).& (3)\end{array}$

The feedforward term a_FF is calculated by the CNS through a priori optimization of the statistics of the desired motor behavior with respect to costs defined by the task-specific input parameters θ, which depend on factors, such as environmental constraints and task requirements (Harris and Wolpert, 1998). The neural noise term a_N represents signal-dependent noise whose variance increases with the magnitude of the muscle activations a (Slifkin and Newell, 1999). Deviations caused by unpredictable dynamics are compensated by the feedback term a_FB, which consists of multiple components that depend on delayed afferent sensory information and are produced at different delays (Franklin and Wolpert, 2011):

$\begin{array}{ll}{a}_{\text{FB}}=\{\begin{array}{ll}0\hfill & \forall t_{\text{p}}\in \left[0,{\delta }_{\text{r}}\right]\hfill \\ a_{\text{FB},\text{r}}(\lambda ,\dot{\lambda })\hfill & \forall t_{\text{p}}\in \left]{\delta }_{\text{r}},{\delta }_{\text{v}}\right]\hfill \\ a_{\text{FB},\text{r}}(\lambda ,\dot{\lambda })+a_{\text{FB},\text{v}}(\theta )\hfill & \forall t_{\text{p}}>{\delta }_{\text{v}}\hfill \end{array},& (4)\end{array}$

where a_FB,r and a_FB,v are reflexive and voluntary feedback muscle activations, respectively. The variables δ_r and δ_v are the corresponding delays and t_p = t − t₀ is the time after the onset of the unpredictable dynamics at t₀. In the interval ]δ_r, δ_v], the feedback muscle activations a_FB only consist of reflexive feedback muscle activations a_FB,r, which are affected by neural conduction delays that depend on the length and type of the nerve fiber. The fastest reflexive feedback is produced by the short-latency monosynaptic stretch reflex, with a delay δ_r,s in the order of 10 − 40 ms (Matthews, 1991). Slightly slower reflexive feedback is produced by the cortical component of the long-latency stretch reflex, with a delay δ_r,l in the order of 30 − 70 ms (Thilmann et al., 1991). As the feedback behavior in this interval does not recruit higher-order cognitive processes, the reflexive feedback muscle activations a_FB,r only depend on the muscle lengths λ and the respective derivatives $\dot{\lambda }$. For all t_p > δ_v, the voluntary feedback muscle activations a_FB,v, which are defined by the task-specific input parameters θ, additionally contribute to the feedback muscle activations a_FB. For voluntary feedback in response to haptic motion perception, the transmission of the associated proprioceptive sensory information to the CNS is subject to delays in the order of 100 ms. Conduction delays of descending motor commands from the motor cortex to the arm muscles are ~ 15 ms (Merton and Morton, 1980). Consequently, the minimum delay δ_v of the voluntary feedback muscle activations a_FB,v is in the order of 115 ms.

The internal torques τ_int are composed of a feedforward term τ_FF, a feedback term τ_FB, and a neural noise term τ_N, that represent the equivalents to the respective muscle activations a in (3):

$\begin{array}{ll}{\tau }_{\text{int}}={\tau }_{\text{FF}}(\theta )+{\tau }_{\text{FB}}+{\tau }_{\text{N}}(a).& (5)\end{array}$

Analogous to the feedback muscle activations a_FB in (4), the feedback torques τ_FB consist of multiple components that are produced at different delays (Lakatos et al., 2011):

$\begin{array}{ll}{\tau }_{\text{FB}}=\{\begin{array}{ll}{\tau }_{\text{FB},\text{i}}\hfill & \forall t_{\text{p}}\in \left[0,{\delta }_{\text{r}}\right]\hfill \\ {\tau }_{\text{FB},\text{i}}+{\tau }_{\text{FB},\text{r}}\hfill & \forall t_{\text{p}}\in \left]{\delta }_{\text{r}},{\delta }_{\text{v}}\right]\hfill \\ {\tau }_{\text{FB},\text{i}}+{\tau }_{\text{FB},\text{r}}+{\tau }_{\text{FB},\text{v}}\hfill & \forall t_{\text{p}}>{\delta }_{\text{v}}\hfill \end{array},& (6)\end{array}$

where τ_FB,i, τ_FB,r, and τ_FB,v are intrinsic, reflexive, and voluntary feedback torques, respectively. In the interval [0, δ_r], the feedback torques τ_FB only consist of intrinsic feedback torques τ_FB,i, which result from muscle intrinsic viscoelastic properties and are thus not affected by any delays (Tee et al., 2004). In the subsequent interval ]δ_r, δ_v] and for all t_p > δ_v, the feedback torques τ_FB additionally contain contributions from reflexive feedback torques τ_FB,r and voluntary feedback torques τ_FB,v, respectively. Inserting (6) in (5) yields the internal torques τ_int for all t > δ_v:

$\begin{array}{ll}{\tau }_{\text{int}}={\tau }_{\text{FF}}(\theta )+{\tau }_{\text{FB},\text{i}}+{\tau }_{\text{FB},\text{r}}+{\tau }_{\text{FB},\text{v}}+{\tau }_{\text{N}}(a).& (7)\end{array}$

This general model represents the basis for the derivation of the involuntary impedance model.

2.1.2. Dynamic Regressor Representation

The rigid body dynamics (1) is non-linear with respect to the arm configuration q and the respective derivatives $\dot{q}$ and $\ddot{q}$. Nonetheless, the Lagrangian formalism enables the derivation of the dynamic regressor representation, which is linear with respect to the standard inertial parameters of the system:

$\begin{array}{ll}{\tau }_{\text{int}}+{\tau }_{\text{ext}}=Y(q,\dot{q},\ddot{q})\pi ,& (8)\end{array}$

$\begin{array}{ll}Y(q,\dot{q},\ddot{q})=\frac{\partial (M_q(q)\ddot{q}+C(q,\dot{q})\dot{q})}{\partial \pi },& (9)\end{array}$

where Y is the regressor matrix and $\pi ={[{\pi }_1^{\text{T}},{\pi }_2^{\text{T}}]}^{\text{T}}$ is the standard inertial parameter vector, which is composed of the inertial parameters of the upper arm π₁ and the forearm π₂ (Gautier and Khalil, 1988). In general, the inertial parameters consist of the masses and the mass moments of first and second order (Khalil and Dombre, 2004). However, some parameters do not have any effect on the system dynamics which corresponds to a zero column in the regressor matrix Y. For the multi-joint arm, omission of these parameters yields the reduced inertial parameter vector

$\begin{array}{ll}{\pi }_{\text{r}}={\left[I_{x_2{x}_2,1},I_{x_2{x}_2,2},m_2{l}_{c{x}_3,2},m_2\right]}^{\text{T}},& (10)\end{array}$

in which, according to the parallel axis theorem,

$\begin{array}{ll}{I}_{x_2{x}_2,i}={\tilde{I}}_{x_2{x}_2,i}+m_i{l}_{c{x}_3,i}^2.& (11)\end{array}$

Thus, the inertial properties are defined by the moments of inertia Ĩ_{x2x2, i}, the masses m_i, and the lengths from the joints to the centers of gravity of the links l_{cx3, i}. Comprehensive analysis of the reduced regressor matrix Y_r reveals that further parameter reduction is possible. Appropriate transformations yield the base inertial parameter (BIP) vector $\overline{\pi }$, which contains the minimal set of parameters, and the corresponding base regressor matrix $\overline{Y}$ of the multi-joint arm:

$\begin{array}{ll}\overline{\pi }={\left[I_{x_2{x}_2,1}+m_2{l}_1^2,I_{x_2{x}_2,2},m_2{l}_{c{x}_3,2}{l}_1\right]}^{\text{T}},& (12)\end{array}$

$\begin{array}{ll}\overline{Y}(q,\dot{q},\ddot{q})=\left[\begin{array}{ccc}{\ddot{q}}_1& {\ddot{q}}_1+{\ddot{q}}_2& {\overline{y}}_{1,3}\\ 0& {\ddot{q}}_1+{\ddot{q}}_2& {\overline{y}}_{2,3}\end{array}\right],& (13)\end{array}$

with

$\begin{array}{l}{\overline{y}}_{1,3}=(2{\ddot{q}}_1+{\ddot{q}}_2)\text{cos}(q_2)-(2{\dot{q}}_1{\dot{q}}_2+{\dot{q}}_2^2)\text{sin}(q_2),\\ {\overline{y}}_{2,3}={\ddot{q}}_1\text{cos}(q_2)+{\dot{q}}_1^2\text{sin}(q_2),\end{array}$

where l₁ is the length of the upper arm (Klodmann et al., 2015).

2.1.3. Involuntary Impedance Model

While the sum of the intrinsic feedback torques τ_FB,i and the reflexive feedback torques τ_FB,r can be modeled by a linear system composed of the involuntary impedance components, damping and stiffness (Bennett et al., 1992), the voluntary feedback torques τ_FB,v are highly task-specific and can therefore not be modeled by such a general formulation (Franklin and Wolpert, 2011). Therefore, in this work, in order to guarantee exclusion of voluntary feedback, we limit the duration of the estimation interval to T_est = δ_v. Due to the short latency of the reflexive feedback and the variability of the respective delays, separation of the intrinsic and the reflexive feedback contributions is difficult (Burdet et al., 2013). Thus, we summate the contributions in the involuntary feedback torques

$\begin{array}{ll}{\overline{\tau }}_{\text{FB}}=\{\begin{array}{ll}{\tau }_{\text{FB},\text{i}}\hfill & \forall t_{\text{p}}\in \left[0,{\delta }_{\text{r}}\right]\hfill \\ {\tau }_{\text{FB},\text{i}}+{\tau }_{\text{FB},\text{r}}\hfill & \forall t_{\text{p}}\in \left]{\delta }_{\text{r}},T_{\text{est}}\right]\hfill \end{array}.& (14)\end{array}$

For small deviations, the feedback behavior evoked by the force perturbations can be described by a linearized model obtained by first order Taylor series expansion of the general model (7) about the unperturbed states q^*, ${\dot{q}}^{*}$, ${\ddot{q}}^{*}$, a^*, ${\tau }_{\text{int}}^{*}$, and ${\tau }_{\text{ext}}^{*}$. According to (2), the variational internal torques $\Delta {\tau }_{\text{int}}={\tau }_{\text{int}}^{*}-{\tau }_{\text{int}}$ after the onset of the perturbation depend on the muscle lengths λ, the respective derivatives $\dot{\lambda }$, and the muscle activations a. Due to the confinement to the estimation interval [0, T_est], the variational muscle activations Δa = a^* − a, as defined by (3) and (4), can only consist of reflexive feedback muscle activations a_FB,r. Considering $\Delta a=f(\lambda ,\dot{\lambda })$ as well as λ = J_m(λ)q and $\dot{\lambda }=J_m(\lambda )\dot{q}$, linearization of the internal torques τ_int for t_p ∈ [0, T_est] yields

$\begin{array}{ll}\Delta {\tau }_{\text{int}}=\frac{\text{d}{\tau }_{\text{int}}}{\text{d}\dot{q}}\Delta \dot{q}+\frac{\text{d}{\tau }_{\text{int}}}{\text{d}q}\Delta q,& (15)\end{array}$

in which all variational variables, indicated by a Δ symbol, represent the deviations from the unperturbed states, e.g., Δq = q^* − q, and the total derivatives are defined as

$\begin{array}{ll}\frac{\text{d}{\tau }_{\text{int}}}{\text{d}\dot{q}}=\frac{\partial {\tau }_{\text{int}}}{\partial \dot{q}}+\frac{\partial {\tau }_{\text{int}}}{\partial a}\frac{\partial a}{\partial \dot{q}},\frac{\text{d}{\tau }_{\text{int}}}{\text{d}q}=\frac{\partial {\tau }_{\text{int}}}{\partial q}+\frac{\partial {\tau }_{\text{int}}}{\partial a}\frac{\partial a}{\partial q}.& (16)\end{array}$

Inserting (7) for t_p ∈ [0, T_est] and (14) in (15) yields

$\begin{array}{ll}\Delta {\tau }_{\text{int}}=-D_q(q,\dot{q},a)\Delta \dot{q}-K_q(q,\dot{q},a)\Delta q.& (17)\end{array}$

In this linearized model, the involuntary joint damping D_q and the involuntary joint stiffness K_q, which are both part of the involuntary impedance components¹, are defined as

$\begin{array}{ll}{D}_q(q,\dot{q},a)=-\frac{\text{d}{\overline{\tau }}_{\text{FB}}}{\text{d}\dot{q}},K_q(q,\dot{q},a)=-\frac{\text{d}{\overline{\tau }}_{\text{FB}}}{\text{d}q}.& (18)\end{array}$

Inserting the BIP vector $\overline{\pi }$ of (12) and the base regressor matrix $\overline{Y}$ of (13) in (8) and using the resulting linear model to calculate the variational internal torques Δτ_int yields

$\begin{array}{l}\Delta {\tau }_{\text{int}}=(\overline{Y}(q^{*},{\dot{q}}^{*},{\ddot{q}}^{*})-\overline{Y}(q,\dot{q},\ddot{q}))\overline{\pi }-\Delta {\tau }_{\text{ext}}\\ =\Delta \overline{Y}(q,q^{*},\dot{q},{\dot{q}}^{*},\ddot{q},{\ddot{q}}^{*})\overline{\pi }-\Delta {\tau }_{\text{ext}}.& (19)\end{array}$

Subsequently inserting (17) in (19) yields the involuntary impedance model

$\begin{array}{l}\Delta {\tau }_{\text{ext}}=\Delta \overline{Y}(q,q^{*},\dot{q},{\dot{q}}^{*},\ddot{q},{\ddot{q}}^{*})\overline{\pi }\\ +D_q(q,\dot{q},a)\Delta \dot{q}+K_q(q,\dot{q},a)\Delta q.& (20)\end{array}$

In some pHRI scenarios, measurements are performed in Cartesian space. In these scenarios, the arm configuration q can be obtained from the arm endpoint configuration x with the inverse kinematics

$\begin{array}{ll}{q}_1=\text{ta}{\text{n}}^{-1}\left(\frac{x_2}{x_1}\right)-\text{ta}{\text{n}}^{-1}\left(\frac{l_2\text{sin}(q_2)}{l_1+l_2\text{cos}(q_2)}\right),& (21)\end{array}$

$\begin{array}{ll}{q}_2=\text{co}{\text{s}}^{-1}\left(\frac{x_1^2+x_2^2-l_1^2-l_2^2}{2l_1{l}_2}\right),& (22)\end{array}$

where l₂ is the length of the forearm and the Cartesian origin is located in the shoulder joint. Furthermore, the Jacobian J can be used to transform the external endpoint forces u_ext to the external torques

$\begin{array}{ll}{\tau }_{\text{ext}}=J^{\text{T}}(q)u_{\text{ext}}.& (23)\end{array}$

The problem considered in this work consists of the estimation of the involuntary impedance components, i.e., BIP vector $\overline{\pi }$, damping D_q, and stiffness K_q, in the estimation interval [0, T_est]. This is to be achieved given the perturbed observations {x, u_ext}, which result from application of force perturbations during multi-joint arm movements, and requires estimation of the unperturbed dynamics $\left\{q^{*},{\dot{q}}^{*},{\ddot{q}}^{*},{\tau }_{\text{ext}}^{*}\right\}$.

2.2. Involuntary Impedance Estimation

In this section, we first present the estimation of the unperturbed dynamics via feedback jerk isolation and the perturbation design. Subsequently, we present the involuntary impedance estimation.

2.2.1. Feedback Jerk Isolation

According to the minimum jerk principle (Flash and Hogan, 1985), the CNS optimizes the arm endpoint trajectory in a point to point movement through minimization of the total endpoint jerk. We take advantage of this by designing the perturbation such that the dominant frequencies of the jerks of the evoked feedback behavior lie above those of the unperturbed movements. Due to this design, we are able to estimate the evoked feedback behavior in the form of the variational jerks $\Delta \stackrel{.}{x}$ through application of a high pass filter to the jerks $\stackrel{.}{x}$ of the perturbed movement. In order to achieve maximum pass band flatness and fast stop band roll-off, we perform this feedback jerk isolation with a Butterworth filter that is implemented as a digital biquad filter with filter order n_HP = 10. It is applied bi-directionally for zero phase distortion and the cut-off frequencies f_c,HP are defined based on the energy spectral densities (ESDs) ψ of the jerks $\stackrel{.}{x}$ of the unperturbed and perturbed movements (Stein, 2000):

$\begin{array}{ll}\psi (f)={\displaystyle |{\displaystyle {\int }_{-\infty }^{\infty }}{e}^{-2\pi ift}\stackrel{⃛}{x}(t)dt{|}^2.}& (24)\end{array}$

We use the set Ψ, which contains the ESDs ψ of the jerks $\stackrel{.}{x}$ of all unperturbed and perturbed movements, to calculate the averaged ESDs of the unperturbed movements ψ_UP, the perturbed movements ψ_P, and the resulting feedback behavior ψ_FB = ψ_P − ψ_UP. The cut-off frequencies f_c,HP are defined as the frequencies f at which ψ_FB > ψ_UP, i.e., the averaged ESDs of the feedback behavior ψ_FB become larger than those of the unperturbed movements ψ_UP. The high pass filtered jerks ${\stackrel{⃛}{x}}_{\text{HP}}$ yield the estimated variational jerks $\Delta \stackrel{⃛}{\widehat{x}}$ and the estimated unperturbed jerks ${\stackrel{⃛}{\widehat{x}}}^{*}$. Integration then yields the estimated variational kinematics $\{\Delta \widehat{x},\Delta \dot{\widehat{x}},\Delta \ddot{\widehat{x}}stretchy=\text{'}false\text{'}\}$ and the estimated unperturbed kinematics $\{{\widehat{x}}^{*},{\dot{\widehat{x}}}^{*},{\ddot{\widehat{x}}}^{*}\}$. For simplicity, from this point on, we refer to the averaged ESDs ψ_UP, ψ_P, and ψ_FB simply as ESDs.

The apparatus in our experiments is controlled by means of an admittance control scheme. Therefore, the external forces u_ext can be calculated with

$\begin{array}{ll}{u}_{\text{ext}}=u_{\text{pert}}-u_{\text{adm}}-M_{\text{handle}}\ddot{x},& (25)\end{array}$

$\begin{array}{ll}{u}_{\text{adm}}=M_{\text{adm}}\ddot{x}+D_{\text{adm}}\dot{x},& (26)\end{array}$

where u_pert is the perturbation force, u_adm is the admittance force, M_handle is the handle inertia, and M_adm and D_adm are the admittance inertia and damping, respectively. Inserting the estimated unperturbed kinematics ${\dot{\widehat{x}}}^{*}$ and ${\ddot{\widehat{x}}}^{*}$ in (25) and (26) yields the estimated unperturbed external forces ${û}_{\text{ext}}^{*}$ and the estimated variational external forces $\Delta {\widehat{u}}_{\text{ext}}$. The estimated unperturbed arm configuration ${\widehat{q}}^{*}$ and the respective derivatives ${\dot{\widehat{q}}}^{*}$ and ${\ddot{\widehat{q}}}^{*}$ are calculated using the inverse kinematics (21) and (22) and the estimated unperturbed external torques ${\widehat{\tau }}_{\text{ext}}^{*}$ are calculated using the transformation (23). Finally, the estimated unperturbed dynamics $\{{\widehat{q}}^{*},{\dot{\widehat{q}}}^{*},{\ddot{\widehat{q}}}^{*},{\widehat{\tau }}_{\text{ext}}^{*}\}$ yield the estimated variational dynamics $\{\Delta \widehat{q},\Delta \dot{\widehat{q}},\Delta \ddot{\widehat{q}},\Delta {\widehat{\tau }}_{\text{ext}}\}$. Figure 3 presents a block diagram that illustrates the complete procedure of the feedback jerk isolation.

FIGURE 3

Figure 3. Block diagram of the feedback jerk isolation. For every variational variable Δχ, the respective unperturbed variable χ^* = Δχ + χ is calculated with the respective perturbed variable χ, and vice versa.

2.2.2. Perturbation Design

The perturbation is designed to meet two essential criteria. First, the dominant frequencies of the jerks of the evoked feedback behavior should lie above those of the unperturbed movements. Second, in order to minimize movement interference and corrective oscillations, it should not only produce deviations, but also move the hand back toward the unperturbed states. In order to fulfill these criteria, we design the perturbation acceleration profile ẍ_pert through concatenation of two sinusoidal functions

$\begin{array}{ll}{\xi }_{\text{p}}=\{\begin{array}{ll}\frac{1}{2}{A}_{\text{p}}\left({\xi }_{\text{p},\text{sin}}+1\right)\hfill & \forall t_{\text{p}}\in [0,\frac{T_{\text{p}}}{2}]\hfill \\ -\frac{1}{2}{A}_{\text{p}}\left({\xi }_{\text{p},\text{sin}}+1\right)\hfill & \forall t_{\text{p}}\in [\frac{T_{\text{p}}}{2},T_{\text{p}}]\hfill \end{array}& (27)\end{array}$

with

$\begin{array}{ll}{\xi }_{\text{p},\text{sin}}=sin\left(\frac{4\pi t_{\text{p}}}{T_{\text{p}}}+\frac{3\pi }{2}\right),& (28)\end{array}$

where A_p and T_p are amplitude and duration. The first of the two functions ξ_p,1 with amplitude A_p,1 and duration T_p,1 is responsible for the deviation from the unperturbed states. The second function ξ_p,2 with amplitude A_p,2 and duration T_p,2 supplies the retracting movement toward the unperturbed states. In order to minimize hardware oscillations due to the perturbation, we define T_p,2 = (3/2)T_p,1 and scale A_p,2 such that all perturbation profiles are zero at the end of the perturbation. Lastly, we combine the perturbation acceleration profile ẍ_pert and velocity profile ẋ_pert with the admittance inertia M_adm and damping D_adm of our apparatus to obtain the perturbation force profile u_pert. It is defined by the amplitude of its first peak A_pert and its duration T_pert = T_p,1 + T_p,2. In the dynamic movement task in this work, A_pert = 40 N and T_p,1 = 70 ms result in T_pert = 175 ms. The resulting normalized perturbation force profile u_pert and the corresponding kinematics profiles $\{x_{\text{pert}},{\stackrel{.}{x}}_{\text{pert}},{\ddot{x}}_{\text{pert}},{\stackrel{.}{x}}_{\text{pert}}\}$ are illustrated in Figure 4. Given the normalized perturbation force profile u_pert, the perturbation force

$\begin{array}{ll}{u}_{\text{pert}}=A_{\text{pert}}{u}_{\text{pert}}{\left[cos{\varphi }_{\text{pert}},sin{\varphi }_{\text{pert}}\right]}^{\text{T}},& (29)\end{array}$

in which the perturbation angle ϕ_pert is defined by the set of perturbation angles Φ_pert.

FIGURE 4

Figure 4. Normalized perturbation force profile u_pert and kinematics profiles $\{x_{\text{pert}},{\stackrel{.}{x}}_{\text{pert}},{\ddot{x}}_{\text{pert}},{\stackrel{⃛}{x}}_{\text{pert}}\}$ of the dynamic movement task in the simulation and the experiment. The first part of the perturbation with duration T_p,1 is responsible for the deviation from the unperturbed states and the second part of the perturbation with duration T_p,2 supplies the retracting movement toward the unperturbed states.

2.2.3. Involuntary Impedance Estimation

Due to the limited duration of the estimation interval T_est, we assume that the involuntary impedance components D_q and K_q are constant for t_p ∈ [0, T_est]. Concatenation of the elements of the BIP vector $\overline{\pi }$ and the resulting involuntary impedance parameters ${\overline{D}}_q$ and ${\overline{K}}_q$ in the unknown model parameter vector

$\begin{array}{r}\zeta =[{\overline{\pi }}_1,{\overline{\pi }}_2,{\overline{\pi }}_3,{\overline{D}}_{q,11},{\overline{D}}_{q,12},{\overline{D}}_{q,21},{\overline{D}}_{q,22},\\ {{\overline{K}}_{q,11},{\overline{K}}_{q,12},{\overline{K}}_{q,21},{\overline{K}}_{q,22}]}^{\text{T}}& (30)\end{array}$

allows expression of (20) by the identification model

$\begin{array}{ll}A\zeta =b,& (31)\end{array}$

where A is the observation matrix and b is the output vector. For an interval with a total of N samples

$\begin{array}{ll}A={\left[X^{\text{T}}(1),X^{\text{T}}(2),\dots ,X^{\text{T}}(N)\right]}^{\text{T}},& (32)\end{array}$

$\begin{array}{ll}b={\left[\Delta {\widehat{\tau }}_{\text{ext}}^{\text{T}}(1),\Delta {\widehat{\tau }}_{\text{ext}}^{\text{T}}(2),\dots ,\Delta {\widehat{\tau }}_{\text{ext}}^{\text{T}}(N)\right]}^{\text{T}},& (33)\end{array}$

in which the independent variable matrix X is defined as

$\begin{array}{ll}X=\frac{\partial (\Delta \overline{Y}(q,{\widehat{q}}^{*},\dot{q},{\dot{\widehat{q}}}^{*},\ddot{q},{\ddot{\widehat{q}}}^{*})\widehat{\overline{\pi }}+{\overline{D}}_q\Delta \dot{\widehat{q}}+{\overline{K}}_q\Delta \widehat{q})}{\partial \zeta }.& (34)\end{array}$

The estimated model parameter vector $\widehat{\zeta }$ is given by

$\begin{array}{ll}\widehat{\zeta }={(A^{\text{T}}A)}^{-1}{A}^{\text{T}}b.& (35)\end{array}$

Due to the limited duration of the estimation interval T_est and the likewise limited duration of the perturbation T_pert, the perturbation in (29) does not possess sufficient richness of frequency components to guarantee persistent excitation (Söderström and Stoica, 1989). A common approach for the compensation of such effects is the a priori estimation of the inertial parameters (Gomi and Kawato, 1997; Darainy et al., 2007; Wong et al., 2009a,b) which only marginally influences the estimated impedance parameters (Gomi and Osu, 1998). As the inertial parameters possess constant values that are independent of the dynamics, the elements of the BIP vector $\overline{\pi }$ can be estimated a priori in a static postural task.

The intrinsic feedback behavior of the multi-joint arm is governed by spring-like characteristics that result from the elastic properties of the individual muscles. Thus, the intrinsic feedback forces possess zero curl and the intrinsic stiffness is symmetric (Shadmehr and Arbib, 1992). The reflexive feedback forces may possess non-zero curl which can only result from heteronymous inter-muscular reflex arcs with unequal activation gains (Hogan, 1985). However, as the resulting antisymmetric stiffness elements are significantly smaller than the corresponding symmetric ones, the reflexive feedback behavior is nonetheless governed by predominantly spring-like characteristics (Mussa-Ivaldi et al., 1985; Artemiadis et al., 2010). As the estimation interval in this work is significantly shorter than the ones in Mussa-Ivaldi et al. (1985) and Artemiadis et al. (2010), we assume that the stiffness ${\overline{K}}_q$ is symmetric for t_p ∈ [0, T_est]. As the potential energy of the vector field of the involuntary feedback forces must increase in response to deviations from the unperturbed states, the symmetric stiffness ${\overline{K}}_q$ must also be positive definite. In order to incorporate both symmetry and positive definiteness of the stiffness ${\overline{K}}_q$ in the estimation, we apply the Cholesky decomposition ${\overline{K}}_q={\overline{L}}_q{\overline{L}}_q^{\text{T}}$, where ${\overline{L}}_q$ is a lower triangular matrix (Horn and Johnson, 2012). With this decomposition and the a priori determined values of the estimated BIP vector $\widehat{\overline{\pi }}$, the unknown model parameter vector ζ reduces to

$\begin{array}{ll}\overline{\zeta }={\left[{\overline{D}}_{q,11},{\overline{D}}_{q,12},{\overline{D}}_{q,21},{\overline{D}}_{q,22},{\overline{L}}_{q,11},{\overline{L}}_{q,12},{\overline{L}}_{q,22}\right]}^{\text{T}}& (36)\end{array}$

and the elements of the output vector b are replaced by

$\begin{array}{ll}\overline{b}=\Delta {\widehat{\tau }}_{\text{ext}}-\Delta \overline{Y}(q,{\widehat{q}}^{*},\dot{q},{\dot{\widehat{q}}}^{*},\ddot{q},{\ddot{\widehat{q}}}^{*})\widehat{\overline{\pi }}.& (37)\end{array}$

Due to the non-linearity introduced by ${\overline{K}}_q={\overline{L}}_q{\overline{L}}_q^{\text{T}}$, we can no longer use the identification model in (31). Instead, we apply non-linear least squares analysis to minimize ${\Vert f(\overline{\zeta })\Vert }^2$ with

$\begin{array}{ll}f(\overline{\zeta })={\overline{D}}_q\Delta \dot{\widehat{q}}+{\overline{L}}_q{\overline{L}}_q^{\text{T}}\Delta \widehat{q}-\overline{b}.& (38)\end{array}$

to determine the estimated model parameter vector $\widehat{\overline{\zeta }}$, i.e., estimated damping ${\widehat{\overline{D}}}_q$ and stiffness ${\widehat{\overline{K}}}_q$.

2.3. Simulation

We validate our method with simulated data of a neuromechanical model of the human arm, which generates transversal movements of a two-link, six-muscle arm through calculation of muscle activities (Franklin et al., 2008). By simulating each movement twice, once with and once without perturbation, we validate the estimation of the unperturbed dynamics. Furthermore, we use the inertial characteristics and the intrinsic impedance to validate the estimation of the involuntary impedance parameters.

We use the neuromechanical model to simulate a static postural task and a dynamic movement task. As the static postural task does not include an underlying movement, we are able to apply the do-not-intervene-voluntarily paradigm. While the do-not-intervene-voluntarily paradigm does not guarantee exclusion of voluntary feedback or constant unperturbed states (Sainburg, 2015), it is widely used and accepted as a plausible approximation of both of these assumptions (Mussa-Ivaldi et al., 1985; Gomi and Kawato, 1997; Osu and Gomi, 1999). Thus, given the absence of an underlying movement, and the fact that the elements of the BIP vector $\overline{\pi }$ are all constant, we are able to use an estimation interval duration T_est that is longer than the delay of voluntary feedback δ_v in the static postural task. As the dynamic movement task does include an underlying movement, it is performed without the do-not-intervene-voluntarily paradigm and we use the limited estimation interval duration T_est = δ_v.

2.3.1. Human Arm Model

The rigid body dynamics of the neuromechanical model are determined by (1) and the correlations between the internal torques τ_int and the muscle tensions m are determined by (2). The muscle tensions

$\begin{array}{ll}m=m_a+m_{\text{imp}},& (39)\end{array}$

where m_a and m_imp represent the muscle tensions due to muscle activation and mechanical impedance, respectively. The muscle tensions due to muscle activation m_a are assumed to be identical to the muscle activations a in (3) and the muscle tensions due to mechanical impedance

$\begin{array}{ll}{m}_{\text{imp}}=D_m{\dot{e}}_m+K_m{e}_m& (40)\end{array}$

$\begin{array}{ll}{D}_m=K_m/12,K_m=K_0+K_{1}a,& (41)\end{array}$

where D_m, K_m, and e_m represent damping, stiffness, and tracking errors with respect to the desired trajectory at the muscle level and K₀ and K₁ contain intrinsic stiffness parameters. The feedforward muscle activations a_FF are calculated a priori based on the inverse kinematics and dynamics of the desired movement. The feedback muscle activations a_FB are modeled by proportional-derivative (PD) control that depends on the tracking errors e_m, the respective derivatives ${\dot{e}}_m$, and a simulated feedback delay δ_s = 60 ms. The signal-dependent noise a_N is calculated with zero mean Brownian motion, which provides movements similar to those observed in previous experiments (Burdet et al., 2001). It is generated by a normally distributed random variable, which is amplified by a parameter α_N = 12.5 and filtered by a fifth order low pass Butterworth filter with a cut-off frequency f_c,N = 2 Hz.

In the implementation of the neuromechanical model, there is no distinction between reflexive and voluntary feedback. Instead, both feedback mechanisms are combined and the simulated feedback delay δ_s is defined to lie between the respective delays δ_r and δ_v. As the simulated feedback delay δ_s is shorter than the duration of our estimation interval T_est, the feedback behavior within the estimation interval is influenced by voluntary contributions. In order to avoid this and enable the validation of the estimated involuntary impedance parameters ${\widehat{\overline{D}}}_q$ and ${\widehat{\overline{K}}}_q$ through the simulated intrinsic impedance components, we use a variation, in which the simulated feedback delay δ_s is defined equal to the delay δ_v of voluntary feedback, i.e., δ_s = δ_v = 115 ms. Although the physiological correctness of the variation is slightly reduced compared to the original implementation, it nonetheless represents a plausible simulation of human behavior. More importantly, it represents a means for validation of our method and for comparison of its performance to those of existing methods. In order to represent the effects of the do-not-intervene-voluntarily paradigm, voluntary feedback is removed in the simulations of the static postural task.

2.3.2. Simulation Design

As the definition of the signal-dependent noise a_N is based on a normally distributed random variable, the static task and the dynamic task are each simulated ten times and the respective results are averaged. In both tasks, the simulated manipulandum inertia and damping are defined according to the admittance inertia M_adm = diag{5, 5} kg and damping D_adm = diag{20, 20} Ns/m of our apparatus.

Static task. In the static task, the arm maintains a total of five different positions that are distributed evenly across the horizontal plane: $x_{\text{P}1}={\left[0,0.35\right]}^{\text{T}}\text{ m}$, $x_{\text{P}2}={\left[-10,0.45\right]}^{\text{T}}\text{ m}$, $x_{\text{P}3}={\left[0,0.55\right]}^{\text{T}}\text{ m}$, $x_{\text{P}4}={\left[10,0.45\right]}^{\text{T}}\text{ m}$, and $x_{\text{P}5}={\left[0,0.45\right]}^{\text{T}}\text{ m}$ (see Figure 5). In order to obtain perturbations that are similar to those of the static estimations in Gomi and Kawato (1997), we define the perturbation profile duration T_p,1 = 160 ms which results in a total perturbation duration T_pert = 400 ms. We change the interaction mode from closed-loop to open-loop which means that applied forces do not contribute to the movement of the manipulandum. As this means that we are essentially applying a position perturbation, we are able to define the amplitude of the perturbation position profile x_pert, which we set to 8 mm. The perturbation angles ϕ_pert are defined by the set Φ_pert = {(π/12)k | k = 1 − 24 \ {6, 12, 18, 24}}, which contains 20 angles. Each perturbation angle ϕ_pert is executed once in every position and the order of the resulting perturbations is randomized. Single execution of each of the 20 perturbation angles ϕ_pert in each of the five positions results in a total of 100 trials. The duration of the estimation interval T_est = 400 ms and the estimated BIP vector $\widehat{\overline{\pi }}$ is obtained by calculating the least squares solution for the complete data set of all 100 trials.

FIGURE 5

Figure 5. Schematic representation of the human arm task space, the static task, and the dynamic task. In the static task, the arm maintains a total of five different positions that are distributed evenly across the horizontal plane. In the dynamic task, it performs point to point movements along the sagittal axis.

Dynamic task. In the dynamic task, the arm performs two-dimensional point to point movements along the sagittal axis from $x_{\text{start}}={\left[0,0.30\right]}^{\text{T}}\text{ m}$ to $x_{\text{goal}}={\left[0,0.55\right]}^{\text{T}}\text{ m}$ (see Figure 5). The duration of the movements T_mov is set to 2 s. The feedforward muscle activations a_FF are calculated with the inverse kinematics and dynamics of a positional data set, which is provided by the authors of Franklin et al. (2008). It contains positional data of 50 two-dimensional point to point arm movements, which possess identical start and goal positions as the movements in our task. The perturbations are designed to generate sufficiently large deviations in as short a time frame as possible, but nonetheless be physically plausible and kinesthetically renderable under hardware limitations. The perturbation amplitude A_pert = 40 N and the perturbation profile duration T_p,1 = 70 ms which corresponds to a total perturbation duration T_pert = 175 ms. The perturbation is initiated when the hand reaches x₂ = 0.4 m which equals a distance of 0.1 m along the axis of the principal movement. The perturbation angles ϕ_pert are defined by Φ_pert, which is identical to the static task. Each perturbation angle ϕ_pert is executed three times and the order of the resulting perturbations is randomized. Three repetitions of each of the 20 perturbation angles ϕ_pert result in a total of 60 trials, of which each consists of an unperturbed and a perturbed movement. The duration of the estimation interval T_est = 115 ms and the estimated damping ${\widehat{\overline{D}}}_q$ and stiffness ${\widehat{\overline{K}}}_q$ are obtained by calculating the least squares solutions for data sets of 20 trials each. Grouping the complete data set of all 60 trials into data sets of 20 trials each, starting from the first 20 trials and moving the respective data set along trial by trial until the last 20 trials, provides a total of 41 individual least squares solutions.

2.4. Experiment

In order to evaluate the performance of the proposed method for real data, we conduct an experiment with 12 participants. The experiment design is identical to the simulation design, apart from a randomized time interval before the onset of the perturbation in the static task and a randomized distribution of unperturbed and perturbed trials in the dynamic task.

2.4.1. Apparatus and Experiment Design

The apparatus is presented in Figure 6. It consists of two linear servo motor driven single rail stages (Copley Controls Thrusttube Module) that are mounted on top of each other in orthogonal orientation. Together, they span a 2-DoF workspace of ±0.20 m and are each equipped with an encoder that yields position data with a precision of 1 μm. Additionally, six motion capture cameras (Qualisys Miqus) are available for tracking of passive markers. A vertical handle and a 6-DoF force-torque sensor (JR3-67M25) are mounted on top of the upper servo motor driven cart to measure the forces in the horizontal plane. A custom-built seat with shoulder belts, a sling attached to the ceiling, and a wrist orthosis are available for limitation of the movements of the participants. Visual feedback is provided on a screen behind the apparatus and implemented with the Psychophysics Toolbox (Brainard, 1997). The position of the cart is visualized by a dot and the workspace safety boundaries are visualized in the form of a boundary box.

FIGURE 6

Figure 6. Reenactment of an individual interacting with the apparatus (during the dynamic task of the experiment). Informed written consent for the publication of this image was obtained from the individual.

Haptic interaction by means of participant force input is enabled by the admittance control scheme defined in (26), where M_adm = diag{5, 5} kg and D_adm = diag{20, 20} Ns/m. The parameters of the admittance control scheme are chosen to guarantee natural interaction with the apparatus and sufficient attenuation of force-torque sensor noise. Precise rendering of the position is ensured by a high gain PD controller, which is implemented in Matlab/Simulink and executed on a Linux system with RT-preempt real-time kernel (Ubuntu 14.04, 3.14.3-rt4-prt). The sample rate is set to f_s = 4 kHz and inputs to the Thrusttube Modules are downsampled to 2 kHz due to hardware limitations. The signals are filtered using a seventh order Savitzky-Golay polynomial filter with a cut-off frequency of f_c,SG = 50 Hz.

In order to avoid predictability of the perturbations and adaptation of the participants, a randomized time interval with T_random ∈ [1, 3] s is included before the onset of the perturbations in the static task. For the same reasons and to obtain both unperturbed and perturbed movements, the ratio of perturbed trials to total amount of trials is set to 33 % in the dynamic task. In order to avoid effects of fatigue, the resulting 180 trials of the dynamic task are performed in three consecutive blocks. Each of these blocks contains a total of 60 trials that consist of 20 perturbed and 40 unperturbed randomly distributed trials. In both tasks, the visual feedback of the position of the cart is deactivated during the perturbations.

2.4.2. Participants and Experiment Procedure

A total of 12 participants (9 males, 3 females) with between-subject mean (SD) age of 26.92 (3.40) years, height of 174.83 (7.38) cm, and weight of 69.43 (9.82) kg volunteered to participate in this experiment. All 12 participants had normal or corrected-to-normal vision. One of the participants was left handed and performed the experiment with the non-dominant hand. The remaining 11 participants were right handed and participated with the dominant hand. Informed written consent was obtained from all participants before the experiment, which was conducted according to the principles in the Declaration of Helsinki. The research ethics were obtained from the ethics committee at the Technical University of Munich.

The participants were seated in front of the apparatus and their upper body was restrained by two shoulder belts to fix the position of the shoulder. Their upper arm was supported by a sling attached to the ceiling to constrain all motions to the horizontal plane and reduce effects of fatigue. In order to avoid wrist motions, the participants were wearing a wrist orthosis. Two passive motion capture markers were placed on shoulder and elbow to measure the lengths of the upper arm and the forearm.

In the beginning of the static task, the cart automatically moved to the first position at the bottom of the workspace. After the application of all 20 perturbations, it automatically moved to the next position. This procedure was repeated until all 100 perturbations were completed. Analogous to the simulation, each perturbation angle ϕ_pert was executed once in each of the positions and the order was randomized. The participants were informed about the procedure, including the perturbations and the randomized time interval before the onset of the perturbation. They were told that their objective was to naturally grasp the handle on top of the cart. Additionally, they were instructed to not voluntarily react to the perturbations in any way and not prepare for them in any kind of preemptive manner, e.g., by co-contraction.

In the beginning of the dynamic task, the cart automatically moved to the start position at the bottom of the workspace. As soon as it reached the start position, the admittance turned on. As soon as the participants reached the goal position, represented by a dot at the top of the workspace, the admittance turned off and the cart automatically moved back down to the start position. This procedure was repeated until all 60 trials with 20 randomly distributed perturbed trials were completed. In total, three of these blocks of trials were completed for a total of 180 trials. Analogous to the simulation, each perturbation angle ϕ_pert was executed three times and the order was randomized. The participants were informed about the procedure, including the perturbations and their random distribution. They were told that their objective was to naturally grasp the handle on top of the cart and move it from the start position to the goal position. They were told that the duration of the movement should be ~ 2 s. In order to help them adjust their movement velocity, a beep sound occurred after 2 s. The participants were informed that they were allowed to voluntarily react to the perturbations. Additionally, they were instructed to not prepare for the perturbations in any kind of preemptive manner, e.g., by co-contraction.

3. Results

In this section, we first present the results of the simulation, which consist of the validation of the feedback jerk isolation and the involuntary impedance estimation. Subsequently, we present the results of the evaluation with experimental data. The results of the simulation are either presented in within-, between-, or cross-simulation mean (SD). The cross-simulation mean (SD) is obtained by calculating the mean of the within-simulation means and the mean of the within-simulation SDs. Analogous correlations apply for the results of the experiment, which are presented in within-, between-, or cross-subject mean (SD). As the majority of the results are cross-simulation/cross-subject mean (SD) results, for simplicity, from this point on, the respective results are referred to simply as mean (SD) results.

3.1. Validation Feedback Jerk Isolation

In this section, we validate the feedback jerk isolation by analyzing its performance for different movement velocities and variations of neural noise. Furthermore, we compare the results to those obtained by application of the methods in Gomi and Kawato (1997) and Erden and Billard (2015). In order to ensure equal conditions and enable performance comparisons without effects of voluntary feedback, all methodologies are applied to the same sets of simulated data and use the same duration of the estimation interval T_est = 115 ms. The estimation of the unperturbed dynamics $\{q^{*},{\dot{q}}^{*},{\ddot{q}}^{*},{\tau }_{\text{ext}}^{*}\}$ is validated by analysis of the estimated variational dynamics $\{\Delta \widehat{x},\Delta \dot{\widehat{x}},\Delta \ddot{\widehat{x}},\Delta {\widehat{u}}_{\text{ext}}\}$. The estimation accuracy is assessed using the normalized root mean square errors (NRMSEs) in the estimation interval [0, T_est]:

$\begin{array}{ll}\text{NRMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\left(\frac{{\chi }_i-{\widehat{\chi }}_i}{{\nu }_i}\right)}^2},& (42)\end{array}$

where ν is the normalizing value, χ is the real value, $\widehat{\chi }$ is the estimated value, and n = 2 is the dimensionality. For the estimated variational dynamics $\{\Delta \widehat{x},\Delta \dot{\widehat{x}},\Delta \ddot{\widehat{x}},\Delta {\widehat{u}}_{\text{ext}}\}$, the normalizing value ν is given by the maximum real value in the estimation interval [0, T_est]. In order to analyze the performance for different movement velocities, the duration of the movements T_mov is changed (1, 3 s). For analysis of the performance for different variations of neural noise, the parameterization of the zero mean Brownian motion is changed in terms of the cut-off frequency f_c,N (1, 3 Hz) and the amplitude α_N (5, 20).

3.1.1. Filter Configuration

Figure 7A shows the mean results of the ESDs ψ_UP and ψ_P. In the unperturbed movements, the energy of the principal movement along the x₂ axis is distinguishable by a peak in the respective ESD ψ_UP,2. A secondary, much lower peak represents effects of neural noise. As the point to point movements do not require movements along the x₁ axis, the respective ESD ψ_UP,1 is significantly lower and represents effects of neural noise. The energy of both axes of the unperturbed movements decreases to marginally low values for high frequencies. The opposite applies for the energy of the perturbed movements, which increases to significantly higher values for high frequencies. The respective ESDs ψ_P are much higher than those of the unperturbed movements and have multiple peaks in the high frequency range.

FIGURE 7

Figure 7. Filter configuration based on simulated data. (A) Mean results of the ESDs ψ_UP and ψ_P of the jerks $\stackrel{⃛}{x}$ of the unperturbed and perturbed movements. (B) Calculation of the cut-off frequencies f_c,HP, defined as the frequencies f at which the ESDs ψ_FB of the jerks $\stackrel{⃛}{x}$ of the evoked feedback behavior become larger than the ESDs ψ_UP of the jerks $\stackrel{⃛}{x}$ of the unperturbed movements. For clarity, the calculation is exemplarily illustrated for a single simulation and both graphs are limited to the low frequency ranges. The respective cut-off frequencies f_c,HP are indicated by vertical lines.

Figure 7B shows the ESDs ψ_UP, ψ_P, and ψ_FB of a single simulation in two axis-specific graphs to demonstrate the calculation of the cut-off frequencies f_c,HP. The cut-off frequency f_c,HP,2 of the x₂ axis is sufficiently high to lie above the energy of the principal movement. Concurrently, it is sufficiently low to lie below the energy of the feedback behavior. Due to the absence of movements along the x₁ axis, the cut-off frequency f_c,HP,1 is slightly lower. The between-simulation mean (SD) results of the cut-off frequencies f_c,HP of all ten simulations are f_c,HP = [1.89 (0.08), 2.24 (0.04)] Hz.

3.1.2. Results

The mean results of the NRMSEs in Figure 8 show that the NRMSEs all increase with t_p. The gradual increase of the NRMSEs of the estimated variational accelerations $\Delta \ddot{\widehat{x}}$ results from the integration of the estimated variational jerks $\Delta \stackrel{⃛}{\widehat{x}}$, which include estimation errors that originate from the calculation of the high pass filtered jerks ${\stackrel{⃛}{x}}_{\text{HP}}$, i.e., the separation of the unperturbed and the variational behavior. The increases of the NRMSEs of the estimated variational velocities $\Delta \dot{\widehat{x}}$ and positions $\Delta \widehat{x}$ result from the consecutive integrations of the estimated variational accelerations $\Delta \ddot{\widehat{x}}$. The slightly larger NRMSEs of the estimated variational external forces $\Delta {\widehat{u}}_{\text{ext}}$ are caused by the estimation errors of the estimated unperturbed velocities ${\dot{\widehat{x}}}^{*}$ and accelerations ${\ddot{\widehat{x}}}^{*}$. Despite the increase of the NRMSEs with t_p, the maximum values of the estimated variational kinematics $\{\Delta \widehat{x},\Delta \dot{\widehat{x}},\Delta \ddot{\widehat{x}}\}$ are all below 6 % and the NRMSE of the estimated variational external forces $\Delta {\widehat{u}}_{\text{ext}}$ is below 8 %.

FIGURE 8

Figure 8. Validation of feedback jerk isolation with simulated data. Mean results of the NRMSEs of the estimated variational dynamics $\{\Delta \widehat{x},\Delta \dot{\widehat{x}},\Delta \ddot{\widehat{x}},\Delta {\widehat{u}}_{\text{ext}}\}$ in the estimation interval [0, T_est].

3.1.3. Comparison

In order to compare the results to those obtained with the methods in Gomi and Kawato (1997) and Erden and Billard (2015), we average the NRMSEs over the complete estimation interval [0, T_est]. For simplicity, from this point on, we refer to the resulting averaged NRMSEs simply as NRMSEs. The mean (SD) results of the NRMSEs for the original simulation are listed in Table 1A and those for different durations of the movements T_mov, cut-off frequencies f_c,N, and amplitudes α_N are listed in Tables 1B–D, respectively. The abbreviations FJI, GOM, and ERD represent our feedback jerk isolation and the methods in Gomi and Kawato (1997) and Erden and Billard (2015), respectively. In the discussion of the NRMSEs of the alternative simulations, the changes in configurations and results are evaluated relative to the original simulation. In order to facilitate the comparisons of the NRMSEs of the original and alternative simulations, Tables 1B–D also include the NRMSEs of the original simulation, and the contents of all three tables are additionally illustrated with intermediate values in Figure 9.

TABLE 1

Table 1. Validation of feedback jerk isolation with simulated data.

FIGURE 9

Figure 9. Validation of feedback jerk isolation with simulated data. Mean results of the NRMSEs of the estimated variational dynamics $\{\Delta \widehat{x},\Delta \dot{\widehat{x}},\Delta \ddot{\widehat{x}},\Delta {\widehat{u}}_{\text{ext}}\}$, averaged for the estimation interval [0, T_est], over different movement durations T_mov, cut-off frequencies f_c,N, and amplitudes α_N.

The NRMSEs in Table 1A show that FJI achieves high estimation accuracy and outperforms both GOM and ERD. The magnitudes and differences of the individual NRMSEs of FJI are in accordance with those of the NRMSEs presented in Figure 8. The performance difference of FJI and GOM increases from a small difference in $\Delta \ddot{\widehat{x}}$ to a large difference in $\Delta \widehat{x}$ and the performance difference in $\Delta {\widehat{u}}_{\text{ext}}$ is marginally larger than that of $\Delta \ddot{\widehat{x}}$. These correlations are plausible, due to the dependency on the averaged unbiased dynamics of all perturbed movements. While forces and accelerations are of similar magnitude, velocities and especially positions vary between movements. As a consequence, the accuracy of the averaged unbiased dynamics decreases significantly from accelerations down to positions. Similar correlations apply to the performance difference of FJI and ERD, but with significantly larger differences, especially for $\Delta \dot{\widehat{x}}$ and $\Delta \widehat{x}$. The magnitudes of the NRMSEs of ERD are plausible, due to the dependency on the constant values of the offsets at perturbation onset. While these constant values represent accurate approximations of the unperturbed states of the quasi-static movements of the manual welding task in Erden and Billard (2015), they are incapable of accurately representing the kinematics of the unperturbed states of dynamic point to point movements. For this reason, we exclude this method from the remainder of the validation of the feedback jerk isolation and the involuntary impedance estimation.

The NRMSEs in Table 1B show that the performance of FJI is decreased for a higher movement velocity (resulting from a shorter duration of the movement) and marginally increased for a lower movement velocity (resulting from a longer duration of the movement). A higher movement velocity results in an increased cut-off frequency of the high pass filter which causes an increased information loss in the isolation of the feedback jerk. As a result, the estimation accuracy is decreased. The marginal increase in performance for a lower movement velocity suggests that the corresponding NRMSEs are close to the smallest possible errors, which are caused by unavoidable frequency overlap of unperturbed movements and feedback behavior, neural noise, and imperfect properties of a non-ideal high pass filter. The performance of GOM is slightly increased for a higher movement velocity and slightly decreased for a lower movement velocity. As the accumulated influence of the neural noise at the onset of the perturbation is negatively correlated with the movement velocity, the deviations from the averaged unbiased dynamics are larger for slow movements than they are for fast movements. As the performance of GOM is directly influenced by these deviations, the estimation accuracy is positively correlated with the different movement velocities. For T_mov = 1 s, GOM outperforms FJI in all NRMSEs except $\Delta \widehat{x}$, for which the performance difference is decreased compared to Table 1A. For T_mov = 3 s, FJI outperforms GOM with slightly increased performance differences compared to Table 1A.

The NRMSEs in Table 1C show that the performance of FJI is slightly decreased for lower and higher cut-off frequencies of the noise, with the latter resulting in slightly larger performance differences than the former. It seems that the lower cut-off frequency of the noise results in such a decreased cut-off frequency of the high pass filter, that it is too close to the frequency content of the unperturbed behavior. For the higher cut-off frequency of the noise, the increased cut-off frequency of the high pass filter results in an increased information loss in the isolation of the feedback jerk. Despite these slight decreases in performance, FJI still achieves high estimation accuracy for both alternative simulations. The performance of GOM is increased/decreased for lower/higher cut-off frequencies of the noise due to smaller/larger deviations from the averaged unbiased dynamics. For f_c,N = 1 Hz, GOM outperforms FJI in all NRMSEs except $\Delta \widehat{x}$, for which the performance difference is decreased compared to Table 1A. For f_c,N = 3 Hz, FJI outperforms GOM with slightly increased performance differences compared to Table 1A.

The NRMSEs in Table 1D illustrate that a lower amplitude of the noise has an almost identical effect on the performance as a lower cut-off frequency of the noise in Table 1C. In contrast, a higher amplitude of the noise has very different effects. For FJI, a higher amplitude of the noise only leads to a marginal decrease of the performance, much smaller than that in Table 1C. For GOM, it leads to a similar decrease of the performance as in Table 1C, except for $\Delta \widehat{x}$, for which the decrease is much larger. For α_N = 5, the performance difference of FJI and GOM is almost identical to that of f_c,N = 1 Hz in Table 1C. For α_N = 20, FJI outperforms GOM with increased performance differences compared to Table 1A, especially for $\Delta \widehat{x}$. The different effects emphasize three issues: (1) For GOM, the estimation accuracy strongly depends on the similarity to the averaged unbiased dynamics. (2) For FJI, it depends more strongly on the frequency of the noise than it does on the amplitude, due to the influence on the cut-off frequency of the high pass filter. (3) The performance difference in GOM due to a higher amplitude of the noise is larger than the performance difference in FJI due to a higher frequency of the noise, especially for $\Delta \widehat{x}$.

3.2. Validation Involuntary Impedance Estimation

In this section, we validate and compare the performance of the involuntary impedance estimation. The estimation accuracy is assessed using the normalized absolute errors (NAEs) of the estimated values. The normalizing values are either given by the real values of the elements of the BIP vector $\overline{\pi }$ or the maximum real values of the elements of the damping ${\overline{D}}_q$ and stiffness ${\overline{K}}_q$, which are obtained by averaging the respective elements for the estimation interval [0, T_est]. For comparability of the results, we transform the inertial parameters in Gomi and Kawato (1997) to the elements of the BIP vector $\overline{\pi }$.

In order to obtain additional means of comparing the overall performance of the methods, we calculate two additional performance criteria that are based on the Akaike information criterion (AIC) and the Bayesian information criterion (BIC). When comparing least squares fitted models, the AIC and the BIC are defined by the residual sum of squares (RSS) of the real values and the estimated model outputs:

$\begin{array}{ll}\text{AIC}=2p+kln(\text{RSS}),& (43)\end{array}$

$\begin{array}{ll}\text{BIC}=ln(k)p+kln(\text{RSS}),& (44)\end{array}$

$\begin{array}{ll}\text{RSS}=\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\displaystyle ({\displaystyle \sum _{j=1}^k}{\left(y_{i,j}-{\widehat{y}}_{i,j}\right)}^2)},& (45)\end{array}$

where y is the real value, ŷ is the estimated model output, p is the number of parameters, k is the number of data points, and n = 2 is the dimensionality (Burnham and Anderson, 2002). As we aim to compare the overall performance of the methods, we define the estimated model output ŷ to be the simulation output that is obtained by inserting the estimated BIP vector $\widehat{\overline{\pi }}$, damping ${\widehat{\overline{D}}}_q$, and stiffness ${\widehat{\overline{K}}}_q$ into the neuromechanical model of the human arm (within the estimation interval [0, T_est]), and the real value y to be the corresponding simulation output of the original simulation.

3.2.1. Results

Table 2A contains the mean (SD) results of the NAEs of the estimated BIP vector $\widehat{\overline{\pi }}$, damping ${\widehat{\overline{D}}}_q$, and stiffness ${\widehat{\overline{K}}}_q$. The NAEs of the elements of the estimated BIP vector $\widehat{\overline{\pi }}$ do not possess SDs, because they are estimated with the complete data set of the static task. The NAEs of all three elements are below 2% which indicates high estimation accuracy of the combination of the dynamic regressor representation with the data of the static task. The NAEs of the elements of the estimated damping ${\widehat{\overline{D}}}_q$ and stiffness ${\widehat{\overline{K}}}_q$ are all approximately equal to or below 10% and demonstrate high estimation accuracy of the non-linear least squares estimation with the data of the dynamic task. The NAEs of the elements ${\widehat{\overline{D}}}_{q,11}$ and ${\widehat{\overline{K}}}_{q,11}$ are slightly increased compared to the remaining elements of the respective matrices. This slight increase is plausible, as these elements represent the contributions of the single-joint shoulder muscles, which are less involved due to less movement of the shoulder joint.

TABLE 2

Table 2. Validation of involuntary impedance estimation with simulated data.

3.2.2. Comparison

The NAEs of the elements of the estimated BIP vector $\widehat{\overline{\pi }}$ of GOM are almost identical to those of FJI. The mean NAEs of the elements of the estimated damping ${\widehat{\overline{D}}}_q$ of GOM are marginally smaller than those of FJI and the opposite applies to the corresponding SDs which essentially makes these NAEs almost identical as well. In contrast, a considerable difference in estimation accuracy is found in the elements of the estimated stiffness ${\widehat{\overline{K}}}_q$, for which both the mean NAEs and the SDs of GOM are larger than those of FJI, with the mean NAEs being approximately twice as large for GOM as they are for FJI. These results are plausible, as the difference in estimation accuracy of the variational dynamics $\left\{\Delta x,\Delta \dot{x},\Delta \ddot{x},\Delta u_{\text{ext}}\right\}$ is largest for the variational positions Δx. According to the involuntary impedance model, the difference in estimation accuracy of the variational angles Δq directly influences the estimated stiffness ${\widehat{\overline{K}}}_q$.

The mean (SD) results of the RSSs, AICs, and BICs of the variational external torques Δτ_ext are presented in Table 2B. FJI outperforms GOM in all three performance criteria. The difference in RSS is especially relevant, as it demonstrates, that the differences in AIC and BIC do not arise solely due to differences in the number of parameters p. Similar to the mean NAEs of the elements of the estimated stiffness ${\widehat{\overline{K}}}_q$ in Table 2A, the mean RSSs are approximately twice as large for GOM as they are for FJI. The differences in RSS, AIC, and BIC demonstrate that (1) the differences in the estimated BIP vector $\widehat{\overline{\pi }}$, damping ${\widehat{\overline{D}}}_q$, and stiffness ${\widehat{\overline{K}}}_q$, represented by the respective NAEs in Table 2A, have a substantial effect on the replicability of the real simulation output and that (2) the involuntary impedance estimation results of FJI yield a more accurate replication of the real simulation output than those of GOM.

3.3. Evaluation

In this section, we evaluate the performance of our method for real data obtained in an experiment with 12 human participants. For simplicity, from this point on, we omit the term estimated when referring to the results of the involuntary impedance estimation, i.e., BIP vector $\widehat{\overline{\pi }}$, damping ${\widehat{\overline{D}}}_q$, and stiffness ${\widehat{\overline{K}}}_q$. In the experiment, every trial is either an unperturbed or a perturbed trial. Consequently, we can not evaluate the estimation accuracy of the unperturbed states. Therefore, we only evaluate the calculation of the cut-off frequencies f_c,HP of the high pass filter.

3.3.1. Feedback Jerk Isolation

Figure 10 shows the mean results of the ESDs ψ_UP and ψ_P, which look very similar to those of the simulation presented in Figure 7A. In the unperturbed movements of the experiment, the energy of the principal movement along the x₂ axis is more widespread than in the simulation. However, the respective ESD ψ_UP,2 nonetheless possesses a peak at a similar frequency. The ESD ψ_UP,1 is significantly lower, due to the absence of movements along the x₁ axis. Due to hardware noise, the energy of both axes of the unperturbed movements does not decrease to values as low as those of the simulation for high frequencies. However, the respective values are nonetheless significantly smaller than those of the perturbed movements, for which the energy increases to significantly higher values for high frequencies. The respective ESDs ψ_P are much higher than those of the unperturbed movements and have multiple peaks in the high frequency range. While the overall energy of the unperturbed movements is slightly decreased, the overall energy of the perturbed movements is increased compared to the simulation. This difference is largely caused by high frequency oscillations in the jerks $\stackrel{⃛}{x}$ of the perturbed movements which result from oscillations due to the perturbations paired with sensor noise and high gain PD control. Nonetheless, the between-subject mean (SD) results of the cut-off frequencies f_c,HP = [1.45 (0.30), 2.09 (0.24)] Hz are very similar to those of the simulation, with both values being slightly lower due to the slight decrease in overall energy of the unperturbed movements. Nonetheless, the cut-off frequency f_c,HP,2 is still sufficiently high to lie above the energy of the principal movement.

FIGURE 10

Figure 10. Filter configuration based on experimental data. Mean results of the ESDs ψ_UP and ψ_P of the jerks $\stackrel{⃛}{x}$ of the unperturbed and perturbed movements.

3.3.2. Involuntary Impedance Estimation

Table 3 lists the within- and cross-subject mean (SD) results of the involuntary impedance estimation. Analogous to the simulation, the elements of the BIP vector $\widehat{\overline{\pi }}$ do not possess SDs, because they are estimated with the complete data set of the static task. The elements of the mean BIP vector $\widehat{\overline{\pi }}$ are similar to the simulated values in the simulation π_sim = [0.1945, 0.0737, 0.0838]. The first element ${\widehat{\overline{\pi }}}_1$ of the BIP vector is decreased compared to the respective value in the simulation. As this element largely depends on the inertial characteristics of the upper arm, the difference could be caused by the comparatively little movement of the shoulder joint which results in less involvement of the upper arm. As the inertial parameters in Gomi and Kawato (1997) are defined by the linearized rigid body dynamics, the respective estimation results can not be used for comparison. Thus, we use the elements of the mean BIP vector $\widehat{\overline{\pi }}$ to calculate the mean inertia ${\widehat{\overline{M}}}_q$, which we transform to Cartesian space to obtain the mean Cartesian endpoint inertia ${\widehat{\overline{M}}}_x$. The elements of this mean Cartesian endpoint inertia ${\widehat{\overline{M}}}_x$ are similar to those reported in existing studies (Tsuji et al., 1995; Chang et al., 2013; Dyck and Tavakoli, 2013). Consequently, the same applies to the Cartesian endpoint ellipse, which is illustrated in Figure 11.

TABLE 3

Table 3. Evaluation of involuntary impedance estimation with experimental data.

FIGURE 11

Figure 11. Evaluation of involuntary impedance estimation with experimental data. Cartesian endpoint ellipses of the mean inertia ${\widehat{\overline{M}}}_x$, damping ${\widehat{\overline{D}}}_x$, and stiffness ${\widehat{\overline{K}}}_x$. The respective Cartesian space matrices are calculated by transformation of the BIP vector $\widehat{\overline{\pi }}$, damping ${\widehat{\overline{D}}}_q$, and stiffness ${\widehat{\overline{K}}}_q$, respectively.

The elements of the mean damping ${\widehat{\overline{D}}}_q$ are similar to the averaged simulated values D_q,sim = [2.42, 1.20, 1.20, 1.42], with the exception of the element ${\widehat{\overline{D}}}_{q,11}$, which is comparatively large. As the element ${\widehat{\overline{D}}}_{q,11}$ represents the contributions of the single-joint shoulder muscles, which are less involved due to less movement of the shoulder joint, it is more difficult to estimate. This correlation is also visible in the increased estimation error in the validation with simulated data (see Table 2A). As the remaining elements are decreased compared to the respective values in the simulation, it is possible that some of the contributions of these elements are allocated to the element ${\widehat{\overline{D}}}_{q,11}$. As the damping results of Gomi and Kawato (1997) are not reported, we can not use them for comparison. However, the elements of the mean damping ${\widehat{\overline{D}}}_q$ are of similar order of magnitude as those of static tasks (Tsuji et al., 1995; Lakatos et al., 2011). The overall increase in magnitude compared to the results of these static tasks is to be expected, as similar correlations are observed for estimates of stiffness (Gomi and Kawato, 1997).

The elements of the mean stiffness ${\widehat{\overline{K}}}_q$ are similar to the averaged simulated values K_q,sim = [29.05, 14.37, 14.37, 17.08], with the elements ${\widehat{\overline{K}}}_{q,12/21}$ being slightly decreased in comparison. Consequently, the elements of the mean stiffness ${\widehat{\overline{K}}}_q$ are also similar to those reported for comparable dynamic tasks (Burdet et al., 2000; Darainy et al., 2007; Wong et al., 2009b), including those reported for the dynamic task in Gomi and Kawato (1997). The difference in size of the SDs of the stiffness ${\widehat{\overline{K}}}_q$ and the damping ${\widehat{\overline{D}}}_q$ fits to the difference in estimation errors in the validation with simulated data (see Table 2A). It could however also indicate that variations in damping during the course of the experiment are generally lower than those in stiffness or that some of the variations in damping are incorrectly interpreted as variations in stiffness by the non-linear least squares estimation.

Figure 11 shows the Cartesian endpoint ellipses of the mean inertia ${\widehat{\overline{M}}}_x$, damping ${\widehat{\overline{D}}}_x$, and stiffness ${\widehat{\overline{K}}}_x$. The shapes and orientations of the ellipses are similar to those in existing studies (Tsuji et al., 1995; Gomi and Osu, 1998; Darainy et al., 2007). Similar to the results for the movements along the sagittal axis in Gomi and Kawato (1997), the major axis of the Cartesian endpoint ellipse of the mean stiffness ${\widehat{\overline{K}}}_x$ is oriented slightly more parallel to the x₂ axis, i.e., the axis of the principal movement. Figure 12 shows the between-subject mean (SD) results of the BIP vector $\widehat{\overline{\pi }}$, damping ${\widehat{\overline{D}}}_q$, and stiffness ${\widehat{\overline{K}}}_q$ for different durations of the estimation interval T_est. In the static task, the elements of the BIP vector $\widehat{\overline{\pi }}$ converge to constant values and do not change for durations T_est > 400 ms. Similar behavior is observable in the dynamic task for the elements of the damping ${\widehat{\overline{D}}}_q$, with a slight decrease for durations T_est > 115 ms. In comparison, the elements of the stiffness ${\widehat{\overline{K}}}_q$ converge at a slower rate and require a longer estimation interval to reach plausible values. Furthermore, the decrease for durations T_est > 115 ms is larger than that of the elements of the damping ${\widehat{\overline{D}}}_q$. This decrease could be caused by the return to the unperturbed states. The longer the estimation interval, the larger the percentage of the variational data with small deviations from the unperturbed states. The larger this percentage becomes, the more influence the respective values of the elements of the stiffness ${\widehat{\overline{K}}}_q$ have on the solution of the non-linear least squares estimation. The decrease could however also be caused by the activation of voluntary feedback contributions.

FIGURE 12

Figure 12. Evaluation of involuntary impedance estimation with experimental data. Between-subject mean (SD) results of the BIP vector $\widehat{\overline{\pi }}$, damping ${\widehat{\overline{D}}}_q$ and stiffness ${\widehat{\overline{K}}}_q$ for different durations of the estimation interval T_est. In order to avoid overlap, the error bars represent ±0.5between-subject SD. The solid vertical lines indicate the durations of the estimation intervals of the static task (T_est = 400ms) and the dynamic task (T_est = 115ms). The estimation intervals begin at the onset of the perturbations. The dashed vertical lines indicate the offset of the perturbations. For the static task, the duration of the estimation interval T_est and the offset of the perturbation, i.e., the solid and the dashed lines, coincide.

4. Discussion

The mean results of the involuntary impedance estimation for different durations T_est in Figure 12 show that the elements of the damping ${\widehat{\overline{D}}}_q$ and the stiffness ${\widehat{\overline{K}}}_q$ reach plausible values for the estimation interval with duration T_est = δ_v = 115 ms. While the elements of the damping ${\widehat{\overline{D}}}_q$ already reach plausible values at T_est ≈ 75 ms, the amount of information necessary for estimation of plausible values of the elements of the stiffness ${\widehat{\overline{K}}}_q$ is not reached until T_est≈110 ms. The mean results in Figure 12 show that the elements of the BIP vector $\widehat{\overline{\pi }}$ converge to plausible values for the estimation interval with duration T_est = 400 ms. These results, in combination with the mean (SD) results in Table 3, which are similar to those reported for comparable dynamic tasks, successfully demonstrate the applicability of our method to real data.

The mean (SD) results of the NRMSEs in Table 1A show that our feedback jerk isolation achieves higher estimation accuracy for the variational dynamics $\{\Delta \widehat{x},\Delta \dot{\widehat{x}},\Delta \ddot{\widehat{x}},\Delta {\widehat{u}}_{\text{ext}}\}$ than the methods reported in Gomi and Kawato (1997) and Erden and Billard (2015). The difference in estimation accuracy is especially large for the variational positions $\Delta \widehat{x}$. As a consequence, the estimation performance of the involuntary impedance estimation is increased for the elements of the stiffness ${\widehat{\overline{K}}}_q$, which is shown by the mean (SD) results of the NAEs in Table 2A. The mean (SD) results of the NRMSEs for different simulation configurations in Tables 1B–D show that the estimation accuracy of the feedback jerk isolation is decreased for a higher movement velocity. This is plausible, as a higher movement velocity results in an increased cut-off frequency of the high pass filter which causes an increased information loss in the isolation of the feedback jerk. Due to similar reasons, a higher cut-off frequency of the neural noise also leads to a decrease of the estimation accuracy. For all of the remaining simulation configurations, changes in the estimation accuracy are marginal. While the estimation accuracy of the method in Gomi and Kawato (1997) is similarly decreased for a higher cut-off frequency of the neural noise, it is contrastly increased for a higher movement velocity. As this method depends on the similarity of the movement to the averaged unbiased dynamics, an increase of the amplitude of the neural noise leads to a significant decrease of estimation accuracy, especially in the variational positions $\Delta \widehat{x}$. In summary, while the feedback jerk isolation is outperformed by the method in Gomi and Kawato (1997) for movements with high velocity and low movement variance due to neural noise, it provides superior estimation performance for movements with moderate to low velocity and moderate to high movement variability due to neural noise, as it is much less affected by a decrease in the similarity of the movements.

According to the main ISO safety standard for robots (DIN EN ISO 10218-1:2011), the maximum robot end-effector velocity during collaboration with a human must not exceed 250 mm/s (Colgate et al., 2008). Some studies on safe physical human-robot collaboration use more conservative values for the maximum robot end-effector velocity, e.g., 150 mm/s in Neranon (2020) and 100 mm/s in Weitschat et al. (2018). In the simulated data used to obtain the NRMSEs in Table 1B, the mean (SD) peak velocities ẋ_peak that correspond to the different durations of the movements T_mov are ẋ_peak(T_mov = 1 s) = 747.8 (6.3) mm/s, ẋ_peak(T_mov = 2 s) = 377.6 (8.9) mm/s, and ẋ_peak(T_mov = 3 s) = 254.0 (9.8) mm/s. Thus, the peak velocities ẋ_peak of those movements, for which the method in Gomi and Kawato (1997) outperforms the feedback jerk isolation, are far above the current constraints for safe pHRI. For those movements with peak velocities ẋ_peak that are similar to or moderately increased compared to the current constraints, the feedback jerk isolation provides superior estimation performance. As the neural noise parameters in the simulation are calculated to provide movement variability similar to that observed in repetitive, straight reaching movements (Burdet et al., 2001) and both neural and kinematic variability are shown to be interrelated individual characteristics (Haar et al., 2017), it is unlikely that realistic pHRI scenarios possess a lower amount of movement variability. Thus, based on these circumstances, we conclude that the feedback jerk isolation is well-suited for our envisaged application in realistic pHRI scenarios.

We approximate the feedback behavior by a linearized model. This approximation is analogously applied in comparable studies (Dolan et al., 1993; Burdet et al., 2000; Darainy et al., 2007) that include deviations of similar size or larger than the ones in our work. In our static task, the position perturbations with amplitudes of 8 mm result in mean (SD) maximum external force deviations of 4.75 (1.80) N. In our dynamic task, the force perturbations result in mean (SD) maximum position deviations of 3.53 (0.34) mm and external force deviations of 9.48 (3.05) N. The mean (SD) results of the NAEs in Table 2A shows that the estimation errors of the elements of the BIP vector $\widehat{\overline{\pi }}$ are almost identical for the dynamic regressor representation and the linearized rigid body dynamics of Gomi and Kawato (1997). This indicates that the variational dynamics $\{\Delta q,\Delta \dot{q},\Delta \ddot{q},\Delta {\tau }_{\text{ext}}\}$ of the static task are sufficiently small to allow for linearization of the rigid body dynamics without loss of estimation accuracy and supports the assumption that the feedback behavior evoked by the force perturbations can be approximated by a linearized model.

Another advantage of perturbations with comparatively small amplitudes is that they are less likely to lead to instability of the movement than those with larger amplitudes. The region of stability of the movement is influenced by a number of factors, including feedforward and voluntary feedback behavior, which are both highly task-specific. As stiffness is positively correlated with internal torques, movements that require larger feedforward torques, e.g., due to increased movement velocity or interaction forces, coincide with larger stiffness (Tee et al., 2004). Thus, for these movements, perturbations with a certain amplitude result in smaller deviations. Assuming that the voluntary feedback behavior can be modeled by stochastic optimal control, the effects of the perturbations after the delay of voluntary feedback depend on the incorporated cost functions, e.g., defined by the tracking error or the smoothness of the movement (Todorov and Jordan, 2002). Furthermore, the region of stability is also heavily influenced by the level of muscle cocontraction, which is known to increase in response to a number of factors, e.g., incomplete or incorrect internal models (Tomi et al., 2008) and increased accuracy requirements (Lametti et al., 2007). In Burdet et al. (2001), participants increase their level of muscle cocontraction in order to successfully compensate effects of unpredictable environmental instabilities. Due to all of these factors, assessment of the region of stability must take into consideration task-specific influences and constraints. In future work, we aim to apply our method to more complex movements and complement corresponding assessments of the region of stability, e.g., in the form of numerical assessments based on Monte Carlo methods.

Multiple studies present impedance estimation methods based on Kalman filters (Deng et al., 2006; Asao et al., 2013) or extended Kalman filters (Roveda et al., 2013). The advantage of these filters is that, given continuous observations of the relevant variables, they allow for estimation of time-varying values of damping and stiffness. However, their application requires a model of the complete system dynamics, including feedforward and voluntary feedback behavior (if the respective method is used for a duration that is longer than the delay of voluntary feedback). Some studies avoid these limitations by assuming that the combination of feedforward and feedback behavior can be modeled by the sum of damping and stiffness (Asao et al., 2013) or just stiffness (Roveda et al., 2013). In Deng et al. (2006), the authors assume that effects of feedforward and voluntary feedback behavior are neglectable due to application of pseudo-random perturbations in combination with a band pass filter. However, the plausibility of these assumptions is not validated, as the method is only applied to a simulation that models the combination of feedforward and feedback behavior by the sum of damping and stiffness. Thus, existing impedance estimation methods based on Kalman filters are extremely limited in possible application scenarios and the application to realistic pHRI would require significantly more complex models. While our method is not able to provide time-varying values of damping and stiffness, it is able to provide accurate estimates within a limited interval without the need for modeling feedforward and voluntary feedback behavior.

5. Conclusion

In this work, we present a method for the estimation of the involuntary impedance components of the human arm during multi-joint movements by application of force perturbations. These perturbations are designed such that the evoked feedback jerk frequency content can be isolated with a high pass filter. We limit the duration of the estimation interval to 115 ms to guarantee exclusion of voluntary feedback. We estimate the inertial parameters in a static postural task and subsequently insert them to estimate the damping and stiffness in a dynamic movement task.The evaluation of the experimental data shows that our method is able to provide plausible involuntary impedance estimates within the limited estimation interval that guarantees exclusion of voluntary feedback. Furthermore, the validation with simulated data shows that it provides superior estimation performance compared to the results obtained by application of existing methodologies within identical conditions. As the difference in estimation accuracy is especially large for the variational positions, the estimation performance is increased for the elements of the stiffness. The analysis of different movement velocities and variations of neural noise shows that the feedback jerk isolation is able to provide superior estimation performance for movements with moderate to low velocity and is much less affected by an increase in movement variability. We conclude that our method allows for involuntary impedance estimation in experiments that emulate realistic pHRI, without the need to include the do-not-intervene-voluntarily paradigm or comparable constraints within the respective dynamic movement tasks. This enables the acquisition of valuable information regarding involuntary impedance modulation strategies, which could be used for assessment of stability and adaptation of robot control behavior during pHRI. In future work, we aim to apply the method to more complex movements, inspired by realistic pHRI scenarios. Furthermore, we intend to complement our current apparatus by EMG sensing modalities, in order to obtain additional insights into the involved muscle activations.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

Ethics Statement

The studies involving human participants were reviewed and approved by Ethics committee at the Technical University of Munich. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author Contributions

HB and SE adapted the simulation and validated the method with simulated data. HB implemented and conducted the experiment, verified the method with experimental data, and wrote the first draft of the manuscript. SE and SH provided the critical revisions. All authors contributed to the conceptualization and methodology, read and approved the manuscript.

Funding

This research was supported by the European Union's Horizon 2020 research and innovation programme REHYB under grant agreement n° 87176 and the ERC Starting Grant Control based on Human Models (con-humo) under grant agreement n° 337654.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

1. ^Due to the inclusion of reflexive feedback, according to the terminology established in Latash and Zatsiorsky (1993), the matrices D_q and K_q represent apparent impedance components. In this work, in order to highlight the exclusion of voluntary feedback, we refer to them as involuntary impedance components. For simplicity, we refer to the individual matrices D_q and K_q as damping and stiffness instead of involuntary joint damping and stiffness.

References

Ajoudani, A., Fang, C., Tsagarakis, N. G., and Bicchi, A. (2015). “A reduced-complexity description of arm endpoint stiffness with applications to teleimpedance control,” in 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (Hamburg: IEEE), 1017–1023. doi: 10.1109/IROS.2015.7353495

CrossRef Full Text | Google Scholar

Artemiadis, P. K., Katsiaris, P. T., Liarokapis, M. V., and Kyriakopoulos, K. J. (2010). “Human arm impedance: characterization and modeling in 3d space,” in 2010 IEEE/RSJ International Conference on Intelligent Robots and Systems (Taipei: IEEE), 3103–3108. doi: 10.1109/IROS.2010.5652025

PubMed Abstract | CrossRef Full Text | Google Scholar

Asao, T., Suzuki, S., and Kotani, K. (2013). “Estimation of driver's steering intention by using mechanical impedance,” in International Conference on Human Interface and the Management of Information (Las Vegas: Springer), 3–11. doi: 10.1007/978-3-642-39209-2_1

CrossRef Full Text | Google Scholar

Bennett, D., Hollerbach, J., Xu, Y., and Hunter, I. (1992). Time-varying stiffness of human elbow joint during cyclic voluntary movement. Exp. Brain Res. 88, 433–442. doi: 10.1007/BF02259118

PubMed Abstract | CrossRef Full Text | Google Scholar

Börner, H., Endo, S., and Hirche, S. (2018). “Estimation of involuntary impedance in multi-joint arm movements,” in 2nd IFAC Conference on Cyber-Physical & Human-Systems (CPHS) (Miami, FL: IFAC).

Google Scholar

Brainard, D. H. (1997). The psychophysics toolbox. Spatial Vis. 10, 433–436. doi: 10.1163/156856897X00357

PubMed Abstract | CrossRef Full Text | Google Scholar

Burdet, E., Franklin, D. W., and Milner, T. E. (2013). Human Robotics: Neuromechanics and Motor Control. Cambridge, MA: MIT Press. doi: 10.7551/mitpress/9007.001.0001

PubMed Abstract | CrossRef Full Text | Google Scholar

Burdet, E., Osu, R., Franklin, D., Yoshioka, T., Milner, T., and Kawato, M. (2000). A method for measuring endpoint stiffness during multi-joint arm movements. J. Biomech. 33, 1705–1709. doi: 10.1016/S0021-9290(00)00142-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Burdet, E., Osu, R., Franklin, D. W., Milner, T. E., and Kawato, M. (2001). The central nervous system stabilizes unstable dynamics by learning optimal impedance. Nature 414, 446–449. doi: 10.1038/35106566

PubMed Abstract | CrossRef Full Text | Google Scholar

Burnham, K. P., and Anderson, D. R. (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. New York, NY: Springer.

Google Scholar

Buzzi, J., Ferrigno, G., Jansma, J. M., and De Momi, E. (2017). On the value of estimating human arm stiffness during virtual teleoperation with robotic manipulators. Front. Neurosci. 11:528. doi: 10.3389/fnins.2017.00528

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, P. H., Park, K., Kang, S. H., Krebs, H. I., and Hogan, N. (2013). Stochastic estimation of human arm impedance using robots with nonlinear frictions: an experimental validation. IEEE/ASME Trans. Mechatron. 18, 775–786. doi: 10.1109/TMECH.2012.2184767

PubMed Abstract | CrossRef Full Text | Google Scholar

Colgate, E., Bicchi, A., Peshkin, M. A., and Colgate, J. E. (2008). “Safety for physical human-robot interaction,” in Springer Handbook of Robotics (Heidelberg: Springer), 1335–1348. doi: 10.1007/978-3-540-30301-5_58

CrossRef Full Text | Google Scholar

Darainy, M., Malfait, N., Gribble, P. L., Towhidkhah, F., and Ostry, D. J. (2004). Learning to control arm stiffness under static conditions. J. Neurophysiol. 92, 3344–3350. doi: 10.1152/jn.00596.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Darainy, M., Towhidkhah, F., and Ostry, D. J. (2007). Control of hand impedance under static conditions and during reaching movement. J. Neurophysiol. 97, 2676–2685. doi: 10.1152/jn.01081.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, M., Inoue, A., Gomi, H., and Hirashima, Y. (2006). Recursive filter design for estimating time varying multijoint human arm viscoelasticity. Int. J. Comput. Syst. Signal 7, 2–18. Available online at: https://www.semanticscholar.org/paper/Recursive-Filter-Design-for-Estimating-Time-Varying-Deng-Inoue/7eec8b4245c500d6ae50460b6230342dcab1a2b5

Google Scholar

Dolan, J. M., Friedman, M. B., and Nagurka, M. L. (1993). Dynamic and loaded impedance components in the maintenance of human arm posture. IEEE Trans. Syst. Man Cybern. 23, 698–709. doi: 10.1109/21.256543

CrossRef Full Text | Google Scholar

Dyck, M., and Tavakoli, M. (2013). “Measuring the dynamic impedance of the human arm without a force sensor,” in 2013 IEEE International Conference on Rehabilitation Robotics (Seattle: IEEE), 1–8. doi: 10.1109/ICORR.2013.6650349

PubMed Abstract | CrossRef Full Text | Google Scholar

Erden, M. S., and Billard, A. (2015). Hand impedance measurements during interactive manual welding with a robot. IEEE Trans. Robot. 31, 168–179. doi: 10.1109/TRO.2014.2385212

PubMed Abstract | CrossRef Full Text | Google Scholar

Faisal, A. A., Selen, L. P., and Wolpert, D. M. (2008). Noise in the nervous system. Nat. Rev. Neurosci. 9, 292–303. doi: 10.1038/nrn2258

PubMed Abstract | CrossRef Full Text | Google Scholar

Flash, T., and Hogan, N. (1985). The coordination of arm movements: an experimentally confirmed mathematical model. J. Neurosci. 5, 1688–1703. doi: 10.1523/JNEUROSCI.05-07-01688.1985

PubMed Abstract | CrossRef Full Text | Google Scholar

Franklin, D. W., Burdet, E., Tee, K. P., Osu, R., Chew, C.-M., Milner, T. E., and Kawato, M. (2008). CNS learns stable, accurate, and efficient movements using a simple algorithm. J. Neurosci. 28, 11165–11173. doi: 10.1523/JNEUROSCI.3099-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Franklin, D. W., and Wolpert, D. M. (2011). Computational mechanisms of sensorimotor control. Neuron 72, 425–442. doi: 10.1016/j.neuron.2011.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Gautier, M., and Khalil, W. (1988). “On the identification of the inertial parameters of robots,” in 1988 IEEE Conference on Decision and Control (Austin: IEEE), 2264–2269. doi: 10.1109/CDC.1988.194738

CrossRef Full Text | Google Scholar

Gomi, H., and Kawato, M. (1997). Human arm stiffness and equilibrium-point trajectory during multi-joint movement. Biol. Cybern. 76, 163–171. doi: 10.1007/s004220050329

PubMed Abstract | CrossRef Full Text | Google Scholar

Gomi, H., and Osu, R. (1998). Task-dependent viscoelasticity of human multijoint arm and its spatial characteristics for interaction with environments. J. Neurosci. 18, 8965–8978. doi: 10.1523/JNEUROSCI.18-21-08965.1998

PubMed Abstract | CrossRef Full Text | Google Scholar

Haar, S., Donchin, O., and Dinstein, I. (2017). Individual movement variability magnitudes are explained by cortical neural variability. J. Neurosci. 37, 9076–9085. doi: 10.1523/JNEUROSCI.1650-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, C. M., and Wolpert, D. M. (1998). Signal-dependent noise determines motor planning. Nature 394:780. doi: 10.1038/29528

PubMed Abstract | CrossRef Full Text | Google Scholar

Hogan, N. (1985). The mechanics of multi-joint posture and movement control. Biol. Cybern. 52, 315–331. doi: 10.1007/BF00355754

PubMed Abstract | CrossRef Full Text | Google Scholar

Hondori, H. M., and Tech, A. W. (2011). “Smart mug to measure hand's geometrical mechanical impedance,” in 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society (Boston, MA: IEEE), 4053–4056. doi: 10.1109/IEMBS.2011.6091007

PubMed Abstract | CrossRef Full Text | Google Scholar

Horn, R. A., and Johnson, C. R. (2012). Matrix analysis. Cambridge: Cambridge University Press. doi: 10.1017/CBO9781139020411

CrossRef Full Text | Google Scholar

Khalil, W., and Dombre, E. (2004). Modeling, Identification and Control of Robots. Oxford, UK: Butterworth-Heinemann.

Google Scholar

Kim, H. K., Kang, B., Kim, B., and Park, S. (2009). Estimation of multijoint stiffness using electromyogram and artificial neural network. IEEE Trans. Syst. Man Cybern. A Syst. Hum. 39, 972–980. doi: 10.1109/TSMCA.2009.2025021

CrossRef Full Text | Google Scholar

Klodmann, J., Lakatos, D., Ott, C., and Albu-Schäffer, A. (2015). A closed-form approach to determine the base inertial parameters of complex structured robotic systems. IFAC Pap. Online 48, 316–321. doi: 10.1016/j.ifacol.2015.05.021

CrossRef Full Text | Google Scholar

Lakatos, D., Petit, F., and Van Der Smagt, P. (2011). “Conditioning vs. excitation time for estimating impedance parameters of the human arm,” in 2011 IEEE-RAS International Conference on Humanoid Robots (Bled: IEEE), 636–642. doi: 10.1109/Humanoids.2011.6100872

CrossRef Full Text | Google Scholar

Lakatos, D., Rüschen, D., Bayer, J., Vogel, J., and van der Smagt, P. (2013). “Identification of human limb stiffness in 5 dof and estimation via EMG,” in Experimental Robotics (Québec: Springer), 89–99. doi: 10.1007/978-3-319-00065-7_7

CrossRef Full Text | Google Scholar

Lametti, D. R., Houle, G., and Ostry, D. J. (2007). Control of movement variability and the regulation of limb impedance. J. Neurophysiol. 98, 3516–3524. doi: 10.1152/jn.00970.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Latash, M. L., and Zatsiorsky, V. M. (1993). Joint stiffness: myth or reality? Hum. Mov. Sci. 12, 653–692. doi: 10.1016/0167-9457(93)90010-M

CrossRef Full Text | Google Scholar

Masia, L., and Squeri, V. (2015). A modular mechatronic device for arm stiffness estimation in human-robot interaction. IEEE/ASME Trans. Mechatron. 20, 2053–2066. doi: 10.1109/TMECH.2014.2361925

CrossRef Full Text | Google Scholar

Matthews, P. B. (1991). The human stretch reflex and the motor cortex. Trends Neurosci. 14, 87–91. doi: 10.1016/0166-2236(91)90064-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehta, B., and Schaal, S. (2002). Forward models in visuomotor control. J. Neurophysiol. 88, 942–953. doi: 10.1152/jn.2002.88.2.942

PubMed Abstract | CrossRef Full Text | Google Scholar

Merton, P., and Morton, H. (1980). Stimulation of the cerebral cortex in the intact human subject. Nature 285:227. doi: 10.1038/285227a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mussa-Ivaldi, F. A., Hogan, N., and Bizzi, E. (1985). Neural, mechanical, and geometric factors subserving arm posture in humans. J. Neurosci. 5, 2732–2743. doi: 10.1523/JNEUROSCI.05-10-02732.1985

PubMed Abstract | CrossRef Full Text | Google Scholar

Neranon, P. (2020). Implicit force control approach for safe physical robot-to-human object handover. Indonesian J. Electric. Eng. Comput. Sci. 17, 615–628. doi: 10.11591/ijeecs.v17.i2.pp615-628

CrossRef Full Text | Google Scholar

Osu, R., and Gomi, H. (1999). Multijoint muscle regulation mechanisms examined by measured human arm stiffness and EMG signals. J. Neurophysiol. 81, 1458–1468. doi: 10.1152/jn.1999.81.4.1458

PubMed Abstract | CrossRef Full Text | Google Scholar

Palazzolo, J. J., Ferraro, M., Krebs, H. I., Lynch, D., Volpe, B. T., and Hogan, N. (2007). Stochastic estimation of arm mechanical impedance during robotic stroke rehabilitation. IEEE Trans. Neural Syst. Rehabil. Eng. 15, 94–103. doi: 10.1109/TNSRE.2007.891392

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, H., O'Neill, G., and Artemiadis, P. (2014). On the effect of muscular cocontraction on the 3-d human arm impedance. IEEE Trans. Biomed. Eng. 61, 2602–2608. doi: 10.1109/TBME.2014.2323938

PubMed Abstract | CrossRef Full Text | Google Scholar

Perreault, E. J., Kirsch, R. F., and Crago, P. E. (2004). Multijoint dynamics and postural stability of the human arm. Exp. Brain Res. 157, 507–517. doi: 10.1007/s00221-004-1864-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Piovesan, D., Pierobon, A., DiZio, P., and Lackner, J. R. (2013). Experimental measure of arm stiffness during single reaching movements with a time-frequency analysis. J. Neurophysiol. 110, 2484–2496. doi: 10.1152/jn.01013.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Roveda, L., Vicentini, F., and Tosatti, L. M. (2013). “Deformation-tracking impedance control in interaction with uncertain environments,” in 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems (Tokyo: IEEE), 1992–1997. doi: 10.1109/IROS.2013.6696621

CrossRef Full Text | Google Scholar

Sainburg, R. L. (2015). Should the equilibrium point hypothesis (eph) be considered a scientific theory? Motor Control 19, 142–148. doi: 10.1123/mc.2014-0056

PubMed Abstract | CrossRef Full Text | Google Scholar

Shadmehr, R., and Arbib, M. A. (1992). A mathematical analysis of the force-stiffness characteristics of muscles in control of a single joint system. Biol. Cybern. 66, 463–477. doi: 10.1007/BF00204111

PubMed Abstract | CrossRef Full Text | Google Scholar

Shin, D., Kim, J., and Koike, Y. (2009). A myokinetic arm model for estimating joint torque and stiffness from EMG signals during maintained posture. J. Neurophysiol. 101, 387–401. doi: 10.1152/jn.00584.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Slifkin, A. B., and Newell, K. M. (1999). Noise, information transmission, and force variability. J. Exp. Psychol. Hum. Percept. Perform. 25:837. doi: 10.1037/0096-1523.25.3.837

PubMed Abstract | CrossRef Full Text | Google Scholar

Söderström, T., and Stoica, P. (1989). System Identification. Hemel Hempstead, UK: Prentice-Hall.

Google Scholar

Stein, J. Y. (2000). Digital Signal Processing: A Computer Science Perspective. New York, NY: John Wiley & Sons, Inc. doi: 10.1002/047120059X

CrossRef Full Text | Google Scholar

Tee, K. P., Burdet, E., Chew, C.-M., and Milner, T. E. (2004). A model of force and impedance in human arm movements. Biol. Cybern. 90, 368–375. doi: 10.1007/s00422-004-0484-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Thilmann, A., Schwarz, M., Töpper, R., Fellows, S., and Noth, J. (1991). Different mechanisms underlie the long-latency stretch reflex response of active human muscle at different joints. J. Physiol. 444, 631–643. doi: 10.1113/jphysiol.1991.sp018898

PubMed Abstract | CrossRef Full Text | Google Scholar

Todorov, E., and Jordan, M. I. (2002). Optimal feedback control as a theory of motor coordination. Nat. Neurosci. 5:1226. doi: 10.1038/nn963

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomi, N., Gouko, M., and Ito, K. (2008). “Impedance control complements incomplete internal models under complex external dynamics,” in 2008 IEEE International Conference of the Engineering in Medicine and Biology Society (Vancouver, BC: IEEE), 5354–5357. doi: 10.1109/IEMBS.2008.4650424

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsuji, T., Morasso, P. G., Goto, K., and Ito, K. (1995). Human hand impedance characteristics during maintained posture. Biol. Cybern. 72, 475–485. doi: 10.1007/BF00199890

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, T., Dordevic, G. S., and Shadmehr, R. (2001). Learning the dynamics of reaching movements results in the modification of arm impedance and long-latency perturbation responses. Biol. Cybern. 85, 437–448. doi: 10.1007/s004220100277

PubMed Abstract | CrossRef Full Text | Google Scholar

Weitschat, R., Ehrensperger, J., Maier, M., and Aschemann, H. (2018). “Safe and efficient human-robot collaboration part I: estimation of human arm motions,” in 2018 IEEE International Conference on Robotics and Automation (ICRA) (Brisbane: IEEE), 1993–1999. doi: 10.1109/ICRA.2018.8461190

CrossRef Full Text | Google Scholar

Wong, J., Wilson, E. T., Malfait, N., and Gribble, P. L. (2009a). The influence of visual perturbations on the neural control of limb stiffness. J. Neurophysiol. 101, 246–257. doi: 10.1152/jn.90371.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, J., Wilson, E. T., Malfait, N., and Gribble, P. L. (2009b). Limb stiffness is modulated with spatial accuracy requirements during movement in the absence of destabilizing forces. J. Neurophysiol. 101, 1542–1549. doi: 10.1152/jn.91188.2008

PubMed Abstract | CrossRef Full Text | Google Scholar