Complex Stiffness Model of Physical Human-Robot Interaction: Implications for Control of Performance Augmentation Exoskeletons

TL;DR

This paper introduces a complex stiffness model for human joint dynamics in exoskeletons, incorporating hysteretic damping, validated through experiments, and demonstrates how this model can enhance control strategies for performance augmentation.

Contribution

The paper proposes a novel complex stiffness model with hysteretic damping for human joints, validated experimentally, and shows how it can improve exoskeleton control performance.

Findings

Hysteretic damping improves modeling accuracy over classical viscous damping.

A linear relationship exists between hysteretic damping and real stiffness.

A fractional order controller exploits hysteretic damping for better performance.

Abstract

Human joint dynamic stiffness plays an important role in the stability of performance augmentation exoskeletons. In this paper, we consider a new frequency domain model of the human joint dynamics which features a complex value stiffness. This complex stiffness consists of a real stiffness and a hysteretic damping. We use it to explain the dynamic behaviors of the human connected to the exoskeleton, in particular the observed non-zero low frequency phase shift and the near constant damping ratio of the resonant as stiffness and inertia vary. We validate this concept by experimenting with an elbow-joint exoskeleton testbed on a subject while modifying joint stiffness behavior, exoskeleton inertia, and strength augmentation gains. We compare three different models of elbow-joint dynamic stiffness: a model with real stiffness, viscous damping and inertia, a model with complex stiffness and inertia, and a model combining the previous two models. Our results show that the hysteretic damping term improves modeling accuracy, using a statistical F-test. Moreover this improvement is statistically more significant than using classical viscous damping term. In addition, we experimentally observe a linear relationship between the hysteretic damping and the real part of the stiffness which allows us to simplify the complex stiffness model as a 1-parameter system. Ultimately, we design a fractional order controller to demonstrate how human hysteretic damping behavior can be exploited to improve strength amplification performance while maintaining stability.

Tables

Table 1. Table I : Experiment Parameters

Exp	$\alpha$	Load (kg)	Grip (kg)	Bias (Nm)	Amplitude (Nm)	Frequency (rad/s)
I.1	$1$	$0.6$	$10$	$0$	$2\alpha$	$2-20$
I.2	$1$	$2.3$
I.3	$1$	$4.5$
I.4	$2$	$4.5$
I.5	$4$	$4.5$
II.1	$1$	$0.6$	$14$	$4\alpha$	$2\alpha$	$3-30$
II.2	$1$	$2.3$
II.3	$1$	$4.5$
II.4	$2$	$4.5$
II.5	$4$	$4.5$
III.1	$1$	$0.6$	$27$	$8\alpha$	$2\alpha$	$4-40$
III.2	$1$	$2.3$
III.3	$1$	$4.5$
III.4	$2$	$4.5$
III.5	$4$	$4.5$

Table 2. Table II : List of Symbols

Symbol	Meaning
${\tau }_d$	Actuator desired torque
${\tau }_s$	Actuator actual torque
$G_{SEA}$	Transfer function from ${\tau }_d$ to ${\tau }_s$
${\tau }_c$	Human-exoskeleton interaction torque
${\tau }_g$	Exoskeleton gravity torque
${\theta }_e$	Joint angular position
$K_h$	Human joint real stiffness parameter
$C_h$	Human joint hysteretic damping parameter
$B_h$	Human joint viscous damping parameter
$M_h$	Moment of inertia of human forearm
$M_e$	Moment of inertia of exoskeleton
$M_{h-e}$	Moment of inertia of human with exoskeleton
$M_{h-e/\alpha }$	Moment of inertia of human with attenuated exoskeleton
$S_h$	Joint dynamic stiffness of human
$S_{h-e}$	Joint dynamic stiffness of human with exoskeleton
$S_{h-e/\alpha }$	Joint dynamic stiffness of human with attenuated exoskeleton
$P_{\alpha },C_{\alpha }$	Plant and controller transfer functions of augmentation
${\omega }_{SEA}$	Natural frequency of $G_{SEA}$
${\omega }_h,{\zeta }_h$	Natural frequency and damping ratio of $S_h$
${\omega }_{h-e},{\zeta }_{h-e}$	Natural frequency and damping ratio of $S_{h-e}$
${\omega }_{h-e/\alpha },{\zeta }_{h-e/\alpha }$	Natural frequency and damping ratio of $S_{h-e/\alpha }$

