Quantitative comparison of time-varying system identification methods to describe human joint impedance

Accurate and swift tuning of joint impedance is crucial to perform movement and interaction with our environment. Time-varying system identification enables quantification of joint impedance during movement. Many methods have been developed over the years, each with their own mathematical approach and underlying assumptions. Yet, for the identification of joint impedance, a systematic comparison revealing each method ’ s unique strengths and weaknesses, is lacking. Here, we propose a quantitative framework to compare these methods. The framework is used to review five time-varying system identification methods using both simulated data and experimental data. These methods included three time-domain methods: ensemble, short data segment


Introduction
Our neuromuscular system is challenged to provide effective and energy-efficient control of every movement we make, while interacting with a highly dynamic environment.One way we accomplish this is by appropriately setting our joint's mechanical properties.Altering our joint's mechanical properties can be achieved by modulation of the mechanical properties of the muscle and tendon via muscle contraction or by modulating properties of the reflexive pathways (Kearney & Hunter, 1990).The joint's mechanical properties are often quantified as the joint's impedance, a metric describing the dynamic relationship between an imposed joint displacement and the restoring joint torque.A joint's impedance undergoes continuous and rapid adaptations during movement.Examples include the modulation of ankle impedance during the gait cycle (Lee & Hogan, 2015) and elbow impedance during reaching movements (Popescu et al., 2003).These adaptations in joint impedance are crucial to maintain stable interaction with the environment with minimal energy consumption (Liu & Todorov, 2007).The importance of proper modulation of impedance becomes even more evident in individuals that have central or peripheral nervous system damage.Spasticity following stroke (Meskers et al., 2009;Mirbagheri et al., 2008), or freezing episodes when suffering Parkinson's Disease (Johnson et al., 1993;Prochazka et al., 1997) reveal how improper control of joint impedance during movement can arise due to malfunctioning of the human control system.
Accurately quantifying joint impedance requires the use of external perturbations.Early work without external perturbations used the concept of quasi-stiffness during movement by exploiting the derivative of the torque-angle relationship with respect to angle (Latash & Zatsiorsky, 1993).Rouse and colleagues showed that quasi-stiffness cannot provide an accurate estimate of joint stiffnessthe static component of the impedanceduring movement, as joint torques and angles are rapidly varying (Rouse et al., 2013).Hence, joint impedance cannot be determined directly from measured joint angle and torque during normal movement.
System identification techniques can provide accurate descriptions of joint impedance by use of external perturbations (Kearney & Hunter, 1990;Kearney et al., 1997;van der Helm et al., 2002).System identification has been widely used to investigate joint impedance during different postural tasks and under varying experimental conditions.In these studies, the joint's operating point is systematically varied (quasi-stationary) and joint impedance is considered time-invariant during every recording condition, enabling the use of time-invariant system identification.This approach has been used to demonstrate that joint stiffness varies with joint torque, joint angle, and muscle activation (Hunter & Kearney, 1982;Kearney & Hunter, 1982;Mirbagheri et al., 2000).However, time-invariant techniques cannot be employed when assessing joint impedance during movement.During movement, joint torque, joint angle, and muscle activation vary throughout the movement and the joint does not stay in one operating point.Extrapolating joint impedance from postural conditions to impedance during movement is inappropriate as, for example, joint impedance is lower during movement than posture (Bennett et al., 1992).Hence, joint impedance during movement can only be quantified appropriately when considering the system as time-varying, where the time-varying variable is e.g., joint angle.The challenge has been to find a time-varying system identification method that provides an accurate description of joint impedance during movement.
An abundance of time-varying system identification methods exist, in both time-and frequency domain, for the identification of joint impedance.However, it is unclear which method is best suited for different applications.Each method provides distinct opportunities, using a unique mathematical approach, to tackle the specific challenges faced while dealing with human data.Challenges include: the poor signal-to-noise (SNR) ratio of human kinematic and kinetic data, limited recording time available, e.g., due to participant fatigue, and lack of a priori knowledge regarding to joint impedance during movement.One of the first methods applied to human data was a non-parametric ensemble-based method, which required very little a priori information about the expected joint dynamics.However, the method requires hundreds of realizations of input and output data with the same underlying time-varying behavior and a strict timing between the timevarying behavior and the input signal (Lacquaniti et al., 1982;Soechting et al., 1981).The latter is difficult to achieve in human experiments and therefore mathematical solutions were sought to provide a robust estimate of joint impedance without the need for time-locking the input signal and time-varying behavior across trials (Lortie & Kearney, 2001;MacNeil et al., 1992).Reliable estimates of joint impedance without the need for a large amount of data can be obtained using various assumptions.These include the assumption that joint impedance varies little inside a short time window (Ludvig & Perreault, 2012); or assuming that the time-varying joint impedance is proportional to e.g., muscle activity or joint position, and this dependence can therefore be included a priori by basis or scheduling functions (Golkar et al., 2017;Sobhani Tehrani et al., 2013, 2014).Moreover, when using parametric time-varying methods the joint impedance is assumed to match a predefined model structure, e.g., that of a mechanical system with inertia, viscosity, and stiffness (Bennett et al., 1992;Kirsch & Kearney, 1997;Lacquaniti et al., 1993).However, such a second-order model may not always be sufficient to describe the joint impedance accurately (Sobhani Tehrani et al., 2017).Hence, the various assumptions underlying time-varying system identification methods directly affect their ability to reliably identify time-varying joint impedance.As a result, the decision on which method to adopt when is not straightforward.
In this study, we perform a systematic comparison of five timevarying methods that are used to quantify joint impedance.First, a summary is provided of the mathematical principles behind each method.Then the identification performance of each method is assessed using simulation data obtained by simulating time-varying impedance of the ankle joint.Finally, the methods are used to quantify joint impedance from experimental data recorded from the ankle joint.Together, the collected information is used to present a framework based on which researchers can decide which method to use, and how to benchmark any newly developed methods.

Algorithms
The included methods were selected such as to obtain a wide-range of time-varying methods with unique properties.A detailed description of each method can be found in the cited work.In this work, we focus on the identification of the joint impedance and do not attempt to separate the dynamics into their intrinsic and reflexive contributions as others have done under time-invariant conditions (Kearney et al., 1997;van der Helm et al., 2002).For all methods we use an angular input u(t) and torque output y(t), while the time-varying dynamics are captured by either time-varying impulse response functions (IRF) h(t, τ) or time-varying frequency response functions H(s, t).

