A Complex Stiffness Human Impedance Model with Customizable Exoskeleton Control
LLOGO
IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. XX, NO. XX, XXXX 20XX 1
A Complex Stiffness Human Impedance Modelwith Customizable Exoskeleton Control
Binghan He, Huang Huang, Gray C. Thomas,
Member, IEEE , and Luis Sentis,
Member, IEEE
Accepted for publication in Transactions on Neural Systems and Rehabilitation Engineering (TNSRE) ©2020 IEEE. Personal use of this material is permitted. Permission fromIEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating newcollective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. DOI: 10.1109/TNSRE.2020.3027501
Abstract — The natural impedance, or dynamic relation-ship between force and motion, of a human operator candetermine the stability of exoskeletons that use interaction-torque feedback to amplify human strength. While humanimpedance is typically modelled as a linear system, ourexperiments on a single-joint exoskeleton testbed involving10 human subjects show evidence of nonlinear behavior: alow-frequency asymptotic phase for the dynamic stiffnessof the human that is different than the expected zero, and anunexpectedly consistent damping ratio as the stiffness andinertia vary. To explain these observations, this paper con-siders a new frequency-domain model of the human jointdynamics featuring complex value stiffness comprising areal stiffness term and a hysteretic damping term. Using astatistical F-test we show that the hysteretic damping termis not only significant but is even more significant thanthe linear damping term. Further analysis reveals a lineartrend linking hysteretic damping and the real part of thestiffness, which allows us to simplify the complex stiffnessmodel down to a 1-parameter system. Then, we introduceand demonstrate a customizable fractional-order controllerthat exploits this hysteretic damping behavior to improvestrength amplification bandwidth while maintaining stabil-ity, and explore a tuning approach which ensures that thisstability property is robust to muscle co-contraction foreach individual.
Index Terms — Human impedance, human performanceaugmentation, exoskeletons.
I. I
NTRODUCTION W HILE the concept of a personal augmentation deviceor exoskeleton has a long history [1], [2], [3], a systemwhich delivers on the dream of transparent interaction, of“feeling like the system is not there,” through amplificationof sensed human interaction forces is still an ambitious goalof force control technology today [4], [5], [6], [7]. Unlike otherassistive exoskeletons that help perform predictable behaviors
This work was supported by the U.S. Government and NASA SpaceTechnology Research Fellowship NNX15AQ33H. The authors thank themembers of the Human Centered Robotics Lab, The University of Texasat Austin, and Apptronik Systems for their support. (Correspondingauthor: Binghan He.)
Binghan He and Huang Huang are with the Department of MechanicalEngineering, The University of Texas at Austin, Austin, TX 78712 USA(e-mail: [email protected]; [email protected])Gray C. Thomas is with the Department of Electrical and ComputerEngineering, University of Michigan, Ann Arbor, MI 48109 USA (e-mail:[email protected])Luis Sentis is with the Department of Aerospace Engineering andEngineering Mechanics, The University of Texas at Austin, Austin, TX78712 USA (e-mail: [email protected]) [8], [9], provide rehabilitation therapy [10], [11], or adjustnatural dynamics in a helpful way [12], human amplificationexoskeletons [4], [13] assist users, whether patients or healthypeople, by amplifying their strength (and power) throughfeedback control. But this type of feedback control introducesa risk of instability (a risk that was first noted in the fieldof impedance control for physical human–robot interactions[14]). Since the exoskeleton is in a feedback interconnectionwith the human, a model of the human’s dynamic behaviorplays a critical role in determining the stability of the closedloop amplification exoskeleton system [15], [16].Time and gait-phase varying models of the human jointimpedance [17], [18], [19], [20], [21] have been pursued bythe bio-mechanics and wearable robotics communities to betterreplicate human lower-limb behavior. Among all differentkinds of impedance model of an individual human joint, per-haps the most popular one is the mass-spring-damper model—with the additional non-linearity that the spring stiffness ofthe human joint can be modified by both voluntary musclecontractions or external torques exerted on the joint [22].Several studies demonstrated a linear relationship between thestiffness (not to be confused with quasi-stiffness [23]) of thehuman (found by fitting a linear mass-spring-damper modelfor a single joint) and an external torque [24], [25], [26],[27], [28]. Joint damping has also been shown to increase withmuscle contractions [29] and external torques [30]. A linearrelationship between damping and external torque has alsobeen reported for the same human joints, but it is statisticallyweaker than the strong linear relationship between stiffnessand external torques [24], [28]. However, it is not clear fromthe literature that a linear relationship between damping andstiffness in human joints does exist.Yet inconsistencies in the variable mass-spring-dampermodel remain [31], and the empirical observation that arelatively consistent damping ratio is maintained in some joints(notably the human elbow [15] and arm [32]) even as jointstiffness and inertia vary is one such anomaly. Frequencydomain identification of the ankle joint impedance [24], [25]also shows a consistent damping ratio within the range from . to . . This damping ratio consistency on the ankle isalso supported by the fact that the ankle damping ratio does nothave significant change with large variations of mean externaltorques exerted on the subjects [30]. For upper limbs, a multi-joint impedance study on human arms [32] showed that thedamping ratio of the minimally damped mode for the 2-D a r X i v : . [ c s . R O ] S e p IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. XX, NO. XX, XXXX 20XX endpoint impedance in the transverse plane is distributed witha mean of . and a standard deviation of . . Althoughthis could be explained as humans adapting their damping tostabilize movement [33], a more detailed explanation of howhumans achieve this consistency remains unclear.One potential solution is to use hysteretic damping mod-els. Hysteretic damping can correctly explain the behaviorsshown in [24], [25], [26], [27], where the phase plots ofthe human stiffness have non-zero values at low frequencies.In particular, Ref. [27, Fig. 6] shows that the human elbowdynamic stiffness has a phase shift around ° for a wide rangeof low frequencies, thus contradicting the viscous dampinghypothesis. This type of phase behavior is explained in thefield of structural mechanics by defining a hysteretic dampingterm whose damping coefficient is proportional to the inverseof the frequency [34].This paper introduces and validates a complex stiffnessmodel for human elbow joint dynamics. The primary modelvalidation experiment uses statistical F-tests to compare threedynamic stiffness models: a linear mass, spring, and viscousdamper model, a nonlinear complex-stiffness-spring and massmodel (that is, a spring, mass, and hysteretic damper model),and a combined model with mass, spring, and both viscousand hysteretic damping. This hysteretic damping explains theconsistent damping ratio of the human–exoskeleton resonantpeak even as the stiffness and exoskeleton inertia change—which is not well explained by the linear model. And italso explains the low frequency phase lag in human stiffness(previously observed in [24], [25], [26], [27]).Using this new model, this paper introduces a customizablefractional-order controller designed to take full advantage ofthe low-frequency phase lag for each individual. Based onresults from the previous test, a customized fractional orderis chosen for each of three subjects such that the behavioris nearly oscillatory (marginally stable). The subjects thenchange their co-contraction level to illustrate the phenomenonof co-contraction induced instability and subject dependent co-contraction relationships with stability. The three subjects spanthe range of observed co-contraction–stability relationships.This paper builds significantly on our earlier conferencepresentation [35]. First, this study investigates the applicabilityof the model in the more general population ( N = 10 ),whereas [35] only supported the model in one single individual(
N = 1 ). Second, this study proposes a novel power lawrelationship that is necessary to describe the range of behav-iors observed in the wider population, whereas [35] relatedstiffness and hysteretic damping with a naive linear model.Third, this study highlights the novel finding that there arefundamental differences in the stiffness–hysteretic dampingrelationships between subjects, and the critical influence thishas on the problem of designing tests to pre-certify safety inexoskeleton controllers, whereas [35] had no inter-subject data.Finally, this study employs a physical implementation, tests thethree most extreme subjects, and provides the first empiricalvalidation that the combination of fractional order control andtuning to the hysteretic damping model can improve dynamicalamplification (at
10 rad / s ) by ∼ , whereas [35]only proposed the concept of using fractional-order controllers (a) (b)Apparatus Experimental Setup6-Axis Force SensorSEAEncoder Hand GripLoadSpring Trigger Fig. 1. Experimental apparatus: the series elastic P0 exoskeletonfeaturing an ATI Mini40 force sensitive cuff and a P170 Orion air cooledseries elastic actuator module acting through a simple 3 bar linkage.During all experiments, subjects apply forces to a adjustable handgrip to regulate their elbow stiffness. A spring trigger is only used forperturbation during the loop shaping experiments in Sec. IV-C. Se τ s θ e Se τ g I M e τ c B h I M h C K h + C h j α − G SEA ( s ) α τ d Bias Amplification − Gravity Compensation − +Excitation Fig. 2. Block diagram consisting of amplification, gravity compensationand experimental perturbation. The dynamics of the human with theexoskeleton are expressed as a bond graph with the effort sources τ s , τ c and τ g . to exploit the hysteretic damping characteristics without anyexperimentation.The experimental protocol for model identification andcontroller implementation was approved by the InstitutionalReview Board (IRB) at the University of Texas at Austin understudy No . . The informed consent forms weresigned by all subjects. II. M
ODELING M ETHODS
A. Apparatus
For this study we employed the P0 series elastic elbow-joint exoskeleton (Apptronik Systems, Inc., Austin, TX), asshown in Fig. 1. This exoskeleton has a moment of inertiaof 0.1 kg · m with no load on it but allows for attachingadditional weights to it. A load, attached 0.45 m from theexoskeleton joint, is pictured in Fig. 1.(b). The contact force f c between the human and the exoskeleton is measured bya six-axis force/torque sensor situated below the white 3Dprinted “cuff” (which includes the adjustable strap that clampsthe forearm). This force torque signal is cast as a torque ( τ c )using the motion Jacobian J of the sensor frame ( τ c = J T f c ).Rubber pads are adhered to the inside surfaces of the cuff andthe cuff strap to improve user comfort. Joint position θ e isdirectly measured by a dedicated encoder at the exoskeletonjoint. The series elastic actuator (SEA) has a spring forcecontrol bandwidth of
10 Hz and provides high fidelity actuatortorque τ s tracking using the disturbance observer of [36].In parallel with an excitation chirp command (which essen-tially performs system identification of the human subject), agravity compensation controller, a human amplification con-troller, and a bias torque comprise the desired actuator torque E et al. : A COMPLEX STIFFNESS HUMAN IMPEDANCE MODEL WITH CUSTOMIZABLE EXOSKELETON CONTROL 3 signal. The gravity compensation controller takes the measure-ment of θ e to calculate and compensate the gravity torque τ g acting on the exoskeleton system. The human amplificationcontroller takes the measurement of τ c and multiplies it by aterm equal to negative α − , where α ≥ is an amplificationfactor. With the assistance of actuator torques produced fromthe amplification command, the human contact forces with theexoskeleton are amplified by the factor α . B. Models
We use three models describing human-exoskeleton inter-actions in our statistical tests. As preliminaries, we first define K h as the (real-valued) human elbow-joint apparent stiffness, H h as the human elbow-joint hysteretic damping, B h as thehuman elbow-joint viscous damping, M h as the moment ofinertia of the human, and M e as the moment of inertia of theexoskeleton. With the amplification control specified by thefactor α , the subject feels an attenuated inertia M e /α from theinteraction with the exoskeleton. Therefore, we also define theperceived inertia M h - e / α (cid:44) M h + M e /α at the elbow joint.The first model is a passive linear model with viscousdamping and stiffness: S h - e / α ( s ) = M h - e / α s + B h s + K h . (M1)Replacing the viscous damping in (M1) by a hysteretic damp-ing we arrive at our second model: S h - e / α ( s ) = M h - e / α s + H h j + K h , (M2)where a complex stiffness appears. Finally, to generalize (M1)and (M2), we consider a third model with both viscous andhysteretic damping: S h - e / α ( s ) = M h - e / α s + B h s + H h j + K h . (M3)In order to take advantage of the clean human cuff sensorsignal, we express these models in terms of the dynamicstiffness of the human alone, S h ( s ) = τ c ( s ) /θ e ( s ) , using thefollowing three equalities to learn the model parameters of(M1)–(M3) respectively: S h ( s ) = M h s + B h s + K h , (1) S h ( s ) = M h s + H h j + K h , (2) S h ( s ) = M h s + B h s + H h j + K h . (3)The original transfer function can be recovered by adding inthe exoskeleton inertia term S h - e / α ( s ) = S h ( s ) + α M e s . Byre-casting the parameter estimation problem as the problem ofestimating this re-creation of S h - e / α ( s ) , we can take advantageof the clean sensor data and avoid various corrupting effects inthe τ s ( s ) signal. Since the actual exoskeleton’s dynamics arebypassed, the potential influence of unmodeled exoskeletondamping on the estimated parameters is eliminated. C. Experimental Protocol for Modeling
The modeling study consists of nine perturbation experi-ments with healthy subjects between the ages of - ,where subjects A-E are females and subjects F-J are males. TABLE IM
ODELING E XPERIMENT P ARAMETERS
Exp α Load Grip Bias Amplitude Frequency Range(kg) (kg) (Nm) (Nm) (rad / s to rad / s)1 1 4 . × to 2 × . . × to 3 × . . × to 4 × . The nine experiments are separated into three groups ofthree experiments. The three experiments in each group areconducted with a . load and an α value of (corre-sponding to no amplification), , and . The gravity torque ofthe load and the exoskeleton itself are cancelled out by thegravity compensation feature of the controller, while the totalinertia is attenuated by a factor of α due to the amplificationfeature.Each of the three experimental groups are differently per-turbed to achieve variation in elbow stiffnesses. Because thestiffness is determined by both muscle co-contraction andcontraction to resist an external torque, we induce variation instiffness by having each subject squeeze an adjustable forcehand grip and by applying a bias torque from the actuator. Thethree experimental groups are divided into pairs of grippingforces and bias torques. The first group uses a
10 kg grippingforce and a bias torque. The second group uses
14 kg and . And the third group uses
27 kg and .Each of the nine perturbation experiments includes ten - sec periods. The bias torque is gradually added during thefirst of each period while the subject raises the forearmto around a ° angle from the resting position and starts tosqueeze the hand grip. Then, a sinusoidal perturbation signalis added for the next
10 sec . After the sinusoidal perturbationsignal finishes, the bias torque is gradually subtracted foranother with the subject bringing the arm back to theresting position and relaxing the hand. To avoid fatigue, thesubject rests for the next
40 sec in each period.In order to capture the natural frequency of the human elbowjoint wearing the exoskeleton, we set different values of theperturbation frequency for the different groups of experimentspreviously described. The three experimental groups use / s , / s , and / s for perturbation in the firsttime period. For other time periods, we set the perturbationfrequencies to be . times the frequency of the previousperturbation.The amplitude of the sinusoidal perturbation signal is set to . However, after the perturbation frequency is higherthan the natural frequency of S h - e / α ( s ) , the inertia effect M h - e / α starts dominating the dynamic response and thereforethe angle of displacement θ e becomes less and less sensitiveto the torque excitation. Thus, starting at the perturbationtime period for each experiment, we increase the amplitude ofthe perturbation signal by . times in order to increase the IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. XX, NO. XX, XXXX 20XX P h s . ( d e g ) M ag . ( d B ) P h s . ( d e g ) M ag . ( d B ) M ag . ( d B ) P h s . ( d e g ) (b) Subject B : S h-e/ α (s) v . s . Exp . h-e/ α (s) v . s . Exp . h-e/ α (s) v . s . Exp . h-e/ α (s) v . s . Exp . h-e/ α (s) v . s . Exp . h-e/ α (s) v . s . Exp . h (s) v . s . Exp . h (s) v . s . Exp . h (s) v . s . Exp . ω (rad / s) ω (rad / s) ω (rad / s)Exp . . . . . . . BSub . F CuffSub . BSub . F CuffSub . BSub . FExp . . . . . . . . . . . . Fig. 3. Bode plots for S h - e / α ( s ) showing all experiments for subject B in (a)-(c) and for subject F in (d)-(f). Bode plots for S h ( s ) in (g)-(i) showingthe mean and standard error across each experimental group for subjects B and F. The gray dots in (g)-(i) show the dynamic stiffness of the cuffidentified through a superposition test. The dash lines show the fitted curves using M3 and (3) . sensitivity to the torque excitation.In the end, we identify the three models of S h - e / α for all subject experiments using linear regression in the frequencydomain obtained from time domain data. The parameters ofall nine experimental settings are summarized in Tab. I. D. Statistical Analysis
Since we split each experiment into ten - sec periodswith - sec resting time within each period, the responseof S h - e / α ( s ) to the sinusoidal perturbation in each periodcompletely dies out before the next period. Therefore, for thepurposes of statistical testing, we can safely assume statisticalindependence between any two single-frequency data pointsin each experiment.Regarding the - sec sinusoidal perturbation within eachperiod, only the data from the second - sec part of theperturbation is used for calculating each frequency domainsample. Because the first - sec perturbation time is greaterthan the settling time for all S h - e / α ( s ) identified in ourexperiments, the output response reaches sinusoidal steady-state before entering the second - sec perturbation time period.For each experiment on each subject, we calculate theresidual square sum ( RSS ) for all three models, denotedas R sub-expM1 , R sub-expM2 and R sub-expM3 respectively, where sub = A , B , · · · , J and exp = 1 , , · · · , are the indices of subjects and experiments. For i = 1 , , , let us define R subMi (cid:44) (cid:88) exp=1 R sub-expMi , R expMi (cid:44) J (cid:88) sub=A R sub-expMi , (4) R allMi (cid:44) (cid:88) exp=1 J (cid:88) sub=A R sub-expMi . (5)In order to compare the significance of B h s and H h j in thehuman-exoskeleton interaction model, we conduct F-tests foreach of the two three-parameter models (M1 and M2) againstthe generalizing four-parameter model (M3). Our F-statisticaccounts for frequency domain data. For i = 1 , , F subMi-M3 = R subMi − R subM3 R subM3 · (2n − · n exp (4 − · n exp , (6) F expMi-M3 = R expMi − R expM3 R expM3 · (2n − · n sub (4 − · n sub , (7) F allMi-M3 = R allMi − R allM3 R allM3 · (2n − · n exp · n sub (4 − · n exp · n sub , (8)where n sub = 10 is the number of subjects, n exp = 9 is thenumber of experiments per subject, n = 10 is the numberof complex value samples in the frequency domain, and thefactor of two represents statistical independence between thereal and imaginary parts of each sample. E et al. : A COMPLEX STIFFNESS HUMAN IMPEDANCE MODEL WITH CUSTOMIZABLE EXOSKELETON CONTROL 5 III. M
ODELING R ESULTS
A. Frequency Domain Results
The frequency data for the two most representative subjectsare shown in Fig. 3. Fig. 3.(a)-(f) shows the Bode plots of S h - e / α ( s ) . Similarly to [24, Fig. 4], [25, Fig. 3], [26, Fig. 2],the phase for each experiment shows a non-zero value (near ° for subjects B and F) at low frequencies. This type ofphase shift is very different from the phase shift values usuallydescribed by linear systems with viscous damping where thephase shift approaches zero as ω → .Since S h ( s ) is unaffected by changes of M h - e / α , we com-pute the statistics for all three experiments in each experimen-tal group. Fig. 3.(g)-(i) shows the mean and standard error foreach experimental group. Similarly to [27, Fig. 6], the phaseshift in each experimental group changes very little across awide range of frequencies before it reaches the second orderzero at the natural frequency ω h of S h ( s ) .Tab. II shows the mean and standard error for the phaseshift of S h ( s ) in each experimental group across differentfrequencies. The data for the last three frequencies is excludedfrom the calculation due to the effect of the second order zeroat ω h . B. Model Comparison Results
We now focus on the statistical significance analysis pre-sented in Fig. 4. Fig. 4.(a) shows a subject-wise comparisonof the significance of the terms B h s and H h j that we use inM3. A critical F-statistic value of . is calculated for . false-rejection probability with (9 , degrees of freedom.The results show that the values of F subM1-M3 for all subjectsare higher than the critical F-statistic value. In particular, thevalues of F subM1-M3 for subjects B, D, and F exceed . Theseresults prove that the existence of H h j in M3 significantlyimproves modeling accuracy of S h - e / α for all subjects. Thevalues of F subM2-M3 are mostly below the critical F-statistic valueexcept for subjects A and C. Another observation is that thevalue of F subM2-M3 is lower than the value of F subM1-M3 for mostof the subjects except for subject A.Fig. 4.(b) shows an experiment-wise comparison of themodels. A critical F-statistic value of . is calculatedfor . false-rejection probability with (10 , degreesof freedom. The results show that the values of F expM1-M3 for all subjects are much higher than the critical F-statisticvalue. These results prove that the existence of H h j in M3significantly improves modeling accuracy of S h - e / α for allstiffness and inertia settings. The values of F expM2-M3 are mostlylower than the critical F-statistic value except for the threeexperiments with amplification factor α = 4 (Exp. , , ).Also, we can see a clear increment of the values of F expM2-M3 (i.e. the significance of B h s in M3) as α gets higher.Regarding the significance of the terms B h s and H h j usedin M3 over all subjects and all experiments, a critical F-statisticvalue of . is calculated for a . false-rejection probabilitywith (90 , degrees of freedom. The value of F allM1-M3 ismuch larger than . while the value of F allM2-M3 is onlyslightly above . . Although the effect of B h s cannot becompletely ignored based on the results of these F-tests, we TABLE IIO
BSERVED P HASE S HIFTS
Subject Exp . . . . E . Mean S . E . Mean S . E . A 27 . . . . . .
5B 27 . . . . . .
0C 16 . . . . . .
7D 34 . . . . . .
0E 17 . . . . . .
2F 33 . . . . . .
0G 23 . . . . . .
6H 19 . . . . . .
2I 18 . . . . . .
4J 19 . . . . . . . . . . . .
2A B C D E F G H I J10 . . . . . . . . . . . . . . . . . . . . (a) F-test v . s . Subject
FM1-M3FM2-M3 . . . . . . . . . . . . . . . . . . . . (b) F-test v . s . Experiment
FM1-M3FM2-M3 -2 -2 Fig. 4. Bar charts on log scale show first, all F subMi-M3 in (a), secondall F expMi-M3 in the first nine columns of (b), and third F allMi-M3 in thelast column of (b), for i = 1 , . The solid line appears on a bar if theF-statistic value is over the critical F-statistic value with a false-rejectionprobability of 0.05. can claim that the term H h j has much more significance thanthe term B h s as used in M3. C. Complex Stiffness Results
In our single-subject pilot study [35], a linear regressionis applied to describe the relationship between H h and K h identified from M2 and M3. However, based on our frequencydomain results, the phase shifts of S h - e / α ( s ) for some subjectsare not consistent over different stiffness values. Therefore, alinear relationship between H h and K h is not always ensuredfor all subjects. Instead, we apply linear regression betweenthe base logarithms of H h and K h and use it to identifya power law between these two parameters. Since the valueof H h is not guaranteed to be positive from the parameteridentification of M3, we only calculate the power law between H h and K h of M2.Tab. III shows the identified parameter values of H h and K h using M2, with a coefficient of determination ( R ) in therange of . ∼ . . We define β and β as the interceptand slope of the linear regression equation between the base logarithms of H h and K h . From the parameter identification IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. XX, NO. XX, XXXX 20XX
TABLE IIIC
OMPLEX S TIFFNESS P ARAMETERS FOR M2 Exp Parameter SubjectA B C D E F G H I J Average1 Kh (Nm / rad) 12 .
68 28 .
67 17 .
76 16 .
88 11 .
55 13 .
41 17 .
95 17 .
37 18 .
85 13 .
77 16 . Hh (Nm / rad) 5 .
26 14 .
09 5 .
22 12 .
38 3 .
92 7 .
26 4 .
48 5 .
22 5 .
30 2 .
56 5 . R .
99 0 .
98 0 .
97 0 .
98 0 .
96 1 .
00 0 .
97 0 .
95 0 .
95 0 .
96 - Kh (Nm / rad) 16 .
05 21 .
43 12 .
62 16 .
85 14 .
72 9 .
49 19 .
39 14 .
62 20 .
16 10 .
78 15 . Hh (Nm / rad) 8 .
08 11 .
57 4 .
79 9 .
58 4 .
59 5 .
98 8 .
98 5 .
42 6 .
55 4 .
43 6 . R .
95 0 .
99 0 .
96 0 .
98 0 .
97 0 .
98 0 .
97 0 .
96 0 .
97 0 .
96 - Kh (Nm / rad) 10 .
16 18 .
59 10 .
88 10 .
60 12 .
87 8 .
83 12 .
33 12 .
41 22 .
70 10 .
03 12 . Hh (Nm / rad) 6 .
60 9 .
70 4 .
63 7 .
19 5 .
72 7 .
04 7 .
05 5 .
14 9 .
14 4 .
35 6 . R .
97 0 .
98 0 .
88 0 .
90 0 .
90 0 .
99 0 .
96 0 .
95 0 .
96 0 .
97 - Kh (Nm / rad) 28 .
69 45 .
01 30 .
81 39 .
08 35 .
16 31 .
57 41 .
56 34 .
24 62 .
06 27 .
66 36 . Hh (Nm / rad) 9 .
95 29 .
42 14 .
01 30 .
98 7 .
13 18 .
18 13 .
72 7 .
97 15 .
69 8 .
57 13 . R .
96 0 .
96 0 .
94 0 .
98 0 .
97 0 .
98 0 .
91 0 .
96 0 .
96 0 .
96 - Kh (Nm / rad) 26 .
97 32 .
64 18 .
81 27 .
25 36 .
50 23 .
94 41 .
65 47 .
88 36 .
65 17 .
92 29 . Hh (Nm / rad) 15 .
25 21 .
69 7 .
51 21 .
78 9 .
74 16 .
28 13 .
54 23 .
64 10 .
60 7 .
26 13 . R .
93 0 .
96 0 .
95 0 .
96 0 .
93 0 .
99 0 .
89 0 .
94 0 .
93 0 .
97 - Kh (Nm / rad) 24 .
23 32 .
94 25 .
85 29 .
05 26 .
09 20 .
93 24 .
41 24 .
21 26 .
29 14 .
55 24 . Hh (Nm / rad) 14 .
23 20 .
52 11 .
31 18 .
23 8 .
48 14 .
48 9 .
93 10 .
15 12 .
11 7 .
01 12 . R .
90 0 .
97 0 .
97 0 .
98 0 .
95 0 .
99 0 .
92 0 .
97 0 .
88 0 .
95 - Kh (Nm / rad) 45 .
12 73 .
23 66 .
65 54 .
75 63 .
99 63 .
31 78 .
08 55 .
26 108 .
33 60 .
11 65 . Hh (Nm / rad) 13 .
52 49 .
62 34 .
33 45 .
87 11 .
70 32 .
12 20 .
20 25 .
01 23 .
61 12 .
83 23 . R .
94 0 .
98 0 .
96 0 .
94 0 .
92 0 .
97 0 .
97 0 .
97 0 .
98 0 .
95 - Kh (Nm / rad) 52 .
45 55 .
48 63 .
19 59 .
58 59 .
37 41 .
71 67 .
47 69 .
80 68 .
61 41 .
52 57 . Hh (Nm / rad) 17 .
52 43 .
15 31 .
05 39 .
03 16 .
03 21 .
37 19 .
20 28 .
97 14 .
61 15 .
99 22 . R .
91 0 .
98 0 .
97 0 .
98 0 .
95 0 .
99 0 .
97 0 .
97 0 .
94 0 .
97 - Kh (Nm / rad) 40 .
87 65 .
45 45 .
33 46 .
37 46 .
67 39 .
03 49 .
83 63 .
44 67 .
79 33 .
68 48 . Hh (Nm / rad) 20 .
08 38 .
08 17 .
76 24 .
44 14 .
52 22 .
11 17 .
62 33 .
32 18 .
52 13 .
92 20 . R .
93 0 .
98 0 .
96 0 .
98 0 .
92 0 .
95 0 .
92 0 .
97 0 .
95 0 .
93 - β . − . − . − . − . − .
01 0 . − . − . − . − . Power Law β .
73 1 .
21 1 .
12 1 .
03 0 .
70 0 .
85 0 .
70 1 .
10 0 .
73 0 .
84 0 . R .
77 0 .
97 0 .
97 0 .
94 0 .
86 0 .
96 0 .
83 0 .
91 0 .
87 0 .
73 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . Kh (Nm / rad) Kh (Nm / rad) Kh (Nm / rad) H h ( N m / r a d ) H h ( N m / r a d ) H h ( N m / r a d ) (a) Subject B (b) Subject F (c) Subject Average . . . . . . . . . Fig. 5. Here we show plots of H h versus K h for M2 on a log scale, with the results for subjects B and F in (a)-(b) and the geometric averageacross all subjects in (c). The dash lines and the ellipsoids show the linear regression results and the co-variances on the log scale. results using M2, a very strong linear relationship betweenlogarithms is observable across all subjects, with an R valuein the range of . ∼ . . Fig. 5.(a)-(b) show the regressionresults of subjects B and F.The last three rows of Tab. III show the identified powerlaw parameters. The damping ratio and low-frequency phaseshift of M2 can be expressed as ζ h - e / α = c h / , φ h - e / α = tan − ( c h ) , (9)where c h is a hysteretic damping loss factor [34] expressed as c h (cid:44) H h /K h = 10 β · K β − h (10)obtained by substituting the power law H h = 10 β · K β h .The last column of Tab. III shows the geometric average (i.e.arithmetic average of the logarithms) of the complex stiffnessparameters across all subjects. We apply a linear regression tothe logarithms of these average values and identify a power law of β = − . , β = 0 . , and R = 0 . (Fig. 5.(c)).As the subject average stiffness increases from . to . / rad , the value of ζ h - e / α decreases from . to . ascalculated using (9) and is within a 1-standard deviation rangeof the damping ratio of the minimally damped mode of thehuman arm ( . ± . ) described in [32].As in [35], the correlation between H h and K h can be intro-duced into M2 to reduce it to a 1-parameter complex stiffnessmodel. Adopting this reduced model allows simplifying (2) to S h ( s ) = τ c /θ e = M h s + K h (1 + c h j ) , (11)and the dynamic stiffness of the human coupled with theexoskeleton S h - e ( s ) becomes, S h - e ( s ) = τ s /θ e = M h - e s + K h (1 + c h j ) , (12)where M h - e = M h + M e is the combined inertia between thehuman and the exoskeleton. Similarly to (9), the damping ratio E et al. : A COMPLEX STIFFNESS HUMAN IMPEDANCE MODEL WITH CUSTOMIZABLE EXOSKELETON CONTROL 7 and low-frequency phase shift of S h ( s ) and S h - e ( s ) can alsobe expressed as c h / and tan − ( c h ) . IV. L
OOP S HAPING M ETHODS
The amplification feedback we discuss in this section is thesame as the direct amplification feedback shown in Fig. 2 inwhich the amplification command is − τ c multiplied by α − .But instead of a constant value of α across all frequencies,we introduce a frequency dependent amplification transferfunction α ( s ) = k p · F ( s ) + 1 , where k p is a proportional gainand F ( s ) is a fractional order controller customized accordingto the complex stiffness behavior displayed by users. A. Proportional Amplification
Based on (11) and (12), the plant transfer function P ( s ) from τ d to τ c can be expressed as P ( s ) = S h ( s ) S h - e ( s ) · G SEA ( s ) = M h s + K h (1 + c h j ) M h - e s + K h (1 + c h j ) · G SEA ( s ) (13)where the SEA transfer function G SEA ( s ) = τ s /τ d acts asa 2nd order low-pass filter. Because of the high bandwidthof the SEA force controller, the natural frequency ω SEA of G SEA ( s ) is much greater than the natural frequencies ω h - e = (cid:112) K h /M h - e and ω h = (cid:112) K h /M h of S h - e ( s ) and S h ( s ) .Considering the frequency domain properties from low tohigh frequencies, P ( s ) has a pair of conjugate poles at ω h - e ,then a pair of conjugate zeros at ω h and then another pairof conjugate poles at ω SEA (Fig. 6). Between ω h - e and ω h , S h - e ( s ) is dominated by its inertia effect and the magnitudeof P ( s ) decreases while the phase decreases from °. On theother hand, S h ( s ) is still dominated by the complex stiffnessand prevents the phase moving below tan − ( c h ) − °.If we apply a very large value of k p , the gain crossoverof P ( s ) falls beyond ω SEA . The phase margin with suchcrossover is very close to zero because of the orderSEA dynamics. Also, the closed loop behavior amplifies thehigh frequency sensor noise from the actual signal from τ c (which is usually de-noised by a low-pass filter beyondthe frequency of ω SEA that makes the closed loop evenmore likely to be unstable). Similarly to [15], the crossoverfrequency cannot be placed between ω h and ω SEA becausemultiple crossovers could easily occur. Besides the multiplecrossovers, this frequency range is also outside of the testedfrequency ranges of Exp. - . The unmodeled dynamics fromthe human and cuff will cause additional stability issues ifa crossover is placed there. Instead, a new crossover can besafely placed at the frequency between ω h - e and ω h using asmaller k p (Fig. 6). As a rule of thumb, k p can be set as k p = ( ω gc /ω h ) = ( M h - e /M h ) ,ω gc (cid:44) (cid:113) K h / ( M h - e · M h ) , (14)where the crossover ω gc of k p · P ( s ) is exactly in the middlebetween ω h - e to ω h in the log scale. ω h - e ω h ω SEA · l og ( k p ) Potential Bandwidth Multi-Crossover NoiseMagnitude (dB) ω h - e ω h ω SEA − − Potential Bandwidth Multi-Crossover NoisePhase (deg) k p · P ( s ) F ( s )k p · P ( s ) · F ( s ) Fig. 6. These conceptual bode plots show P ( s ) with its poles(crosses) and zeros (circles). The various regions are color-coded: themodel is trustworthy in the green region, the blue region reflects themulti-crossover behavior which makes an amplification controller designunreliable, and the yellow region is dominated by sensor noise from τ c .A fractional-order controller F ( s ) complements a proportional controller k p by boosting the low-frequency amplification. B. Fractional Order Amplification
In [15], an additional integral term is added to the propor-tional gain k p to boost the amplification at low frequencieswhile maintaining the same crossover frequency between ω h - e and ω h . However, a PI controller has a − ° phase at lowfrequency, which can result in loss of stability if the zero ofthe PI controller is too close to the crossover.In [35], a fractional order controller was proposed to takeadvantage of the complex stiffness model, F ( s ) = k f · s − f , (15)where f is the fractional order (that is, a non-integer powerof s ) of F ( s ) and k f is a gain which allows tuning themagnitude of F ( s ) in the frequency domain. The fractionalorder controller in (15) has its magnitude decreasing − · f dB per decade and its phase staying at − · f degrees at allfrequencies.By multiplying (15) by the proportional gain k p , the gainof the controller is increased at low frequency and reducedat high frequency (for further de-noising the measurement of τ c ). If k f is tuned to make F ( s ) have the exact same crossoverfrequency as k p · P ( s ) , we will obtain the magnitude bode plot k p · P ( s ) · F ( s ) rotated from k p · P ( s ) with pivot at the point ofthe gain crossover frequency (Fig. 6). Since the exact crossoverfrequency of k p · P ( s ) varies with the value of K h , k f can beset as k f = ˆ ω fgc , ˆ ω gc (cid:44) (cid:113) ˆ K h / ( M h - e · M h ) , (16)where ˆ ω gc is chosen as a nominal crossover frequency of k p · P ( s ) with ˆ K h (cid:44) ( ¯ K h · ¯ K h ) / being the geometric meanbetween the lower bound ¯ K h and the upper bound ¯ K h .Because of the non zero phase shift associated with thecomplex stiffness behavior, a positive phase margin can beguaranteed if < f < tan − ( c h ) / . The fractional order IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. XX, NO. XX, XXXX 20XX ˜ τ c θ e ˜ τ c ˜ τ s Subject B : f = 0 . − τ c ( N m ) ˜ τ c ( N m ) − time (s) − τ c ( N m ) Subject B : f = 0 . − . . τ c ( N m ) Subject G : f = 0 . θ e (r a d ) θ e (r a d ) . . . . τ s ( N m ) − Subject G : f = 0 . ˜ τ s ( N m ) − Subject D : f = 0 . Subject D : f = 0 . / sBias : 4 NmGrip : 14 kgMotion : 10 rad / s Fig. 7. (a)-(c) show the responses of θ e (dash red) and ˜ τ c (solid blue) for the post-tuning tests on subjects B, D and G. (d)-(f) show the responsesof ˜ τ s (dash red) and ˜ τ c (solid blue) for the amplification tests on subjects B, D and G. controller can be precisely designed for all subjects based onthe values of β and β shown in Tab. III through the settings f = (cid:40) tan − (10 β · ¯ K β − h ) / − φ/ , if β < , tan − (10 β · ¯ K β − h ) / − φ/ , if β ≥ , (17)where φ > is a user-defined guaranteed phase margin.Differently from [35] where a constant c h is assumed, (17)considers the lowest value of c h of a subject in the stiffnessrange [ ¯ K h , ¯ K h ] .As a fractional-order controller, F ( s ) cannot be imple-mented directly into a computational control process. How-ever, from [35], we can approximate it as the product of many1st order lag filters, F ( s ) = k f p f1 · n (cid:89) i =1 s/ z i s/p i , (18) z i /p i = r zp , for i = 1 , , · · · , n (19) p i /p i − = r pp , for i = 2 , , · · · , n, (20)where n is the number of lag filters and the pole and the zerofor each lag filter are − p i and − z i . We define r zp such thatall lag filters have an equal distance between the pole and thezero, and we define r pp such that there is a constant distancebetween adjacent lag filters (in log frequency space). Theamplification controller in (18) functions as a fractional-orderfilter in the frequency range of [ p , z n ] rad / s . The fractionalorder can be approximated as f ≈ log( r zp ) / log( r pp ) . C. Experimental Protocol for Loop Shaping
Based on (17), we conducted loop shaping experiments onsubjects B, D and G who, respectively, had the highest value,the closest value to , and the lowest value of β acrossall subjects. Our loop shaping study consists of two tuningexperiments and two amplification experiments.The value of M e is .
01 kg · m , which includes a . load at the end of the exoskeleton arm. Although we donot measure the value of M h directly from our subjects, anaverage M h of .
11 kg · m can be obtained from a -subject measurement study presented in [27]. Based on these inertiavalues and (14), we set k p = 3 . .As shown in Tab. III, the human stiffness changes from . to .
33 Nm / rad across all subjects and all ex-periments, which gives us a nominal value of ˆ K h =32 .
96 Nm / rad . Based on (16) and (18), we compute ˆ ω gc ≈
10 rad / s and implement an approximate fractional-order con-troller using lag filters with p = 1 and r pp = 10 . suchthat ˆ ω gc is located at the center of the frequency range definedby [ p , p ] rad / s .The two tuning experiments we perform aim to find out thefractional order of a subject where the minimum phase margin φ is near zero. From (17), we gradually increase the fractionalorder, f , from low value to higher values until the exoskeletonstarts to oscillate. We do that with subjects employing lowand high human stiffness behaviors. The maximum stablevalue of f will be the lower value between the two stiffnesscases. An important advantage is that this tuning strategy doesnot require prior knowledge of the human complex stiffness.Similarly to the modeling experiments previously presented,we regulate the low and high stiffness of a subject by settingthe gripping force as and
27 kg and the bias torque as and .After the tuning experiments outlined above, we subtract . from the marginally stable fractional order, which pro-vides a minimum phase margin φ = 10 . °. Then, we conducttwo amplification experiments both with a gripping force of
14 kg and a bias torque of . These two experiments areconducted using sinusoidal voluntary movements performedby the subjects with frequencies of and
10 rad / s . Thevoluntary sinusoidal movements are guided by showing thesubject a visual signal of the actual joint position θ e and thedesired sinusoidal wave on a screen. The amplification factor α for these sinusoidal voluntary movements can be calculatedfrom the experimental data after the experiments. V. L
OOP S HAPING R ESULTS
In order to study the performance of the tuned amplificationcontrollers for various subjects, we define ˜ τ s (cid:44) τ s + τ g − bias E et al. : A COMPLEX STIFFNESS HUMAN IMPEDANCE MODEL WITH CUSTOMIZABLE EXOSKELETON CONTROL 9 and ˜ τ c (cid:44) τ c − bias . The real-time amplification factor for theproposed amplification controller can be expressed as α ( t ) =˜ τ s ( t ) / ˜ τ c ( t ) + 1 . A. Tuning Results
After gradually increasing the fractional order for lowand high human stiffness behaviors until the exoskeletonstarts to oscillate continuously, we obtain the values f =0 . , . , . for subjects B, D, G. In order to displaythe tuning results concisely, we conduct two post-tuning testsinvolving low and high human stiffness setups. We attach aset of mechanical springs to the tip of the exoskeleton arm(Fig. 1.(b)) and quickly detach it to test the dynamic responseof the controller. The response of θ e and ˜ τ c for the post-tuningtests are shown in Fig. 7.(a)-(c).In Tab. III, we had identified that β for subject B wasgreater than . This explains why the post-tuning test forsubject B applying high stiffness is less oscillatory than thepost-tuning test applying low stiffness. Because the value of β is very close to for subject D, the results for both post-tuningtests are very similar. Similarly, the high stiffness post-tuningtest for subject G is more oscillatory than the low stiffnesspost-tuning test because β < . B. Amplification Results
When we subtract . from f for subjects B, D and G,we get the values . , . and . . The behaviors of ˜ τ s and ˜ τ c for sinusoidal voluntary movements between and / s demonstrate that the exoskeleton is stable for all subjectsusing our proposed custom robust amplification controllers(Fig. 7.(d)-(f)). The values for the gain and the phase shift for ˜ τ s and ˜ τ c during the amplification tests are shown in Tab. IV.Notice that the subjects are able to maintain values between | ˜ τ s ˜ τ c | = 2 . ∼ . with a voluntary motion of
10 rad / s .In our prior research which did not incorporate the proposedcomplex stiffness model [15], the value of | ˜ τ s ˜ τ c | was between . ∼ . at . / s (experimentally validated), and avalue of . at
10 rad / s (theoretically estimated). Therefore,our proposed control strategy shows a ∼ improvementin the magnitude when using a dynamical amplification factor α ( s ) = ˜ τ s ( s ) / ˜ τ c ( s ) + 1 . VI. D
ISCUSSION
The sensor configuration for this experiment measures thedeflection at the exoskeleton’s hinge joint as well as the humantorque using the cuff’s six-axis force/torque sensor. Using thissetup, we conduct a test of superposition for differentiatingbetween the human elbow and cuff impedances. Fig. 3.(g)-(i) and Tab. II show the frequency data and phase shiftvalues of the cuff attached to a rigid object. The phase ofthe cuff is lower than the phases of all human subjects.Furthermore, the magnitude of the cuff stiffness is above ( / rad ), which is significantly higher than thestiffness values of all human subjects as shown in Tab. II. Weconclude that the human impedance becomes the dominantfactor measured in our experiments. TABLE IVO
BSERVED D YNAMICAL A MPLIFICATION
Subject f | ˜ τ s ˜ τ c | (1 rad / s) ∠ ˜ τ s ˜ τ c (1 rad / s) | ˜ τ s ˜ τ c | (10 rad / s) ∠ ˜ τ s ˜ τ c (10 rad / s) B 0 .
60 10 . − . ° . − . ° D 0 .
44 6 . − . ° . − . ° G 0 .
10 3 . − . ° . − . ° While the noticeable phase shift values observed in thisstudy are consistent with the ankle and elbow joint phasevalues reported in [24], [25], [26], [27], some research studiesalso report very small phase shift values for other humanjoints. Yet, our proposed complex stiffness model still holdsfor those results. For example, the damping ratio for the humanknee joint reported in [37, Fig. 5] is . , which results in aphase shift of . ° using equation (9). This kind of smallhysteretic damping characteristic can be easily overlooked.Low-frequency phase shifts are found in muscle spindles[38] and arteries [39] of mammals, suggesting that joint hys-teretic damping could be due to the bio-mechanical propertiesof the human tissue. Therefore, we suspect that the humanneuromuscular system, either through muscle and tendon hys-teresis or through neural hysteretic behavior, is the mechanismbehind our hysteretic damping hypothesis.Because Fig. 3.(g)-(i) shows a consistent phase shift acrossa wide range of low frequencies, it is natural to considerthat the phase behavior of S h ( s ) has already reached a low-frequency asymptote at the lowest tested frequency and itwill not change much at lower frequencies than that. Thisis difficult to experimentally verify because lower frequenciesrequire longer experimental times making it harder for thesubjects. Nonetheless, our lowest tested frequency, / s ( ≈ . ), is below the frequencies reported in references[24], [25], [26], [27], [28]. In addition, our tested frequencyrange covers the frequencies that are important for practicalcontrol system design.In this research, we use sinusoidal perturbations to identifya frequency domain model of human complex stiffness. Ref.[24], [25] use white noise perturbations for system identifica-tion, also revealing a non-zero phase shift at low frequenciesand consistent damping ratios. Because our complex stiffnessmodel is non-causal and nonlinear, it does not have an exactmodel representation in the time domain [40]. For this reason,it is difficult to identify the human complex stiffness behaviorusing impulse, step and ramp perturbations.Although the fractional order part of the proposed ampli-fication controller is only useful for shaping the behaviorin a certain frequency range, it is sufficient to demonstratecrossover at frequencies that could not be robustly stablewith the conventional human joint model. Between this paperand [15], we have a natural comparison between the controldesign problem with and without the hysteretic adjustment tothe human impedance. The result is clear: without complexstiffness, controllers must be designed to cross over beforethe lowest natural frequency resulting from the human andexoskeleton inertia and the softest human stiffness; with themodification, the crossover can exceed this frequency byimplementing an approximate fractional-order controller. VII. C
ONCLUSION
Exoskeletons with feedback human forces must be coupledstable with the natural human impedance to avoid undesiredvibrations. This paper presents a model for human impedanceusing an imaginary stiffness term to fill an energy-dissipationrole similar to damping. The paper also presents experimentswhich demonstrate that this new term is a more significantcontributor to model accuracy than a linear damping term forcyclic motion of the elbow in the -subject cohort we studied.The loop shaping experiments demonstrate the stability andbandwidth of our controller and highlight the importance oftesting both maximum and minimum human stiffness caseswhen tuning the fractional order controller. R EFERENCES [1] N. Yagn, “Apparatus for facilitating walking, running, and jumping,”
USpatent , vol. 420179, 1890.[2] J. B. Makinson, D. P. Bodine, and B. R. Fick, “Machine augmentation ofhuman strength and endurance Hardiman I prototype project,” SpecialtyMaterials Handling Products Operation, General Electric Company,Tech. Rep., 1969.[3] H. Kazerooni and J. Guo, “Human extenders,”
Journal of DynamicSystems, Measurement, and Control , vol. 115, no. 2B, pp. 281–290,1993.[4] H. Kazerooni, “Exoskeletons for human power augmentation,” in
In-telligent Robots and Systems (IROS), 2005 IEEE/RSJ InternationalConference on . IEEE, 2005, pp. 3459–3464.[5] A. M. Dollar and H. Herr, “Lower extremity exoskeletons and active or-thoses: challenges and state-of-the-art,”
IEEE Transactions on Robotics ,vol. 24, no. 1, pp. 144–158, 2008.[6] S. C. Jacobsen and M. X. Olivier, “Contact displacement actuatorsystem,” Sep. 30 2014, US Patent 8,849,457.[7] M. Fontana, R. Vertechy, S. Marcheschi, F. Salsedo, and M. Bergamasco,“The body extender: A full-body exoskeleton for the transport andhandling of heavy loads,”
IEEE Robotics & Automation Magazine ,vol. 21, no. 4, pp. 34–44, 2014.[8] J. Zhang, P. Fiers, K. A. Witte, R. W. Jackson, K. L. Poggensee,C. G. Atkeson, and S. H. Collins, “Human-in-the-loop optimization ofexoskeleton assistance during walking,”
Science , vol. 356, no. 6344, pp.1280–1284, 2017.[9] S. Lee, J. Kim, L. Baker, A. Long, N. Karavas, N. Menard,I. Galiana, and C. J. Walsh, “Autonomous multi-joint soft exosuit withaugmentation-power-based control parameter tuning reduces energy costof loaded walking,”
Journal of Neuroengineering and Rehabilitation ,vol. 15, no. 1, p. 66, 2018.[10] K. Kong, H. Moon, D. Jeon, and M. Tomizuka, “Control of an exoskele-ton for realization of aquatic therapy effects,”
IEEE/ASME Transactionson Mechatronics , vol. 15, no. 2, pp. 191–200, 2010.[11] B. Kim and A. D. Deshpande, “An upper-body rehabilitation exoskeletonharmony with an anatomical shoulder mechanism: Design, modeling,control, and performance evaluation,”
The International Journal ofRobotics Research , vol. 36, no. 4, pp. 414–435, 2017.[12] G. Lv and R. D. Gregg, “Underactuated potential energy shaping withcontact constraints: Application to a powered knee-ankle orthosis,”
IEEETransactions on Control Systems Technology , vol. 26, no. 1, pp. 181–193, 2018.[13] H.-D. Lee, B.-K. Lee, W.-S. Kim, J.-S. Han, K.-S. Shin, and C.-S. Han,“Human–robot cooperation control based on a dynamic model of anupper limb exoskeleton for human power amplification,”
Mechatronics ,vol. 24, no. 2, pp. 168–176, 2014.[14] S. P. Buerger and N. Hogan, “Complementary stability and loop shapingfor improved human–robot interaction,”
IEEE Transactions on Robotics ,vol. 23, no. 2, pp. 232–244, 2007.[15] B. He, G. C. Thomas, N. Paine, and L. Sentis, “Modeling and loopshaping of single-joint amplification exoskeleton with contact sensingand series elastic actuation,” in . IEEE, 2019, pp. 4580–4587.[16] H. Huang, H. F. Cappel, G. C. Thomas, B. He, and L. Sentis, “Adap-tive compliance shaping with human impedance estimation,” in . IEEE, 2020, pp. 5131–5138. [17] E. J. Rouse, L. J. Hargrove, E. J. Perreault, and T. A. Kuiken, “Esti-mation of human ankle impedance during the stance phase of walking,”
IEEE Transactions on Neural Systems and Rehabilitation Engineering ,vol. 22, no. 4, pp. 870–878, 2014.[18] H. Lee and N. Hogan, “Time-varying ankle mechanical impedanceduring human locomotion,”
IEEE Transactions on Neural Systems andRehabilitation Engineering , vol. 23, no. 5, pp. 755–764, 2015.[19] H. Lee, E. J. Rouse, and H. I. Krebs, “Summary of human anklemechanical impedance during walking,”
IEEE journal of translationalengineering in health and medicine , vol. 4, pp. 1–7, 2016.[20] A. L. Shorter and E. J. Rouse, “Mechanical impedance of the ankleduring the terminal stance phase of walking,”
IEEE Transactions onNeural Systems and Rehabilitation Engineering , vol. 26, no. 1, pp. 135–143, 2018.[21] ——, “Ankle mechanical impedance during the stance phase of running,”
IEEE Transactions on Biomedical Engineering , pp. 1–1, 2019.[22] D. Bennett, J. Hollerbach, Y. Xu, and I. Hunter, “Time-varying stiffnessof human elbow joint during cyclic voluntary movement,”
ExperimentalBrain Research , vol. 88, no. 2, pp. 433–442, 1992.[23] E. J. Rouse, R. D. Gregg, L. J. Hargrove, and J. W. Sensinger, “Thedifference between stiffness and quasi-stiffness in the context of biome-chanical modeling,”
IEEE Transactions on Biomedical Engineering ,vol. 60, no. 2, pp. 562–568, 2012.[24] G. Agarwal and C. Gottlieb, “Compliance of the human ankle joint,”
Journal of Biomechanical Engineering , vol. 99, no. 3, pp. 166–170,1977.[25] G. L. Gottlieb and G. C. Agarwal, “Dependence of human anklecompliance on joint angle,”
Journal of Biomechanics , vol. 11, no. 4,pp. 177–181, 1978.[26] G. Zahalak and S. Heyman, “A quantitative evaluation of the frequency-response characteristics of active human skeletal muscle in vivo,”
Journal of Biomechanical Engineering , vol. 101, pp. 28–37, 1979.[27] S. C. Cannon and G. I. Zahalak, “The mechanical behavior of activehuman skeletal muscle in small oscillations,”
Journal of Biomechanics ,vol. 15, no. 2, pp. 111–121, 1982.[28] I. Hunter and R. Kearney, “Dynamics of human ankle stiffness: Variationwith mean ankle torque,”
Journal of Biomechanics , vol. 15, no. 10, pp.747 – 752, 1982.[29] J. Becker and C. Mote, “Identification of a frequency response model ofjoint rotation,”
Journal of Biomechanical Engineering , vol. 112, no. 1,pp. 1–8, 1990.[30] P. Weiss, I. Hunter, and R. Kearney, “Human ankle joint stiffness overthe full range of muscle activation levels,”
Journal of Biomechanics ,vol. 21, no. 7, pp. 539–544, 1988.[31] E. Sobhani Tehrani, K. Jalaleddini, and R. E. Kearney, “Ankle jointintrinsic dynamics is more complex than a mass-spring-damper model,”
IEEE Transactions on Neural Systems and Rehabilitation Engineering ,vol. 25, no. 9, pp. 1568–1580, 2017.[32] E. J. Perreault, R. F. Kirsch, and P. E. Crago, “Multijoint dynamics andpostural stability of the human arm,”
Experimental Brain Research , vol.157, no. 4, pp. 507–517, 2004.[33] T. E. Milner and C. Cloutier, “Compensation for mechanically unstableloading in voluntary wrist movement,”
Experimental Brain Research ,vol. 94, no. 3, pp. 522–532, 1993.[34] R. E. D. Bishop and D. C. Johnson,
The Mechanics of Vibration .Cambridge University Press, 1960.[35] B. He, H. Huang, G. C. Thomas, and L. Sentis, “Complex stiffnessmodel of physical human-robot interaction: Implications for control ofperformance augmentation exoskeletons,” in . IEEE,2019, pp. 6748–6755.[36] N. Paine, S. Oh, and L. Sentis, “Design and control considerations forhigh-performance series elastic actuators,”
IEEE/ASME Transactions onMechatronics , vol. 19, no. 3, pp. 1080–1091, 2014.[37] M. Pope, R. Crowninshield, R. Miller, and R. Johnson, “The static anddynamic behavior of the human knee in vivo,”
Journal of biomechanics ,vol. 9, no. 7, pp. 449–452, 1976.[38] R. Poppele and R. Bowman, “Quantitative description of linear behaviorof mammalian muscle spindles.” journal of Neurophysiology , vol. 33,no. 1, pp. 59–72, 1970.[39] N. Westerhof and A. Noordergraaf, “Arterial viscoelasticity: a general-ized model: effect on input impedance and wave travel in the systematictree,”
Journal of Biomechanics , vol. 3, no. 3, pp. 357–379, 1970.[40] J. A. Inaudi and J. M. Kelly, “Linear hysteretic damping and the hilberttransform,”