Table 3. Table III : Subject Dynamic Stiffness Parameters

Exp	Model	$K_h\left(\frac{Nm}{rad}\right)$	$C_h\left(\frac{Nm}{rad}\right)$	$B_h\left(\frac{Nms}{rad}\right)$	$M_{h-e/\alpha }$$\left(kg{m}^2\right)$	${\omega }_{h-e/\alpha }\left(\frac{rad}{s}\right)$	${\zeta }_{h-e/\alpha }$
I.1	M1	$10.05$	$--$	$1.03$	$0.28$	$5.95$	$0.31$
	M2	$10.05$	$5.89$	$--$	$0.28$	$5.95$	$0.29$
	M3	$10.05$	$4.97$	$0.18$	$0.28$	$5.95$	$0.30$
5pt. I.2	M1	$11.80$	$--$	$1.51$	$0.60$	$4.44$	$0.28$
	M2	$11.80$	$6.68$	$--$	$0.60$	$4.44$	$0.28$
	M3	$11.80$	$5.44$	$0.31$	$0.60$	$4.44$	$0.29$
5pt. I.3	M1	$15.74$	$--$	$2.09$	$1.18$	$3.65$	$0.24$
	M2	$15.74$	$8.33$	$--$	$1.18$	$3.65$	$0.26$
	M3	$15.74$	$10.44$	$-0.60$	$1.18$	$3.65$	$0.26$
5pt. I.4	M1	$13.82$	$--$	$1.46$	$0.60$	$4.78$	$0.25$
	M2	$13.82$	$6.87$	$--$	$0.60$	$4.78$	$0.25$
	M3	$13.82$	$6.01$	$0.21$	$0.60$	$4.78$	$0.25$
5pt. I.5	M1	$12.09$	$--$	$1.22$	$0.28$	$6.59$	$0.33$
	M2	$12.09$	$6.84$	$--$	$0.28$	$6.59$	$0.28$
	M3	$12.09$	$4.26$	$0.52$	$0.28$	$6.59$	$0.32$
II.1	M1	$12.73$	$--$	$1.41$	$0.20$	$7.94$	$0.44$
	M2	$12.73$	$10.18$	$--$	$0.20$	$7.94$	$0.40$
	M3	$12.73$	$5.86$	$0.66$	$0.20$	$7.94$	$0.44$
5pt. II.2	M1	$18.79$	$--$	$1.91$	$0.57$	$5.72$	$0.29$
	M2	$18.79$	$11.77$	$--$	$0.57$	$5.72$	$0.31$
	M3	$18.79$	$11.54$	$0.04$	$0.57$	$5.72$	$0.31$
5pt. II.3	M1	$25.95$	$--$	$3.08$	$1.03$	$5.02$	$0.30$
	M2	$25.95$	$16.75$	$--$	$1.03$	$5.02$	$0.32$
	M3	$25.95$	$15.48$	$0.26$	$1.03$	$5.02$	$0.32$
5pt. II.4	M1	$25.77$	$--$	$2.83$	$0.52$	$7.02$	$0.39$
	M2	$25.77$	$20.49$	$--$	$0.52$	$7.02$	$0.40$
	M3	$25.77$	$16.60$	$0.60$	$0.52$	$7.02$	$0.40$
5pt. II.5	M1	$19.07$	$--$	$1.88$	$0.28$	$8.32$	$0.41$
	M2	$19.07$	$16.27$	$--$	$0.28$	$8.32$	$0.43$
	M3	$19.07$	$15.72$	$0.08$	$0.28$	$8.32$	$0.43$
III.1	M1	$48.15$	$--$	$1.97$	$0.23$	$14.4$	$0.29$
	M2	$48.15$	$25.45$	$--$	$0.23$	$14.4$	$0.26$
	M3	$48.15$	$16.66$	$0.76$	$0.23$	$14.4$	$0.29$
5pt. III.2	M1	$48.60$	$--$	$2.85$	$0.58$	$9.13$	$0.27$
	M2	$48.60$	$25.61$	$--$	$0.58$	$9.13$	$0.26$
	M3	$48.60$	$15.19$	$1.23$	$0.58$	$9.13$	$0.27$