Ensemble Impulse response function (eIRF)
The ensemble impulse response function method was one of the earliest methods applied for identification of time-varying human joint properties (Soechting et al., 1981).The method was later improved for the identification of human joint impedance by MacNeil et al. (1992) and Lortie and Kearney (2001).
Ensemble-based methods use input and output data across an ensemble of realizations to estimate impedance at each time point.This approach requires the collection of data containing many realizations of input and output data while the studied system undergoes the same time-varying behavior.The output y r (t) of a system for each realization r is represented by the convolution of a time-varying IRF h(t, τ) and input u r (t): Assuming the IRF is of finite length and therefore h(t,τ) = 0 for |τ| > τ 0 while the data are recorded digitally with a sample time Δt the output y r (t) can be approximated by a discrete convolution: where i is the time index and n = τ 0 /Δt.Altogether, this will result in a system consisting of R linear equations, R being the number of inputoutput realizations available.In matrix form this set of equations can be described as: where the length of vector Y i is determined by the number of realizations recorded (R), the length of vector H i is equal to p = 2n + 1 assuming a two-sided IRF (implying a non-causality), and U i is a matrix of size p × R. H i contains the IRF for each time point i, and can be found by solving (2) when input and output are known.Different cost functions result in different results, but in this work a solution is found by minimizing the squared error between the recorded and predicted output (MacNeil et al., 1992).This vector can be found by performing a singular value decomposition of U i to find a least-squares approximation of the system's IRF for each time point i.However, in the current study we use an alternative to solve (2), using a correlation approach in combination with a pseudoinverse (Lortie & Kearney, 2001).This solution was found to provide a more robust estimate of the system's dynamics under low SNR values.Hence, (2) can be rewritten as: where φ uu and φ yu are time dependent auto-and cross-correlation functions.In matrix form: (5) H i can be solved by using a singular value decomposition, removing those terms that present a small singular value by a Minimum Description Length cost function, and calculating the pseudoinverse (Lortie & Kearney, 2001).

Short Data segments (SDS)
Ludvig and Perreault (2012) combined both time-invariant and time-varying correlations to estimate the system's dynamics during short data segments.They noticed that the need for many realizations of input-output data with the same underlying system dynamics makes the eIRF method practically difficult to use for physiological systems, e.g., humans, where recording time is limited.Moreover, when the time-variation is slow with respect to the system's dynamics, or contains periods of time-invariant dynamics, the eIRF approach is underutilizing the available data.
Ludvig and Perreault (2012) used the same correlation based method (4) as Lortie and Kearney (2001).For time-varying conditions the cross-correlations can be calculated at each time point i and lag k using: where the same time-varying system dynamics is assumed across each realization r of the total set (R).To effectively average both across time within a short data segment as well as across realizations, the timeinvariant and time-varying estimators of cross-correlations can be combined.The multisegment correlations were defined as: where N is the number of points in each segment and t is the time at the middle of each realization.Combining (4) and (7), h(i, j) is solved by using a pseudoinverse which was found to provide a robust estimate of the system's dynamics under low SNR values (Lortie & Kearney, 2001).

Basis Impulse response function (bIRF)
The basis IRF methods proposed by Guarin and Kearney (2017), and improved in Guarin and Kearney (2018b), combines ensemble and deterministic methods to estimate the time-varying joint impedance.A linear combination of basis functions is used to estimate the time-varying IRF h(t, τ) according to: where {Λ j (t)} j=b j=0 are a set of time-varying basis functions and λ τ,j their coefficients.Subsequently, the output y(t) is the convolution of the input with the time-varying IRF as determined by the time-varying basis functions: For the estimation of the coefficients of the basis functions the algorithm uses a sparse linear identification algorithm.A Bayesian identification algorithm (relevance vector machine (Bishop, 2006)) was chosen to eliminate the basis functions that do not reduce the prediction error.

Ensemble spectral method (ESM)
The ensemble spectral method by Schouten et al. (in preparation), in line with the SDS method, uses short data segments in which dynamics are considered time-invariant.However, rather than using the cross-and auto-correlations their frequency-domain equivalents are used: crossand auto spectral density.Before transferring to the frequency domain, data are segmented around time point t, and each segment is multiplied by a Hanning window.Subsequently, data are transformed to the frequency domain using the fast Fourier Transform.The cross-(S yu (f, t)) and auto spectral densities (S uu (f, t)) are then calculated using: where N is the number of samples in each segment, U(f , t) and Y(f, t) the Fourier transformed input-and output data of the segment around time t, and U * (f, t) the complex conjugate of U(f , t).Taking the mean cross-(S yu (f, t)) and auto spectral densities (S uu (f,t)) across all realizations, the time-varying frequency response function H yu (f, t) is determined by: H yu (f, t) represents the time-varying joint impedance.

Kernel-based regression (KBR)
Another frequency-domain method is the kernel-based regression method, originally described by Lataire et al. (2017).The method has been successfully applied to the identification of joint impedance with a limited amount of data but was found to struggle capturing fast dynamics.Therefore, in our comparison we use the modified KBR method as proposed by Cavallo et al. (2020) which enables identification of faster changes in joint impedance.
For the KBR method we assume that the joint's input and output signals, u(t) and y(t), satisfy a linear differential equation of the form: where a m (t) and b m (t) are the time-varying coefficients which are smooth periodic functions of t, and f ext (t) is an additional smooth error term.The time-varying coefficients are estimated via kernel-based regression.In essence, the estimator is defined as the following minimizer: where E is the discrete Fourier transform of the equation error (i.e., the difference between the left and right hand side of ( 13)), evaluated in K int , and a m and b m are vectorised versions of a m (t) and b m (t) in t = 0,T s , …, (N − 1)T s .K int represents the bins of the frequency band of interest and σ2 E is (an estimate of) the noise variance of E. R(a m , b m ) is a quadratic regularization term, to impose the smoothness and the approximate time-periodicity of the estimates: The parameters a n and b n are decomposed into a n and b n (which are constant vectors to which no regularization is applied) and ȃn and bn (which are time-varying and regularized to impose their smoothness).
For the kernel matrix K the local periodic kernel is used: in which l 1 determines the smoothness of the estimated time-varying coefficients, γ represents the inverse of the amount of regularization applied, defining a bias versus variance trade-off of the estimated coefficients, p the length of the period and l 2 the length scale of the periodicity.