5pt. III.3	M1	$42.23$	$--$	$3.19$	$1.01$	$6.47$	$0.24$
	M2	$42.23$	$23.60$	$--$	$1.01$	$6.47$	$0.28$
	M3	$42.23$	$24.08$	$-0.07$	$1.01$	$6.47$	$0.28$
5pt. III.4	M1	$32.22$	$--$	$2.82$	$0.46$	$8.35$	$0.37$
	M2	$32.22$	$25.36$	$--$	$0.46$	$8.35$	$0.39$
	M3	$32.22$	$20.83$	$0.55$	$0.46$	$8.35$	$0.39$
5pt. III.5	M1	$42.33$	$--$	$2.08$	$0.27$	$12.43$	$0.31$
	M2	$42.33$	$26.50$	$--$	$0.27$	$12.43$	$0.31$
	M3	$42.33$	$27.66$	$-0.11$	$0.27$	$12.43$	$0.31$

Equations

$$S_h={\tau }_c/{\theta }_e=M_h{s}^2+B_{h}s+K_h.$$

$$S_h={\tau }_c/{\theta }_e=M_h{s}^2+C_{h}j+K_h.$$

$$S_h={\tau }_c/{\theta }_e=M_h{s}^2+B_{h}s+C_{h}j+K_h.$$

$$S_{h\text{mathchar}45e\text{mathchar}8239\alpha }=M_{h\text{mathchar}45e\text{mathchar}8239\alpha }{s}^2+B_{h}s+K_h,$$

$$S_{h\text{mathchar}45e\text{mathchar}8239\alpha }=M_{h\text{mathchar}45e\text{mathchar}8239\alpha }{s}^2+C_{h}j+K_h,$$

$$S_{h\text{mathchar}45e\text{mathchar}8239\alpha }=M_{h\text{mathchar}45e\text{mathchar}8239\alpha }{s}^2+B_{h}s+C_{h}j+K_h,$$

$${\zeta }_{h\text{mathchar}45e\text{mathchar}8239\alpha }=\frac{B_h}{2\sqrt{K_h{M}_{h\text{mathchar}45e\text{mathchar}8239\alpha }}}\text{for M1,}$$

$${\zeta }_{h\text{mathchar}45e\text{mathchar}8239\alpha }=\frac{C_h{\omega }_{h\text{mathchar}45e\text{mathchar}8239\alpha }^{-1}}{2\sqrt{K_h{M}_{h\text{mathchar}45e\text{mathchar}8239\alpha }}}=\frac{C_h}{2K_h}\text{for M2, and}$$

$${\zeta }_{h\text{mathchar}45e\text{mathchar}8239\alpha }=\frac{B_h}{2\sqrt{K_h{M}_{h\text{mathchar}45e\text{mathchar}8239\alpha }}}+\frac{C_h}{2K_h}\text{for M3,}$$

$$F_{Mi\text{mathchar}45M3}=\frac{RS{S}_{Mi}-RS{S}_{M3}}{RS{S}_{M3}}(2n-4),\text{for}i=1,2$$

$$={\tan }^{-1}(\frac{C_h}{K_h})={\tan }^{-1}(c_h+\frac{d_h}{K_h})\text{for M2, and}$$

$$={\tan }^{-1}(c_h+\frac{d_h+B_h\omega }{K_h})\mathrm{for}M3,$$

$${\zeta }_{h\text{mathchar}45e\text{mathchar}8239\alpha }$$

$$S_{h\text{mathchar}45e\text{mathchar}8239\alpha }=M_{h\text{mathchar}45e\text{mathchar}8239\alpha }{s}^2+K_h(1+a_{h}s).$$

$${\zeta }_{h\text{mathchar}45e\text{mathchar}8239\alpha }=\frac{a_h{K}_h}{2\sqrt{K_h{M}_{h\text{mathchar}45e\text{mathchar}8239\alpha }}}=\frac{a_h}{2}\cdot {\omega }_{h\text{mathchar}45e\text{mathchar}8239\alpha },$$