Methods
This section describes the simulation study that was used to quantify the performance of the time-varying methods when identifying simulated joint impedance.A preliminary comparison of the eIRF, SDS and bIRF algorithms based on a subset of the same simulation data and performance metrics was presented previously in a conference abstract (van de Ruit et al., 2021).

Modeling human joint impedance
An anti-causal time-varying open-loop model of human joint impedance was used with an angular input u(t), torque output y(t), and measurement noise n(t) (Fig. 1).This model represents a protocol commonly employed in human experiments where participants perform a force task, exerting a (time-varying) joint torque, while receiving angular perturbations (Kirsch & Kearney, 1997).Theoretically, this model is described by: H joint (s, t) is the time-varying human joint impedance and is modelled using a 2nd-order inertia-spring-damper system k: in which I is the joint inertia, b the joint viscosity and k(t) the timevarying joint stiffness.The Laplace variable s equals j2πf (f represents the frequency).Hence, only the intrinsic joint mechanical properties were modeled with time-varying joint stiffness.

Simulation parameters & perturbation signal
The model was implemented in MATLAB 2019b -Simulink 9.7 (The MathWorks, Inc., Natick, Massachusetts, United States).The parameters of H joint (s, t) were chosen to represent the human ankle joint impedance while being seated.Joint inertia and viscosity remained time-invariant (I = 0.02 kg⋅m 2 ; b = 2.2 Nm⋅s/rad), while the joint stiffness k(t) was represented by a square with a period of 2 s, offset of 0.5 s and stiffness levels of 50 and 150 Nm/rad (based on: Ludvig et al., 2011;Mirbagheri et al., 2000).The square wave represents the fastest possible transition in joint properties, thereby enabling to reveal the true identification limits of the time-varying methods.
Two commonly applied perturbation input signals were used: (1) a filtered noise signal (filtered with a low-pass 2nd-order Butterworth with a 5 Hz cutoff) with a mean absolute velocity of 0.20 rad/s; and (2) a pseudo random binary sequence (PRBS) signal with a switching rate of 147 ms and a mean absolute velocity of 0.08 rad/s.Both signals were scaled such that they had a root mean square (RMS) value of 0.5 • (0.0087 rad).In addition to the position input, velocity and acceleration input data are required to determine the model's torque output.Therefore, a 40-Hz low-pass analog Bessel filter (2nd-order) was applied to limit the bandwidth of the position input, and subsequently the integrator states were used to extract velocity and acceleration data from the filtered position.This avoided having to compute the derivatives numerically, which is known to be problematic.Output noise n(t) was added as a 40-Hz low-pass filtered (2nd-order Butterworth) normally distributed noise.The amplitude of the output noise was scaled such as to result in a desired SNR of 10 dB, in line with previous research (Golkar et al., 2017;Guarin & Kearney, 2018b).To evaluate the effect of the amount of data collected on the performance of the time-varying methods, simulations were run for a total of 40, 80, 150, 300 and 600 s (sample frequency: f s = 1000 Hz).Each simulation condition was repeated 100 times with different realizations of the input and noise signals.

Preprocessing of simulation data
Before applying any time-varying methods, angle and torque data were 50-Hz low-pass filtered (30th-order finite impulse response filter) and decimated from a sample frequency of 1000 Hz to 100 Hz.No further preprocessing was performed on the simulation data.

Ensemble impulse response function (eIRF)
For the eIRF method the only a priori defined parameter is the memory, or maximum lag (n), for up to which the IRF is estimated.The maximum lag was chosen as 120 ms (12 samples).In the procedure of establishing the time-varying IRF a singular value decomposition is performed of which the optimal model order is based on the Minimum Description Length cost function.The joint stiffness was extracted by taking the integral of the non-parametric impedance time-varying IRF at each time point, effectively calculating the area-under-the-curve.In discrete-time this is formulated as: where n is the number of lags over which the IRF is calculated, H i (i, t) the matrix containing the IRF calculated for each time-point within a realization and Δt the sample time.

Short Data segments (SDS)
For the SDS method the a priori defined parameters are the maximum lag (n) and the number of samples over which the behavior is considered time-invariant (N).The maximum lag was chosen as 120 ms (12 samples).The length of a segment over which behavior was considered timeinvariant was chosen as 100 ms (10 samples).The joint stiffness was extracted by taking the integral of the non-parametric impedance timevarying IRF at each time point, effectively calculating the area-underthe-curve (20).

Basis Impulse response function (bIRF)
The maximum lag was chosen as 120 ms (12 samples) In addition, the time-varying coefficient of the time-varying IRF Λ j (t) were represented using ten cubic B-splines as basis functions.Joint stiffness was estimated using a Bayesian identification algorithm known as relevance vector machine (Bishop, 2006), multiplying the estimated basis function coefficients with their corresponding basis functions.This algorithm assumes a non-restrictive Gaussian prior centered around zero for each of the ten coefficients associated with one of the cubic B-spline basis functions that together describe the joint stiffness.Simultaneously, ten additional coefficients are estimated to model joint viscosity, and a single constant to represent joint inertia (not reported).To update the mean and variance of the distribution associated with each coefficient in the model a direct optimization approach was used together with the position, torque, and calculated velocity and acceleration.

Ensemble spectral method (ESM)
For the ESM method the only a priori defined parameter is the number of samples over which the behavior is considered time-invariant (N).The segment length was chosen to be 30 samples, 300 ms.This number was chosen as a trade-off between the time window considered time-invariant and achieving a sufficient frequency resolution.Joint stiffness is extracted by taking the real value at 0 Hz in the frequency response function.

Kernel based regression (KBR)
For KBR, the parameters N a and N b -determining the system's order were chosen as 0 as 2, respectively.Therefore, we define the transfer function of the time-varying system as: When evaluated in t * , the function H joint (s, t * ) is the transfer function of the LTI system, obtained by fixing the time-varying parameters to their values at time instant t * , i.e., a n (t * ) and b n (t * ).Joint stiffness was extracted from the frozen transfer function by calculating b0 (t * ) a0(t * ) .The hyperparameter p was set to 2 s based on the 0.5-Hz simulated stiffness modulation.The other hyperparameters were l 1 = 1/e 2 s, l 2 = 120, 000 s and γ = 4.5.The parameters were selected as a trade-off between variability of the output signal due to noise, due to the complexity of the model parameters and due to the period-to-period variability of the model parameters.The term f ext (t) was regularized with squared exponential kernel with high smoothness to capture lowfrequency errors.Simulations of 300 and 600 s were excluded from analysis with the KBR method due to the method's high computational demands and limited benefit for the stiffness estimations to include more data.The identification was performed over a bandwidth of 0-20 Hz.