$$S_{h\text{mathchar}45e\text{mathchar}8239\alpha }=M_{h\text{mathchar}45e\text{mathchar}8239\alpha }{s}^2+K_h(1+c_{h}j).$$

$$S_h={\tau }_c/{\theta }_e=M_h{s}^2+K_h(1+c_{h}j),$$

$$S_{h\text{mathchar}45e}={\tau }_s/{\theta }_e=M_{h\text{mathchar}45e}{s}^2+K_h(1+c_{h}j),$$

$$\frac{{\tau }_{\alpha }}{{\tau }_s}=\frac{(\alpha -1)\cdot S_h+S_{h\text{mathchar}45e}}{S_{h\text{mathchar}45e}}=\alpha \cdot \frac{S_{h\text{mathchar}45e\text{mathchar}8239\alpha }}{S_{h\text{mathchar}45e}}.$$

$$P_{\alpha }(s)=\frac{{\tau }_{\alpha }}{{\tau }_d}=\alpha \cdot \frac{S_{h\text{mathchar}45e\text{mathchar}8239\alpha }}{S_{h\text{mathchar}45e}}\cdot G_{SEA}(s),$$

$$C_{\alpha }(s)=k_f/s^f,$$

$$0<90f$$

$$0<f$$

$$C_{\alpha }(s)$$

$$z_i/p_i$$

$$p_i/p_{i-1}$$

Full text

Accepted for publication in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) ©2019 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. DOI: 10.1109/IROS40897.2019.8968005

I Introduction

While the concept of a personal augmentation device or exoskeleton is an old idea [1, 2, 3], a system which delivers on the dream of transparent interaction, of “feeling like the system is not there,” through augmentation of sensed human interaction forces is still an ambitious goal of force control technology today [4, 5, 6, 7]. Unlike assistive exoskeletons which help complete predictable behaviors [8, 9] or rehabilitation exoskeletons [10, 11] which simulate rehabilitation therapy, human augmentation exoskeletons [4, 12] use non-passive feedback control to amplify the user’s strength. But this type of feedback control brings the system closer to instability. And since the exoskeleton is in a feedback interconnection with the human, a model of the human’s dynamic behavior plays a critical role in determining the stability of an augmentation exoskeleton [13, 14].

Among all different kinds of dynamic model of an individual human joint, perhaps the most popular one is the mass-spring-damper model—with the additional non-linearity that the spring stiffness of the human joint can be modified by both voluntary muscle contractions or external torques exerted on the joint [15]. Several studies demonstrated a linear relationship between the stiffness of the human (found by fitting a linear mass-spring-damper model for a single joint) and an external torque [16, 17, 18]. For modeling the human joint damping, some other studies explored the fact that not only the stiffness but also the damping increases with muscle contractions [19] and external torques [20]. A linear relationship between the damping and the external torque has also been identified for the human ankle joint, but it is statistically weaker than the strong linear relationship between the stiffness of the ankle and the external torques [16, 18]. However, it is not clear from the literature that a linear relationship between the damping and the stiffness of a human joint can be expected in more general cases.

Another way to model the damping in the linear mass-spring-damper model is through the empirical observation that a relatively consistent damping ratio is maintained by the human elbow across different joint stiffnesses [14]. Frequency domain identification of the ankle joint impedance [16, 21] also showed a consistent damping ratio within the range from 0.22 to 0.49. This damping ratio consistency on the ankle is also supported by the fact that the ankle damping ratio does not have significant change with large variations of mean external torques exerted on the subjects [20]. For upper limbs, a multi-joint impedance study on human arms [22] showed that the damping ratio of the minimally damped mode for the 2-D endpoint impedance in the transverse plane is distributed with a mean of 0.26 and a standard deviation of 0.08. Although this could be explained as the effect of humans adapting their damping to stabilize movement [23], a more detailed explanation of how humans achieve this consistency remains unclear.

Hysteretic damping models have seen success in biomechanical modelling before. In [16], experimental results showed a hysteretic relationship between the applied torque and the ankle angle at very low frequencies. Hysteretic damping is shown indirectly in [17] (see Fig. 6 of that paper), where the human elbow stiffness has a phase shift around 25 degrees in a wide range of low frequencies—contradicting the viscous damping hypothesis. This type of phase behavior is explained (in the field of structural mechanics) by defining a hysteretic damping whose damping coefficient is proportional to the inverse of frequency [24]. Models with hysteretic damping have also been adapted to describe the dynamic properties of the whole body of a seated human [25] as well as cockroach legs [26].