Performance quantification
Five metrics were defined that enabled us to compare the identified joint stiffness to the true simulated joint stiffness (based on Ludvig & Perreault, 2012).First, the bias error described the mean estimation error with respect to the simulated joint stiffness: where K(iΔt) is the simulated stiffness and K(iΔt) is the mean estimated stiffness across all simulated trials at timepoint iΔt.Second, the random error quantified the variance of the estimate across simulation trials (noise sensitivity): where K(iΔt, r) is the stiffness estimated at timepoint iΔt for realization r.Third, the total error described the overall error in the identification of joint stiffness: The fourth performance metric was the rise time.The rise time provided a measure of how quickly a time-varying system identification method can adapt to an instantaneous change in simulated stiffness.The rise time was quantified by the difference between the time points at which the mean estimated stiffness first exceeded 60 and 140 Nm/rad respectively.
The fifth performance metric used was the variance-accounted-for (VAF) to quantify how well the time-varying model of joint impedance explains the output torque data.The VAF was calculated for each 2-s realization (r) according to: where y(iΔt, r) is the simulated (noiseless) or recorded output torque for realization r at timepoint iΔt minus the grand mean across all realizations (the voluntary contribution) and ŷ(iΔt, r) is the estimated output torque at timepoint iΔt based on the identified model(s) of joint impedance.VAF was calculated using the aligned and re-sampled torque data, and additionally when torque data were first band-pass filtered by a 3rd order Butterworth filter with a passband of 1-10 Hz.A passband of 1-10 Hz was used to minimize the contribution of inertia, which generally dominates the VAF using full bandwidth data, and provide a better estimate of how well the data are fitted by the joint stiffness and viscosity.A time-varying VAF was also determined by calculating the VAF across each time sample during a 2-s period of time-varying behavior, and only for the aligned and re-sampled torque data.
If the calculated (time-varying) VAF was less than 0% due to a high error variance, VAF was set to 0%.The mean ± 2 standard deviations (TV-)VAF across all realizations was used to evaluate algorithm performance.

Results
An example of simulated and estimated IRFs for the IRF based methods (eIRF, SDS and bIRF) is provided in Fig. 2. Figs. 3 and 4 present the estimation results for all time-varying system identification methods using filtered noise and PRBS input signals, respectively.Under most conditions, the performance of the algorithms did not depend on the input signal and results only varied slightly.Therefore, further results will be discussed irrespective of the type of input signal.
The eIRF method provided an accurate mean estimated joint stiffness (low bias error), however with substantial variance across simulation trials (high random error).Nevertheless, because the time-varying behavior is determined on a sample-by-sample basis, the instantaneous transition between stiffness levels was well captured (low rise time).In contrast, to the eIRF method, the bIRF method struggled to capture the rapidly varying joint impedance.As a result, the bIRF bias error is largest among all methods.Yet, the bIRF method provides a small random error meaning the method is the least sensitive to the variance across the simulation trials and can handle less data better than the eIRF method.The other methods (SDS, EMS, and KBR) all provide estimates with average bias error, random error and rise time.
The sensitivity for the amount of data used for identification of joint impedance was primarily identified from the random error of the estimated joint stiffness.The eIRF method demonstrated the largest sensitivity for the amount of data, with the random error decreasing substantially when more than 40-s of data were used.Despite the distinct properties of each method, all provided an estimate of joint impedance which resulted in VAFs > 85% for all simulations.Whilst on average the VAF was high, the time-varying VAF across the period of time-varying dynamics demonstrates that the estimates of joint impedance are worse for all but the eIRF method, when the stiffness transitions occur.

Participants
Six healthy participants (range: 23-25 years, 4 women) with no selfreported history of neurological or orthopedic ankle problems, participated in the experiment.All participants provided written informed consent to participate in the study, which was approved by the human research ethics committee (HREC) of Delft University of Technology.The experiment and all experiments was performed in accordance with the Declaration of Helsinki.

Experimental setup
An ankle manipulator (the 'Achilles' -MOOG FCS Inc., Nieuw-Vennep, the Netherlands) enables one degree-of-freedom angular perturbations to the ankle in the sagittal plane (see Fig. 5).Participants were comfortably seated in front of a monitor, with their left foot attached to a rigid footplate using Velcro straps to ensure a firm link between the manipulator and the foot (manipulator stiffness ∼ 3 kNm/rad).The lower leg and foot were positioned such as to achieve 120 • knee and 90 • ankle angles.The monitor in front of the participant was used to present a target torque level, and the ankle torque exerted by the participant.

Measurement protocol
Before starting the experiment, the maximum voluntary torque (MVT) for ankle flexion was determined.Consequently, each participant completed 20 trials during which they were asked to voluntarily modulate their ankle flexion torque by following an onscreen target line (0.5-Hz sinusoid scaled between 5 and 50% MVT).During each trial participants received continuous angular perturbations, the same filtered noise (10 trials) or PRBS perturbation (10 trials) as used in the simulation study, which they were asked to ignore.To aid this instruction, the exerted torque was 0.6-Hz low-pass filtered before being fed back to the participant onscreen.Each trial lasted 70 s, including 5 s at the start and end of each trial without any perturbations.The torque on the manipulator's footplate and angle of the footplate were sampled and   3. Performance metrics obtained using a filtered noise angular input signal in simulation.Each column shows the results for one time-varying system identification method.The top panels show estimated joint stiffness (mean ± 2 standard deviations) for a single 2-s realization of time-varying behavior and corresponding time-varying VAF (mean ± 2 standard deviations) when 150-s of data were used for the system identification.The bar graphs display the performance metrics for the estimated joint stiffness compared to the simulated joint stiffness for all simulation durations.The presented VAF was calculated based on the unfiltered torque data.The KBR method was not analyzed for 300-s and 600-s of data.stored at 1024 Hz.

Preprocessing of experimental data
Before applying time-varying system identification methods, angle and torque data were 64 Hz low-pass filtered using (30th-order finite impulse response filter) and decimated from 1024 to 128 Hz.The first and last 5-s of data were removed from each trial.The data were rearranged to separate the remaining 600-s of data for each perturbation signal into 300 realizations of 2 s.Subsequently, a data alignment, detrending and exclusion procedure was performed (described below) to minimize the variance among the recorded torque realizations at each point in time caused by varying tracking performance within and over trials.After exclusion, all realizations identified as outliers were removed and remaining torque realizations were used for the identification of impedance.Before estimating the impedance, the torque data of each period were detrended by subtracting the mean of remaining realizations (effectively the mean voluntary torque contribution).Bootstrapping with replacement was used to select 75 or 150 realizations of input and output data 100 times.