In this paper we study the human stiffness and damping behavior when coupled to an exoskeleton inertia, and test the effectiveness of a hysteretic damping term in the system model. More specifically we compare three models 1) a linear mass, spring, and viscous damper model, 2) a nonlinear complex-stiffness-spring and mass model (that is, a spring, mass, and hysteretic damper model), and 3) a combination model with mass, spring, and both viscous and hysteretic damping. Our results show that there is a statistically significant benefit of the hysteretic damping term (comparing model 1 to model 3 with an F-test), and a less significant benefit for the viscous damping term (comparing model 2 to model 3). This hysteretic damping explains the consistent damping-ratio of the human–exoskeleton resonant peak even as the stiffness and exoskeleton inertia change—which is not well explained by the linear model. And it also explains the low frequency phase lag in human stiffness (previously observed in [17]). Our elbow joint experiments vary parameters which would result in a differing damping ratio if the linear model were true: we change the inertia of the exoskeleton, and (indirectly, using an adjustable exercise hand grip and a bias torque) the stiffness of the human joint. We also test different exoskeleton strength amplification factors, and it does not appear to elicit a different human behavior than when the inertia is simply reduced. One further contribution of the paper is the theorizing of an amplification controller which uses fractional order filtering to exploit the hysteretic damping of the human, offering improved performance over previous strategies.

II Methods

II-A Apparatus

For this study we employed the P0 series elastic elbow-joint exoskeleton from Apptronik Systems, as shown in Fig. 1. This exoskeleton has a moment of inertia of 0.1 $\rm{kg\cdot m^{2}}$ with no load on it, but allows for attaching additional weights to it. A load, attached 0.45 m from the exoskeleton joint, is pictured in Fig. 1.b. The contact force $f_{c}$ between the human and the exoskeleton is measured by a six-axis force/torque sensor situated below the white 3D printed “cuff” (which includes the adjustable strap which clamps the forearm). This force torque signal is cast as a torque ($\tau_{c}$) using the motion Jacobian $J$ of the sensor frame ($\tau_{c}=J^{T}f_{c}$). Rubber pads are adhered to the inside surfaces of the cuff and the cuff strap to improve user comfort. Joint position $\theta_{e}$ is directly measured by a dedicated encoder at the exoskeleton joint. The series elastic actuator (SEA) has a spring force control bandwidth of 10 Hz and provides high fidelity actuator torque $\tau_{s}$ tracking using the force control disturbance observer of [27].

In parallel with an excitation chirp command (which essentially performs system identification of the human subject), a gravity compensation controller, a human augmentation controller, and a bias torque comprise the desired actuator torque signal. The gravity compensation controller takes the measurement of $\theta_{e}$ to calculate and compensate the gravity torque $\tau_{g}$ acting on the exoskeleton system. The human augmentation controller takes the measurement of $\tau_{c}$ and multiplies $\tau_{c}$ by negative $\alpha-1$. With the assistance of actuator torques produced from the augmentation command, the human’s interaction forces with the exoskeleton are amplified by a factor of $\alpha$. This exoskeleton augmentation strategy differs from the one we applied in [14] in the directness of the augmentation feedback.

II-B Experimental Protocol

The experimental protocol was approved by the Institutional Review Board (IRB) at the University of Texas at Austin. It consists of fifteen perturbation experiments with a 28-year old male subject. The experiments are separated into three groups (Exp. I-III) of five experiments. The first three experiments in each group are conducted with loads of 0.6 kg, 2.3 kg and 4.5 kg and an $\alpha$ value of 1 (corresponding to a non-augmentation controller) while the last two experiment in each group are conducted with a load of 4.5 kg and $\alpha$ values of 2 and 4. The mass of the loads and the mass of the exoskeleton have their gravitational bias torque fully compensated through gravity compensation control, while their inertia is attenuated by a factor of $\alpha$ due to the cuff torque feedback.