Alignment of exclusion of realizations
The same alignment and exclusion procedure was followed for both perturbation signal types.First, all recorded realizations of human torque data were aligned to the 0.5 Hz reference torque sine wave that the Fig. 4. Performance metrics obtained using a PRBS angular input signal in simulation.Each column shows the results for one time-varying system identification method.The top panels show joint stiffness (mean ± 2 standard deviations) for a single 2-s realization of time-varying behavior and corresponding timevarying VAF (mean ± 2 standard deviations) when 150-s of data were used for the system identification.The bar graphs display the performance metrics for the estimated joint stiffness compared to the simulated joint stiffness for all simulation durations.The presented VAF was calculated based on the unfiltered torque data.The KBR method was not analyzed for 300-s and 600-s of data.participant had to track.Next, a time shift was sought that minimized the squared error between the recorded torque and torque target: This minimum was sought over a range of time shift from − 0.5-0.5 s and was computed for each 2 s period.Any torque realizations that required a shift magnitude of more than 0.25 s were excluded (number of exclusions based on this criterion for filtered noise: 0-3 realizations; PRBS: 5-40 realizations).Third, the 0.5-Hz frequency component (the frequency of the voluntary applied torque) of the torque was extracted.This was accomplished by removing all frequencies except 0.5 Hz digitally.Any torque realizations that exceeded the mean ± 3 standard deviations of the ensemble of realizations were excluded (number of exclusions based on this criterion for filtered noise: 4-12 realizations; PRBS: 3-10 realizations).Next, the mean voluntary torque contribution was removed by subtracting off the ensemble mean from each realization.The residuals, which consisted primarily of the perturbationevoked torques, were again assessed for outliers.These data were filtered using a 3rd-order 5-Hz low-pass Butterworth filter and realizations were excluded that exceeded the mean ± 3 standard deviations. of all detrended realizations (number of exclusions based on this criterion for filtered noise: 17-29 realizations; PRBS: 15-20 realizations).

Application of time-varying system identification algorithms
For all algorithms, the same settings were used as in the simulation study, unless stated otherwise below.Any differences are driven by the different sample frequencies at which data during simulation and the experiment were acquired, and the different expected behaviorsquare wave versus sinusoidal modulation of joint stiffness.
• For the eIRF, SDS and bIRF methods the maximum lag was chosen 117 ms (15 samples).• For the SDS method the length of a segment over which behavior was considered time-invariant set to 94 ms (12 samples).• The segment length in the ESM method was chosen to be 234 ms ( 30samples).• For the KBR method.the hyperparameters p was set to 2 s based on the 0.5-Hz voluntary torque modulation.In addition, l 1 = 1/e 1 s, l 2 = 30 s and γ = 13.8.Instead of 150-s and 300-s of data, 40-s and 80-s of data (20 and 40 realizations) were used to perform the system identification.

Performance quantification
Joint stiffness estimation performance was described by determining the mean ± 2 standard deviations. of estimated joint stiffness across the 100 bootstrap realizations along with the VAFs.The random error was represented by the standard deviation of the stiffness estimate.Bias and total estimation error could not be determined as for the simulation study, as here the true joint stiffness is unknown.In addition, rise time is not an appropriate measure given no instantaneous changes in stiffness are to be expected.

Results
Data for a representative subject is shown in Fig. 6.After alignment, exclusion and detrending, illustrated in Fig. 7, on average 260 ± 12 realizations (mean ± standard deviations) were included for the filtered  noise perturbation, whereas 255 ± 14 realizations were included for the PRBS perturbation trials.
Fig. 8 provides a comparison of the mean estimated joint stiffness using 150-s of data for all participants and both perturbation signals.
The estimates of stiffness differed amongst the different methods and perturbation signals but demonstrate similar time-varying characteristics.Greater ankle torque is associated with a higher joint stiffness and vice versa.
Further details of the joint stiffness estimate for a single participant are provided in Fig. 9a.The joint stiffness (mean ± 2 standard deviations) is shown for all time-varying system identification methods using both perturbation signals and 150-s of data.The random error (or standard deviation) was found to be in the same order of magnitude for all methods but the KBR.Using the KBR method, the random error is greater when using a PRBS perturbation compared to the filtered noise perturbation.When 300-s of data are used, the random error decreases (Fig. 10b) while the mean estimated joint stiffness remains similar.Overall, the random error is greatest for the eIRF method, and smallest for the bIRF method.
The VAF (Fig. 10a) and time-varying VAF (Fig. 9b) show that using a PRBS input signal results in a 20-30% higher VAF compared to using a filtered noise input signal for all but the KBR method.Using 150-s of data, the mean VAF for the eIRF, SDS, bIRF and ESM method ranged from 74.8 to 87.2% when using a PRBS perturbation and 49.4 to 62.7% when using filtered noise perturbation.Increasing the amount of data to 300 s reduces, or only marginally improves, the VAFs to ranging from 74.7 to 83.0% when using a PRBS perturbation to 49.3 to 57.2% when using filtered noise perturbation.VAFs for the KBR method was 49.6% when using a PRBS perturbation and 67.4% when using a filtered noise perturbation (regardless of the amount of data).However, when filtering the estimated and recorded torque data by a 1-10 Hz band-pass filter, VAF when using the PRBS perturbation is smaller (~10%) than when using the filtered noise.In addition, the time-varying VAF reveals a reduced VAF at the beginning, halfway and end of the realization.

Discussion
In this study we performed a systematic comparison of five timevarying system identification methods that are used to quantify joint impedance while joint position, torque and/or muscle activity are changing.A simulation study was performed to systematically evaluate the performance of each method in estimating joint impedance and reveal each method's strengths and weaknesses.An additional experimental study was used to directly compare estimates of joint stiffness of the methods using human data.Together, the simulation and experimental study revealed distinct properties of each method which can be quantified using a comprehensive set of performance metrics.This result emphasizes the importance of careful selection of the best method for each application and calls for elaborate benchmarking of every newly developed method in the future.