The stiffness of the human elbow is influenced by muscle co-contraction as well as by contraction to resist the bias torque. In order to obtain different values of elbow stiffness for the three experiment groups, both the bias torque component of the controller and the co-contraction are varied. The three groups have, respectively, [math] Nm, $4\alpha$ Nm, and $8\alpha$ Nm of bias torque. Co-contraction is controlled by having the subject squeeze an adjustable force hand grip. The three groups have a 10-kg, a 14-kg, and a 27-kg gripping force respectively. The amplitude of the perturbation chirp signal is set to be $2\alpha$ Nm.

To avoid fatigue of the subject, the duration of each perturbation experiment is set to be 100 seconds. The perturbation is set to be an exponential chirp signal, and the results are typically analyzed in the frequency domain. To sufficiently capture the natural frequency for damping feature identification, we set different ranges of frequency for the chirp signal according to the stiffness values the subject achieved from the bias torque and the gripping force. Frequency ranges of 2-20 rad/s, 3-30 rad/s and 4-40 rad/s are set for the chirp signals for the three experiment groups.

After the chirp perturbation experiments, we transfer the time domain data into the frequency domain and identify the dynamic stiffness model of the subject by linear regression. The parameters of the three experiment groups are summarized in Tab. I.

II-C Models

In our models, we define $K_{h}$ as the human elbow-joint real stiffness, $C_{h}$ as the human elbow-joint hysteretic damping, $B_{h}$ as the human elbow-joint viscous damping, $M_{h}$ as the moment of inertia of the human and $M_{e}$ as the moment of inertia of the exoskeleton. See list of symbols in Tab. II.

A passive linear model of human dynamic stiffness with viscous damping can be expressed as

[TABLE]

If we consider a human model with hysteretic damping (complex stiffness) we have a nonlinear model

[TABLE]

And to generalize the two, we also consider a nonlinear model with both viscous and hysteretic damping

[TABLE]

However, these models are difficult to identify from the experimental $\tau_{c}$ and $\theta_{e}$ values because the natural frequency of the human dynamic stiffness can easily go beyond the range of the frequency for the experiments. With the augmentation controller, the operator feels an attenuated inertia from the exoskeleton. Therefore, we added a nominal attenuated inertia of $M_{e}/\alpha$ to the frequency domain data of $\tau_{c}/\theta_{e}$ for the model identification. In essence, we desensitize our identification to errors far above the natural frequency of the human spring and the exoskeleton inertia. Combining this additional term with (1), (2) and (3), the three models of human-exoskeleton interaction can be expressed as

[TABLE]

where $M_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}=M_{h}+M_{e}/\alpha$ is the perceived inertia at the human joint.

We also calculate the damping ratio $\zeta_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ of $S_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$, as a measure of the degree of oscillation at the resonant zero-pair. Because M2 and M3 have the $C_{h}j$ term which provides a damping effect in addition to $B_{h}s$, we define the damping ratio of each model using the imaginary part of the transfer function evaluated at the resonance:

[TABLE]

where $\omega_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}=\sqrt{K_{h}/M_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}}$ is the natural frequency of $S_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$.

II-D Statistical Analysis

In order to compare the significance of $B_{h}s$ and $C_{h}j$ in the human-exoskeleton interaction model, we calculate the residual square sum (RSS) for all three models, denoted $RSS_{M1}$, $RSS_{M2}$ and $RSS_{M3}$ respectively. For each experiment, we conduct F-tests for each of the two three-parameter models (M1 and M2) against the generalizing four-parameter model (M3). Our F-statistic accounts for complex value data,

[TABLE]

where $n$ is the number of complex value samples at the frequency domain and the real and imaginary parts of each sample are statistically independent. The significance of $B_{h}s$ and $C_{h}j$ then will be evaluated by comparing this F statistic against a critical F statistic threshold based on a 0.05 false-rejection probability.

We split the 100 seconds of time domain data for each experiment into 10 sequences. For each of the 10 second sequences, only the data from the first 5.78 seconds is used for calculating the frequency domain sample. The remainder period of 4.22 seconds is greater than the 2% settling time for all the 2nd order dynamics of $S_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ identified in the experiments. By this method we can safely assume statistical independence between the 10 single-frequency data points comprising our estimate of the frequency response function for the purposes of statistical testing.

III Results