Comparison of time-varying system identification methods to assess timevarying ankle impedance
The results demonstrate the unique properties of the five timevarying system identification methods.An overview of these properties is provided Table I.The ensemble based eIRF method provides the highest VAFs and smallest total joint stiffness estimation error, by providing time-invariant estimation of joint impedance on a sample-bysample basis (Lortie & Kearney, 2001).However, this can only be achieved given that a few hundreds of realizations under realistic noise conditions are available.Recording hundreds of realizations with consistent behavior may be an issue in experiments with healthy participants, where variability is at the heart of our motor behavior, and Fig. 7. Illustration of the data alignment, detrending and exclusion procedure using a representative data set (input signal: filtered noise).The top row shows all 300 realizations of the recorded ankle torque (left) and all aligned and detrended torque realizations used for time-varying system identification (right).The bottom row visualizes the intermediate processing step.Step 1: All torque realizations are aligned to the participants target torque (black line) -and realizations excluded which require a time shift of more than 0.25 s.
Step 2: The 0.5 Hz torque contribution was extracted and realizations exceeding mean ± 3 standard deviations (indicated by the solid and dashed black lines) of all realizations excluded.Step 3: All realizations were mean subtracted and realizations exceeding mean ± 3 standard deviations (indicated by the solid and dashed black lines) of all realizations excluded.processes like attention and muscle fatigue play a major role (Zhang & Rymer, 2001), let alone in participants with a compromised motor system.The estimation performance of the eIRF method quickly drops with reduced amounts of data or with increased variance in time-varying behavior across realizations, making all other methods more appropriate when a large dataset cannot be obtained.
The ensemble based SDS, bIRF and ESM methods reduce sensitivity for the amount of data by adding direct or implicit assumptions about the time-varying dynamics.The SDS and ESM method consider behavior to be time-invariant during a short time-window, whereas the bIRF uses (time-varying) basis functions to model the time-varying joint impedance (Guarin & Kearney, 2017, 2018a;Ludvig & Perreault, 2012).The reduced need for data compared to the eIRF method comes at the cost of losing the ability to capture fast changes in joint impedance.Inappropriate definition of window size (SDS, ESM) or choice of basis function (bIRF) directly affects the maximum rate of change in joint impedance that can be identified.Hence, the identified time-varying joint impedance is primarily driven by the dynamics that can be captured across the windows or by the (sum of) the basis functions.For example, the cubic B-splines used in the bIRF method are suboptimal to capture the rapid change in joint stiffness in the simulation study but fit the expected periodic behavior in the experimental study well (see Appendix A for simulation results of sinusoidal time-varying behavior).Identifying rapid changes in joint impedance would require the use of Haar wavelets which are square-shaped functions that would capture the simulated steps in joint impedance better (Guarin & Kearney, 2018b).
In the context of data requirements, the KBR method is truly unique.The time-varying joint impedance could be captured across single realizations of data.However, as the method imposes an approximate periodicity of the time-varying behavior, additional realizations will improve the estimate as information can transfer from one period to the other.The number of realizations still improving estimator accuracy is determined by the width of the kernel (parameter l 2 ).The resulting estimated model is an analytical model which can easily be used for further simulation.Unfortunately, the identification procedure is computationally heavy and requires optimization and assumptions of many underlying parameters.The model orders (N a and N b ) and hyperparameters need to be determined a priori, which is a non-convex mixed-integer optimization problem.Hence, simply assuming a secondorder system to represent joint impedance as often done, may not suffice (Sobhani Tehrani et al., 2017).
The methods compared in this study were chosen based on their availability and primary purpose: to identify joint impedance.Methods not used include the (short segment) structural decomposition subspace method (SDSS) (Jalaleddini et al., 2017a(Jalaleddini et al., , 2017b)), the subspace, linear parameter varying, parallel cascade (LPV-PC) method (Golkar et al., 2017), skirt decomposition method (Lataire et al., 2012;van de Ruit et al., 2020) and a time-frequency method together with modal analysis (Piovesan etal., 2009(Piovesan etal., , 2012)).The SDSS and LPV-PC method are specifically designed to separate intrinsic and reflexive contributions to joint impedance which is outside the scope of this work.The SDSS method allows to directly obtain a parametric description of the joint parameters without making a priori assumptions about the system order, as the order is estimated as part of the method.The method is currently only capable of capturing slow changes in joint impedance, which contrasts the LPV-PC method, that can capture rapid changes in joint impedance.The LPV-PC method requires a scheduling function a priorirepresenting the assumed time-varying behaviorto achieve this.The correct scheduling function may be difficult to obtain and may change throughout an experiment.The skirt decomposition method was not included as it is only applicable when using sum-of-sine perturbations signals whilst the   time-frequency method with modal analysis does work with transient (impulse) perturbations rather than continuous perturbations.

Factors not considered in this study
Both the simulation and experimental study used an open-loop design.The simulations and experiments used an angular perturbation as an input signal where the resulting torque does not affect the joint angle.This does not accurately reflect time-varying joint impedance during our daily life activities, where most of our movements are made while our joints interact with a compliant environment.This effectively means that any torque generated in response to an angular perturbation will result in a change in joint angle.The presence of this interaction means that joint angle and torque are measured in a closed loop system (de Vlugt et al., 2001;van der Helm et al., 2002).Application of the open loop methods used in this study to data recorded in closed-loop conditions would result in biased parameter estimates (Kearney & Hunter, 1990;Ludvig & Kearney, 2009).All methods used in this study can be rewritten to be applicable to closed-loop data (e.g., for the SDS method: Ludvig & Kearney, 2009;Ludvig et al., 2017).While the modifications necessary to the methods to enable use of closed-loop data are not expected to have a substantial effect on the unique properties of the methods as presented for open-loop data in this paper, we are unable to conclude how each of these methods would perform when modified to work in closed-loop situations.
In this study we focused on the identification of joint stiffness, without separating intrinsic and reflexive contributions.Reflexes may contribute up to ~50% to the total joint stiffness, depending on the muscle contraction level, perturbation characteristics and the joint studied (Mrachacz-Kersting & Sinkjaer, 2003;Sinkjaer et al., 1988;Stein & Kearney, 1995;Thomson & Chapman, 1988).Hence, during the experimental study, the time-varying ankle torque will have resulted in a time-varying reflex contribution within the trials.Moreover, the difference in mean absolute velocity between the perturbation signals may have resulted in a different reflex contribution across trials.Reflexes have been shown to reduce with mean absolute velocity and therefore, trials with the PRBS perturbation may have a smaller contribution of reflexes due to the lower mean absolute velocity of the perturbation signal (Stein & Kearney, 1995).The identification of time-varying reflexes requires more advanced and complex techniques, that may or may not require the recording of electromyography (EMG) from the muscles around the joint and additional assumptions on unknown parameters that describe the nonlinear reflex dynamics.Whilst an elegant and successful parallel cascade model has been described to distinguish (time-varying) intrinsic and reflex contributions to joint stiffness using ensemble time-domain based methods (eIRF, SDS and bIRF) (Golkar et al., 2017;Guarin & Kearney, 2017;Jalaleddini et al., 2017;Kearney et al., 1997;Ludvig et al., 2011), this remains an unsolved problem for the frequency-domain methods used in this study.Although, reflex and intrinsic contributions were left unseparated to enable comparison across all methods, one should be aware that any reflex activity present in the current data may have biased the stiffness estimates (Guarin & Kearney, 2018b).
Several additional study design decisions were made to keep the work manageable.Joint inertia and joint viscosity were considered time-invariant (simulation study) or unidentified (experimental study).We chose to only vary and estimate joint stiffness as stiffness modulation is most pronounced and does not require additional parameter estimation procedures having their own pitfalls.Joint inertia is generally found to vary little across movement phases, whereas joint viscosity co-varies with joint stiffness during movement and between joint operating points (Bennett et al., 1992;Lee & Hogan, 2015;Mirbagheri et al., 2000).Nonetheless, the identification of inertia and viscosity could be done by fitting a parametric model to the time-varying impedance IRF (eIRF, SDS) or frequency response function (ESM, KBR) and by Bayesian identification (bIRF).In simulation, we confirmed that covarying joint viscosity with joint stiffness at an order of magnitude of 1-3 Nms/rad neither does affect the identified properties of the methods nor presents different results for the estimation of joint stiffness.In addition to focusing on joint stiffness, we only compared the performance of the methods using one SNR level (10 dB) motivated by findings from experimental data (Mirbagheri et al., 2000).Whereas different SNR levels were not compared, the sensitivity of the methods to different levels of SNR can be inferred from the different amounts of data used in the simulation study performed.The use of more data for identification, reduces noise due to averagingimplicit to the ensemble-based methods.