III-A Phase Shift

In the frequency domain results of $S_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ (Fig. 3), the phase starts (at low frequencies) from a value between $25\degree$ to $45\degree$ instead of zero and changes very little across all the frequencies before it reaches the second order zero at $\omega_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ for each experiment. This type of phase shift is very different from the phase shift usually experienced by a linear system with a constant time delay or a constant damping property in which the phase shift approaches zero in the limit as $\omega\rightarrow 0$. As shown in Fig. 4, this phase shift is clearly visible even in time domain comparisons of $\tau_{c}$ and $\theta_{e}$. The data show that the human joint motion $\theta_{e}$ is not perfectly sinusoidal—it stops following the trend of the torque after they both reach their peak values and “waits” before following the torque $\tau_{c}$ in its descent. At low frequencies, these peaks seem especially flat.

III-B Model Comparisons

The results of the identified parameters (Tab. III-C) show that the three models give the same values of $K_{h}$, $M_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ and consequently $\omega_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ to two decimal places for each experiment. This is because the difference between the three models is restricted to the imaginary part of $S_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ while $K_{h}$ and $M_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ are the coefficients of the real part of $S_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$. Although the identified values of $B_{h}$ and $C_{h}$ are quite different between the three models, the values of $\zeta_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ are still very close for each experiment. This means that the three models give very similar values for the slope of the phase at the resonant frequency $\omega_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$.

From M1 to M3, the values of $B_{h}$ have been reduced considerably. This means that M3 uses the $C_{h}j$ term to replace part of the $B_{h}$ term in M1 while maintaining a similar phase behavior at the frequency $\omega_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$. From M2 to M3, the values of $C_{h}$ have been reduced except for Exp. I.3, III.3 and III.5 in which M3 gives a negative value for $B_{h}$. These negative value of $B_{h}$ is because there is no lower bound constraint on the value of $B_{h}$ during the frequency domain regression for M3. Although a negative value of $B_{h}$ brings non-passivity to a linear mass-spring-damper system in the common sense, the $C_{h}j$ term in M3 enforces the dynamics of $S_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$ to remain passive across the range of frequencies in our experiments.

The results from the F-tests (Fig. 5) relate to the significance of $B_{h}s$ and $C_{h}j$ in M3. Based on the 20 statistically independent data values for each experiment, a critical F-statistic value of 4.49 is calculated for 0.05 false-rejection probability. The results show that values of $F_{M1\mathchar 45\relax M3}$ for all the experiments are much higher than the critical F-statistic value, with the values of $F_{M1\mathchar 45\relax M3}$ in Exp. II.3 and II.5 exceeding 100 (c.f. the critical value of 4.49). This proves that the existence of the $C_{h}j$ term in M3 significantly improves modeling accuracy of $S_{h\mathchar 45\relax e\mathchar 8239\relax\alpha}$. The values of $F_{M2\mathchar 45\relax M3}$ are mostly below the critical F-statistic value except for Exp. I.5, II.1, III.1 and III.2. The other observation is that the value of $F_{M2\mathchar 45\relax M3}$ is always much lower than the value of $F_{M2\mathchar 45\relax M3}$ for all experiments. Although the effect of the $B_{h}s$ term cannot be completely ignored based on the results of these F-tests, we can claim that the $C_{h}j$ term is still much more significant than the $B_{h}s$ term in M3.

III-C Linear Regression between $C_{h}$ and $K_{h}$

Because the $C_{h}j$ term is created to describe the phase shift effect from the complex human stiffness in M2 and M3, we suspect that the identified value of $C_{h}$ has a linear relation with the value of $K_{h}$. Therefore, we apply linear regression between the values of $C_{h}$ and $K_{h}$ identified from M2 and M3 (Fig. 6). Compared with M3, the linear regression result with M2 shows a stronger linear relationship with a much higher coefficient of determination ($R^{2}$). The regression equation identified from the M2 parameters also has a smaller value of bias from the origin of the $C_{h}\mathchar 45\relax K_{h}$ plane compared with the regression equation identified from the M3 parameters. Intuition leads us to expect low bias in the regression equation, since a nonzero value of $C_{h}$ when the value of $K_{h}$ is zero could not be explained as hysteretic spring behavior.