Use the VAF with care
During this study, some important practical lessons were learned.In the experimental study, the VAF of the recorded torque and the standard deviation of the estimated joint stiffness were used to validate the identification results.The VAFs of the experimental study (Figs.9b and  10a) would suggest a better estimation performance of the method using a PRBS rather than a filtered noise perturbation.This can be explained by the larger high frequency content of the PRBS signal, which gets amplified in the measured joint torque due to the 2nd order nature of joint impedance and the inability of the human to effectively suppress these higher frequencies.The VAF is therefore largely determined by the ability of the model to correctly estimate the highest frequencies and as the torque response to a PRBS perturbations contains greater power at these higher frequencies, the VAF will always be higher.A solution may be found by removing frequencies > 10 Hz (the ankle's natural frequency), thereby reducing the contribution of inertia on the VAF.This VAF better reflects how well the method estimates the signal content determined by the joint viscosity and stiffness.
There is the seemingly contradictory finding in the experimental • No a priori assumptions on system dynamics or expected time-varying behavior are required • Many realizations are required with the same underlying time-varying behavior.

SDS
• Time-varying behavior is assumed stationary across realizations.
• Time-varying behavior is considered time-invariant across small time window -which influences the fastest dynamics that can be captured.
• No a priori assumptions on system dynamics or expected time-varying behavior are required.bIRF • Time-varying behavior is assumed stationary across realizations.
• Time-varying behavior is assumed to be represented by a set of a priori chosen basis functions which dictates the dynamics of the captured behavior.ESM • Time-varying behavior is assumed stationary across realizations.
• Time-varying behavior is considered time-invariant across small time window -which influences the fastest dynamics that can be captured.KBR • Time-varying behavior is assumed to be represented by a set of a priori chosen basis functions which dictated the dynamics of the captured behavior.• Order of system dynamics should be chosen a priori.compared to the simulation data that use of more data does not improve VAF.In the simulation data the VAF is calculated between the estimated torque output and simulated, noiseless, torque output.Therefore, the VAF will go towards 100% when more data are used.In case of the experimental data, we only have the measured (noisy) torque recordings because of which the VAF is limited by the SNR of the recorded data.More data may thereby not change the VAF much whereas the random error (equivalent to the standard deviation) of the estimated joint stiffness is markedly reduced.Hence, the VAF is a poor measure to draw conclusions on the amount of data required to get a reliable estimate of the joint's properties estimated, and additional measures are required to reflect how well a system identification method estimates joint properties.
The time-varying VAF demonstrates worse VAF where a rapid change in dynamics occurs (simulation data) or at the start, halfway and end of a realization (experimental data) when ensemble-based methods are used.What all these time points have in common is the strong rate of change in joint dynamics that takes place.Hence, a stronger rate of change of joint dynamics is associated with a larger variance across realizations and thereby results in worse identification results.This also stresses the importance of appropriate alignment procedures to ensure maximal overlap in underlying joint impedance between realizations, while minimizing the variance at every time point.

An assessment framework for time-varying system identification methods
A uniform evaluation of time-varying system identification method properties facilitates comparison, especially for those that wish to apply time-varying system identification as a tool without being experts in the field.Therefore, we suggest that authors that develop a new timevarying system identification method follow our framework in which a simulation and experimental study are employed to demonstrate the strengths and weaknesses of their methods.The simulation study reveals sensitivity to the amount of data used for identification and provides a quantitative description of properties and estimation errors.Moreover, by including extreme conditions e.g., a square wave time-varying stiffness, the true strength and weaknesses of the method can be revealed.Next, the experimental study should be employed as a use case, with actual recorded human data.Ideally, in between the simulation and experimental study, the method is also validated using a mechanical setup with dynamical properties that mimic human joint properties.In addition, we encourage authors to publish datasets which can serve as a benchmark dataset for their newly developed method(s).

Conclusion
The use of (time-varying) system identification provides valuable insights into human movement control (Vlaar & Schouten, 2015) and helps the design of biomimetic prosthetics and orthotics (Hansen et al., 2004) which is of crucial importance to their users (Shepherd et al., 2018).The continuous development and improvement of system identification methods has helped tremendously to aid reliable description of human joint dynamics despite the noisy data recorded during human experiments.Nonetheless, the system identification community needs to recognize that to enable widespread use of their methods, code needs to be made easily accessible, and the method's strength and weaknesses clearly defined.This work serves to provide a clear overview of the unique properties of different time-varying system identification methods available, and the importance of an elaborate description of their properties using bias and random error, and maximum rate of change in joint dynamics.This allows researchers to make a well-justified decision on which method would be most appropriate for their application based on which a priori assumptions can be confidently made and what performance is desired.

Declaration of Competing Interest
The authors declare no conflict of interest.

Appendix A
To facilitate interpretation of the experimental results, where sinusoidal time-varying joint stiffness is expected, additional simulation trials where performed.

Methods
The same anti-causal time-varying open-loop model of human joint impedance and angular PRBS or filtered noise input were used as in the main simulation study.Joint inertia and viscosity remained time-invariant (I = 0.02 kg⋅m 2 ; b = 2.2 Nm⋅s/rad), while the joint stiffness k(t) was represented by a sinusoid with a period of 2 s mean stiffness of 100 Nm/rad and amplitude of 50 Nm/rad.These parameters reflect the time-varying behavior expected from the experimental study.The amplitude of the output noise was scaled such as to result in a desired SNR of 10 dB.Simulations were run for a total of 150 s (sample frequency: f s = 1000 Hz), equivalent to the amount of data used in the experimental study.Each simulation condition was repeated 100 times with different realizations of the input and noise signals.All data processing and time-varying system identification algorithms were applied as described in the simulation study.

Fig.
A.1 presents the estimation results for all time-varying system identification methods using filtered noise and PRBS input signals, respectively.All methods provide a good estimate of the time-varying joint stiffness, with the greatest bias found for the bIRF method and greatest variance across simulation trials for the eIRF method.The lowest bias error is found for the SDS and KBR method, whereas the lowest random error is observed for the bIRF and KBR method.VAFs are > 90% regardless of the method used, with the time-varying VAF showing greater variance for the PRBS input signal and lower joint stiffness values.Performance metrics obtained when simulating sinusoidal time-varying joint stiffness using a filtered-noise and PRBS angular input signal.Each column shows the results for one time-varying system identification method.The top two panels show estimated joint stiffness (mean ± 2 standard deviations) for a single 2 s realization of time-varying behavior using each input signal when 150-s of data were used for the system identification.The middle two panels show the corresponding time-varying VAF (mean ± 2 standard deviations).The bar graphs display the performance metrics for the estimated joint stiffness compared to the simulated joint stiffness for both input signals.The presented VAF was calculated based on the unfiltered torque data.

Fig. 1 .
Fig.1.The implemented time-varying model of human joint impedance H joint (s, t) with angular input u(t), torque output y(t), and additive torque output noise n(t).H joint (s, t) was simulated as a 2nd order model with time-varying stiffness k(t), modelled as a 0.5 Hz square wave.The input signal was either a 5 Hz low-pass filtered noise u 1 (t) or pseudorandom binary sequence u 2 (t).

Fig. 2 .
Fig. 2. Simulated IRF (gray line) and estimated IRF for the IRF based methods (eIRF: blue dashed line; SDS: red dashed line and bIRF: green dashed line) using 150-s of data and a PRBS input signal for one time point at which the joint stiffness (k) is low (50 Nm/rad) and one at which stiffness is high (150 Nm/ rad) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

Fig.
Fig.3.Performance metrics obtained using a filtered noise angular input signal in simulation.Each column shows the results for one time-varying system identification method.The top panels show estimated joint stiffness (mean ± 2 standard deviations) for a single 2-s realization of time-varying behavior and corresponding time-varying VAF (mean ± 2 standard deviations) when 150-s of data were used for the system identification.The bar graphs display the performance metrics for the estimated joint stiffness compared to the simulated joint stiffness for all simulation durations.The presented VAF was calculated based on the unfiltered torque data.The KBR method was not analyzed for 300-s and 600-s of data.

Fig. 5 .
Fig. 5. Schematic of the lower leg and foot positioned such to achieve a 90-deg ankle angle with the robotic manipulator 'Achilles'.Recorded signals are the angular input (ϴ) and torque output (T).The inset presents the onscreen feedback participants received to perform the instructed task: apply a timevarying ankle torque.

Fig. 6 .
Fig.6.Representative 10-s of input and output y(t) data from a single participant for both the filtered noise (u 1 (t) and y 1 (t) -top row) and pseudorandom binary sequence (u 2 (t) and y 2 (t) -bottom row) perturbation signals.

Fig. 8 .
Fig. 8. Mean estimated joint stiffness of all participants (P1-P6) and all timevarying system identification methods when 150 s of data were analyzed.

Fig. 9 .
Fig. 9.Estimated joint stiffness and time-varying VAF for all time-varying system identification methods using 150-s of data from a single participant.Each column represents the results for one time-varying system identification method.The gray line highlights the requested 0.5-Hz sinusoidal time-varying torque behavior the participant was requested to perform.(a) Estimated joint stiffness (mean ± 2 standard deviations) (filtered noisetop panel; PRBSbottom panel).(b) time-varying VAF (mean ± 2 standard deviations).

Fig. 10 .
Fig. 10.VAF and random error using 150-s or 300-s of data and a filtered noise (top panel) or PRBS (bottom panel) perturbation signal (a) VAF (mean ± 2 standard deviations).for all participants) of all time-varying system identification methods.VAF was calculated based on the full bandwidth recorded and estimated torque output (left), and band-pass filtered torque data between 1 and 10 Hz (right).(b) Estimated random error of the estimated ankle stiffness (mean ± 2 standard deviations for all participants).

Fig
Fig. A.1.Performance metrics obtained when simulating sinusoidal time-varying joint stiffness using a filtered-noise and PRBS angular input signal.Each column shows the results for one time-varying system identification method.The top two panels show estimated joint stiffness (mean ± 2 standard deviations) for a single 2 s realization of time-varying behavior using each input signal when 150-s of data were used for the system identification.The middle two panels show the corresponding time-varying VAF (mean ± 2 standard deviations).The bar graphs display the performance metrics for the estimated joint stiffness compared to the simulated joint stiffness for both input signals.The presented VAF was calculated based on the unfiltered torque data.

Table 1
Overview of assets and assumptions of the five time-varying system identification methods compared in this study.Time-invariant estimate of joint dynamics is obtained at each sample time.• Time-varying behavior is assumed stationary across realizations.