PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Biomech Eng. Author manuscript; available in PMC 2010 August 30.
Published in final edited form as:
PMCID: PMC2929930
NIHMSID: NIHMS134713

Optimal Estimation of Dynamically Consistent Kinematics and Kinetics for Forward Dynamic Simulation of Gait

Abstract

Background

Forward dynamic simulation provides a powerful framework for characterizing in vivo loads, and for predicting changes in movement due to injury, impairment or surgical intervention. However, the computational challenge of generating simulations has greatly limited the use and application of forward dynamic models for simulating human gait.

Methods

In this study, we introduce an optimal estimation approach to effciently solve for generalized accelerations that satisfy the overall equations of motion and best agree with measured kinematics and ground reaction forces. The estimated accelerations are numerically integrated to enforce dynamic consistency over time, resulting in a forward dynamic simulation. Numerical optimization is then used to determine a set of initial generalized coordinates and speeds that produce a simulation that is most consistent with the measured motion over a full cycle of gait. The proposed method was evaluated with synthetically created kinematics and forceplate data in which both random noise and bias errors were introduced. We also applied the method to experimental gait data collected from five young healthy adults walking at a preferred speed.

Results and Conclusions

We show that the proposed residual elimination algorithm (REA) converges to an accurate solution, reduces the detrimental effects of kinematic measurement errors on joint moments, and eliminates the need for residual forces that arise in standard inverse dynamics. The greatest improvements in joint kinetics were observed proximally, with the algorithm reducing joint moment errors due to marker noise by over 20% at the hip and over 50% at the low back. Simulated joint angles were generally within 1 deg of recorded values when REA was used to generate a simulation from experimental gait data. REA can thus be used as a basis for generating accurate simulations of subject-specific gait dynamics.

Keywords: Gait, Joint Moments, Residual Forces, Noise Rejection, Numerical Integration, Initial Conditions, Geometric Constraints

I. INTRODUCTION

Forward dynamic models provide a powerful framework for investigating the causal relationship between muscle forces and the movement generated in both normal and pathological gait [1]. In addition, dynamic models are predictive and thus can be used to simulate how movement will change as a result of an intervention that affects the system or neuromuscular control [2]. Despite the benefits, forward dynamics is not nearly as commonly used as inverse dynamics analysis. This is, in part, attributable to the inherent difficulties in generating coordinated simulations of human movement driven at either the muscle or joint level

Traditional inverse dynamics approaches employ an iterative solution of the equations of motion. In gait analysis, measured ground reactions and kinematics are used to calculate proximal joint loadings at each segment, starting at the foot and proceeding proximally [3]. However, this iterative analysis results in a set of residual forces and moments acting on the last segment in the chain that have no physical meaning, i.e. there is no subsequent segment or joint that can produce them. Thus, it is difficult to use inverse dynamics results as a basis for generating a forward simulation without including these non-physical residuals. A promising approach for eliminating residuals is to perform least squares inverse dynamics, in which a set of joint moments are determined that are simultaneously most consistent with measured kinematics and ground reactions [46]. While this approach has been shown to improve joint moment calculations, the estimates are not dynamically consistent over time and thus the resulting accelerations, when integrated, will not reproduce the desired motion.

Dynamic optimization can be used to determine joint moment trajectories or muscle excitation patterns that produce motion that is dynamically consistent over time [7, 8]. However, dynamic optimization suffers from practical limitations. The objective functions used often exhibit multiple local minima, and the convergence behavior and the quality of the solution depend strongly on the accuracy of the initial guess. Furthermore, the inordinate amounts of computation time [9] make the approach often impractical; especially, when data of multiple subjects and trials has to be processed.

The aim of this study was to generate dynamically consistent kinematics and kinetics that satisfy the whole body equations of motion and best agree with measured motion and external forces. We show, using both synthetic and experimental gait data, that the proposed algorithm is able to eliminate residual forces while improving the estimation of joint moments. Thus, the method provides a systematic, automated approach for generating forward dynamic simulations that closely replicate subject-specific gait patterns.

II. METHODS

A. Overall Equations of Motion

The residual elimination algorithm (REA) uses least-squares estimation to solve the overall constrained equations of motion for the generalized accelerations at each time step. The estimated accelerations are then integrated to determine the generalized coordinates and speeds which are used in subsequent time steps. Numerical optimization is employed to determine the initial states for this integration so that the resulting motion is dynamically consistent over time and optimally reproduces a desired movement. Following is a detailed description of the formulation.

The Newton-Euler equations of motion require that the sum of all external forces Fi balances the sum of the mass-acceleration products of the individual body segments [10].

i=1NFi=j=1nmj(r¨jg)
(1)

Similarly, the equations for rotational motion about the origin of the coordinate system can be stated as:

i=1N(ρi×Fi+Mi)=j=1nrj×mj(r¨jg)+j=1n[Ijω˙j+ωj×Ijωj]
(2)

Fi and Mi are external forces and moments respectively, ρi is the location of a point along the external reaction force Fi, rj points to the center of gravity, ωj is the angular velocity, and mj and Ij are the mass and the inertia dyadic of segment j. Gravity is denoted by the vector g. The number of segments is n, the number of external reactions is N. Each external reaction inlcudes three moments and three forces. Thus, in the analysis of gait, N = 1 for single support and N = 2 for double support. Note that the unknown joint moments do not appear in these overall equations of motion, which has computational advantages in comparison to other methods [4, 5] that utilize the full set of equations of motion.

B. Cartesian and Generalized Coordinates

The equations of motion stated above need the individual segment positions, velocities, and accelerations in Cartesian coordinates. However, the current state of a skeletal model is usually described in generalized coordinates q and their derivatives q. and q. These generalized coordinates (that mostly correspond to joint angles) can be mapped to the Cartesian coordinates by a nonlinear transformation t.

[rθ]=t(q)
(3)

The instantaneous relationship between Cartesian speeds and the generalized speeds can be expressed with the Jacobian J. (Note: As the J is depending on the current state q, it would be more appropriate to state J(q) )

[r˙ω]=Jq˙
(4)

Differentiating this expression leads to:

[r¨ω˙]=Jq¨+J.q˙
(5)

C. Linear Equations of Motion in Generalized Coordinates

Equations (1) and (2) can be rearranged and put into linear matrix form.

[m100[1][0]0m10[0]00m10ρ1zρ1y0r1zm1r1ym1ρ1z0ρ1x[1]r1zm10r1xm1I1ρ1yρ1x0r1ym1r1xm10][F1M1r¨1ω˙1]=j=1n[mjgmj(rj×g)+ωj×Ijωj]
(6)

or shorter

Acart(r)[FMr¨ω˙]=fcart(r,ω)
(7)

Where Acart is a matrix of dimension 6 × (n · 6 + N · 6) and fcart(r,ω) a six dimensional vector function. This equation states that the external forces and segment accelerations have to balance the gyroscopic and gravitational forces defined in fcart(r,ω). They are coupled by the Newton-Euler equations that impose a dynamic constraint between accelerations and external reactions. The residual forces that arise during inverse kinematics analyses are a consequence of the violation of this constraint by erroneous measurements and/or model inaccuracies. We will refer to this equation later, when geometric constraints are introduced to handle prescribed Cartesian accelerations. However, for the time being, we focus on the dynamic constraints, which can be satisfied by adaptation of the measured generalized accelerations and external reactions. Using Eq. (5), the constraints can be mapped into the generalized coordinate system.

[FMr¨ω˙]=[FMJq¨+J.q˙]=[[1]000[1]000[J]][FMq¨]+[00000000[J.]][00q˙]
(8)

which can be used to state the equations of motion as a linear function of the generalized coordinates:

A(q)[FMq¨]=f(q,q˙)
(9)

with

A(q)=Acart(r)[[1]000[1]000[J]]
(10)

and

f(q,q˙)=j=1n[mjgmj(rj×g)+ωj×Ijωj]Acart(r)[00000000[J.]][00q˙]
(11)

where the Cartesian coordinates and speeds have been calculated from the generalized co-ordinates using

[rθ]=t(q)and[r˙ω]=Jq˙.

The matrix A is of the dimension 6×(N · 6 + nq), where nq denotes the number of generalized coordinates q.

D. Variation Formulation

As a result of model and measurement errors, the experimental recordings of ground reactions and accelerations will generally not satisfy Eq. (9). Small variations (δq, δF , δM) can be introduced to account for the differences between model, measurement and reality, to fulfill the dynamic constraints, and to eliminate the undesired residual forces at the base segment:

F=F+δFM=M+δMq¨=q¨+δq¨
(12)

where F, M and q¨ represent the desired ground reaction forces, moments and generalized accelerations. The desired ground reactions can simply be set to the net experimentally recorded ground reaction forces (F) and moments (M). The desired accelerations include adjustments to the experimentally derived accelerations (q¨) that account for current errors between the experimental generalized speeds (q˙) and coordinates (q) and the corresponding simulated values:

q¨=q¨+kv(q˙q˙)+kp(qq)
(13)

where kv and kp represent feedback gains. This feedback is important to prevent drift from the desired coordinate trajectories, particularly for low inertia segments whose dynamic behavior does not contribute substantially to the ground reactions.

With the variations expressed in (12), equation (9) can be written as

A(q)[δFδMδq¨]=f(q,q˙)A(q)[FMq¨]
(14)

or simply as the linear equation

Aδ=b
(15)

E. Satisfying the Variation Formulation

As Eq. (15) is an under-determined system of linear equations, generally no unique solution exist to this problem. We introduce a diagonal matrix W of measurement standard deviations, which represents the accuracy of the experimental estimates. Quantities that can be measured more precisely will obtain smaller values than quantities that are prone to errors. With this matrix Eq. (15) can be expanded to:

AWW1δ=b
(16)

Using the right Moore-Penrose Matrix Inverse [11] or ‘Pseudo Inverse’ B+ = BT(BBT)−1, Eq. (16) can be solved for a set of deltas

δ=W(AW)+b,
(17)

that will minimize the Euclidean norm W1δ. As a high measurement standard deviation will lead to a small entry in W−1, it will allow a larger value for the corresponding δ-entry and result in larger changes to the original measurement.

It is important to note that this way of computing consistent accelerations and external reactions differs significantly from a preliminary published version of the REA algorithm [12], where only the six pelvic degrees of freedom were altered.

F. Integration and Initial Conditions

The preceding processing of the kinetic and kinematic data was expressed on a per-frame basis and entirely in terms of accelerations. The Newton-Euler equations return no information about position and velocity. On the contrary, these values are necessary in Eq. (15) to evaluate the matrix A, and to determine the current gravitational and gyroscopic forces. The only feasible way of obtaining these values is via numerical integration of the computed generalized accelerations. Unfortunately, the alterations made to the measured accelerations will cause the integrated motion to deviate from the recorded trajectories - an effect that is further augmented by the inherently unstable dynamics of gait. To circumvent this, numerical optimization is used to determine a set of initial generalized coordinates and generalized speeds that, after solving Eq. (15) for accelerations and numerically integrating, produced a motion that best replicated the recorded motion (Fig. 1).

FIG. 1
Schematic of the Residual Elimination Algorithm. Dynamically consistent accelerations are estimated based on the measured data and numerically integrated. The initial conditions for this integration are optimized to yield maximal agreement of simulated ...

The agreement of recorded and calculated motion can be assessed by the squared and weighted distances between calculated and experimentally derived generalized coordinates:

e(t)iq=wiq(q(t)iq(t)i)2
(18)

The weights wqi in this error function are necessary to define the relation of plevic translation errors, rotation errors, and joint angles. They can also be used to improve the tracking of selected generalized coordinates. Alternatively to expressing this agreement in terms of generalized coordinates, it is possible to define it in Cartesian space. A promising approach in this context is the summation of the euclidian distances between predicted and measured marker positions:

e(t)ir=wirr(t)ir(t)i2
(19)

Tracking markers directly has the advantage of removing the intermediate inverse kinematics step for comparison and closing the bridge to the experimental data. The weights wri in the error function can be used to emphasize the spatial tracking of certain markers.

Both error formulations can be used simultaneously. Integrated over the entire motion, these values express a cost function that maps a set of initial conditions to a measure of ‘agreement’ between the calculated and recorded motion. The complete optimization problem can be stated as follows:

minE(qo,q˙o)=totendi=1nqe(t)iq+i=1nre(t)irdts.t.q¨=q¨+δq¨withδ=W(AW)+b
(20)

G. Removing Measurement Offsets

We expanded the search space to include any static offset between the reference frame of the forceplates and the optical tracking system. This was done by adding a constant offset in each direction to the marker kinematics prior to cost function. These three offset parameters were included in the optimization search space when minimizing Eq. (20).

H. Evaluation with Synthetic and Experimental Data

For all analyses, a 3-dimensional full body model consisting of 12 independent segments was used. The generalized coordinates q corresponded to the 27 degrees of freedom that completely described the model’s state (Tab. I). Three different data sets were used in conjunction with this model for evaluation of REA:

  1. Experimental Data Whole body kinematics were collected on five healthy young females (age 25±2 yrs, height 164±4 cm, mass 57±5 kg) for a full cycle of preferred speed (1.41±0.10 m/s) gait using an 8 camera motion capture system (Motion Analysis Corporation, Santa Rosa, CA, USA). Subjects had no history of major orthopedic diagnosis or pain in the lower back, pelvis, or lower extremity. Each subject gave informed consent in accordance with a protocol approved by the University of Wisconsin’s Health Sciences Institutional Review Board. Ground reaction forces were collected simultaneously with three force plates embedded in the lab floor (Advanced Mechanical Technologies, Newton, MA). The model’s segment lengths and anthropometric parameters were scaled to the estimated joint-to-joint lengths of individual subjects and the subject’s overall height and body mass [13]. The locations of body fixed markers were adjusted with respect to an upright standing calibration trial. Joint axes orientation were not altered. A global optimization based inverse kinematics approach [14] was used at each time step to compute the generalized coordinates of the model that best fit measured locations of 38 anatomical and tracking markers (Fig. 2). The resulting generalized coordinate trajectories were then low-pass filtered (6 Hz), spline-fitted, and the corresponding accelerations determined from the second spline derivative. Force plate data was low-pass filtered at 50 Hz.
    FIG. 2
    The residual elimination analysis was used to determine initial states that minimized the sum of weighted squared distances between model predicted and measured marker positions. The 38 markers that contributed to this objective function were the same ...
  2. Synthetic Data with Noise A forward dynamic simulation of a full cycle of gait was used to synthetically create marker and force plate data. The simulation was based on the experimentally recorded motion of a young healthy adult [12]. The model used for simulation was identical to the one used for analysis. It contained no residual forces and the joint moments were known, thereby providing a ground truth by which to compare the results. The kinematics of the 38 anatomical and tracking markers (Fig. 2) were extracted and normally distributed unbiased white noise was added in all three spatial directions. Noise levels (standard deviation) varied from 3 mm to 15 mm by increments of 3 mm. Ground reactions were not altered. 10 trials were performed at each noise level.
  3. Synthetic Data with Noise and Offset In the third data set, constant offset values were added to all spatial directions of the synthetically created marker data. The offset that emulated calibration errors of the optical tracking system- varied from 0 mm to 25 mm by increments of 5 mm. Additionally, a 9mm sd noise was added to all markers. 10 trials were performed at each offset level. Within REA, measurement standard deviations of external forces and moments were assumed to be negligible (1e-7N and 1e-7Nm respectively). The standard deviation of measured generalized accelerations were considered equal for all joint angles and set to 1rad/s2. As the force plate data allows predominantly conclusions about the acceleration of the base segment, pelvic sd’s where set to higher values to allow bigger variation for the underlying coordinates (10 rad/s2 for rotational accelerations and 1 m/s2 for translational accelerations of the pelvis).
TABLE I
The 3-dimensional full body model consisted of 12 independent segments (Pelvis, thighs, shanks, feet, trunk, upper arm and forearms). The model state was described by 27 generalized coordinates q.

Optimization of the initial conditions and offset values was first performed for the entire set of generalized coordinates and speeds, resulting in 59 optimization variables. To compare accuracy and computational performance, in a second analysis only the offset values and the initial conditions of the 6 pelvic coordinates were optimized and the remaining generalized coordinates and generalized speeds remained fixed to the outcome of the inverse kinematics routine. In both cases, the trajectories of all 38 markers contributed to the objective function of the numerical optimization. Higher weights (w = 5) were put on the anatomical markers on the pelvis and lower extremity. Lower weights (w = 1) were used for the torso markers, upper extremity markers and lower extremity tracking markers (Fig. 2). Feedback gains kp = 100 and kv = 20 were used. Constant offsets between the kinematic and forceplate reference frames were estimated. Finally, inverse dynamic analysis was performed with the original and with the REA generated kinematics to determine the effect of the processing on the computed joint moments.

The data processing was implemented in C, using the Dynamics Pipeline of SIMM (Motion Analysis Corporation, Santa Rosa, CA, USA). The equations of motion were derived with SD/FAST (Parametric Technology Corporation, Waltham, MA). All numerical optimizations were performed using a sequential quadratic programming method (FSQP; AEM Design, Tucker, GA).

III. RESULTS

A. Experimental Data

When performing inverse dynamics analysis on the raw experimental data, substantial residual forces were present with average rms values of 14-25 N in translation and 4-18 Nm in rotation (Tab. II). The proposed REA algorithm eliminated the residual forces and moments by estimating generalized accelerations that balanced the overall equations of motion (Fig. 3). After optimizing for the best initial generalized coordinates and speeds, the mean kinematic changes introduced by REA were generally less than 1 deg for pelvic orientation and joint angles (Tab. III). The greatest discrepancy was with pelvic transverse rotation, with an average difference of 1.3 (± 1.0) deg.

FIG. 3
Traditional inverse dynamics analysis produces substantial residual forces and moments (solid lines) that have no physical meaning, as shown for this sample experimental data set. The proposed algorithm eliminates these residuals completely (dashed lines) ...
TABLE II
Average (± 1 s.d.) rms values of the residual forces and moments at the pelvis for the five subjects. These were obtained when using conventional inverse dynamics analysis of the experimental data. X, Y and Z refer to the anterior, superior and ...
TABLE III
Average (± 1 s.d.) root-mean-squared differences between measured kinematics (obtained by inverse kinematics analysis) and the altered kinematics (obtained by integrating over the accelerations estimated by REA) for five subjects.

The global optimization inverse kinematics routine produced body configurations that were on average within 10.9 (± 1.2) mm across all the measured marker positions. REA generated joint kinematic patterns that exhibited slightly larger differences with measured marker positions, with average differences of 12.9 (± 1.7) mm. Thus, the additional marker discrepancies necessary to produce dynamically consistent motion averaged only 2.0 (± 0.9) mm for a full gait cycle. Similarly, the pelvic translational trajectories were also well reproduced, with REA predicted values generally remaining less than 5 mm of the trajectories estimated via inverse kinematics (Tab. III).

B. Synthetic Data with Noise

Noise polluted synthetic marker data was used to quantitatively asses the impact of REA on the accuracy of kinematic and kinetic measures. White noise in the range of 6 mm to 9 mm added on to the synthetic data produced residual forces and moments on the same order of magnitude as the experimental data. REA eliminated these residuals and also improved the accuracy of joint moment estimations. For a 9 mm sd noise, the moment estimation error was reduced by approximately 59% at the lower back and by 23% at the hip (flexion/extension moments). Errors remained relatively constant at the knee and the ankle (Fig. 5). The level of improvement in joint moments decreased with the magnitude of noise at the back, hip and knee but was relatively independent of noise level at the ankle. (Tab. IV).

FIG. 5
Average joint moment estimation error with a 9 mm SD noise. Values for the flexion/extension of lower back, hip, knee, and ankle joints. These values are computed with respect to noisless synthetic data
TABLE IV
Average % reduction of the joint moment estimation error as a function of noise. Values for the flexion/extension of lower back, hip, knee, and ankle joints. The average inprovement over all joints (not including pelvic degrees of freedom) is given, as ...

The beneficial effect of REA extended to kinematic measures. After the low-pass filtering removed roughly two thirds of the white noise, the algorithm was able to further reduce marker noise by about 9%. This mainly improved the estimate of the position of the pelvis (on which the algorithm has the biggest impact). This position error was lowered by about 19%. The average error of the remaining generalized coordinates (the joint angles) remained nearly constant (Fig. 6).

FIG. 6
The errors in marker positions (top), pelvic translations (middle), and joint angles (bottom) of the noise polluted data are compared to values after data processing with REA. The algorithm was able to reduce marker noise by about 9% and pelvic position ...

C. Synthetic Data with Noise and Offset

Adding constant offsets to the noise polluted data seriously degraded the quality of standard kinematic and kinetic measures, with the marker kinematic errors increasing linearly with the offset magnitude (Fig. 7). However, the inclusion of offset estimation within the REA algorithm corrected for static offsets resulting in kinematic and kinetic errors that were independent of the offset magnitude.

FIG. 7
If REA is estimating constant offset values between optical tracking system and the force plate reference frame, the algorithm can reduce the adverse effect of calibration errors. In this example, a 9 mm sd noise was overlaid with a variable offset that ...

D. Computational Performance

Despite the large (59 dimensional) search space and the complex objective function, the numerical optimization was well behaved, and consistently converged to an optimal set of initial conditions. Computation time for processing a full gait cycle of experimentally recorded data was 100 minutes on a personal computer with a 1.83 GHz Pentium processor. Lowering the number of estimated accelerations to the 6 dof’s of the pelvis substantially reduced the computation time (to less then 4 minutes) while only slightly degrading the ability to reject noise (white bars in Fig. Fig.5,5, ,6,6, and and77).

IV. DISCUSSION

In this article, we have introduced a residual elimination algorithm for computing kinematic and kinetics measures that are consistent with a whole body dynamic model of gait. Joint moments, computed using this approach, exactly reproduce the underlying motion when used to drive a forward simulation. The computed kinematics and kinetics can therefore be used as a starting point for generating a muscle-actuated simulation [12, 15], to identify how individual joint moments contribute to the overall movement [16], or to predict how movement will change in response to interventions that alter the musculoskeletal or nervous system [2].

The residual elimination algorithm was shown to partially reject marker noise and offsets that can arise from measurement noise, soft tissue movement, and calibration errors. As shown by others, improved inverse dynamic results can be obtained by accounting for the redundant nature of whole body kinematics and ground reaction force measures [4, 6, 17, 18]. For example, the vertical ground reaction force provides a measure of the vertical acceleration of the center of mass, which can be mathematically described as a constraint on the generalized accelerations. In the least-squares process, the need for residual forces acting on the base segment is eliminated. While traditional inverse dynamics results are sometimes used to investigate the coordination of movement [16], the presence of residuals means that part of the movement must be attributed to non-physical forces and thus could conceal the influence of actual forces (e.g. muscle) acting on the system.

The use of weighting factors in the least squares allows for uncertainty in specific measurements to be accounted for. If the measurement error covariance matrix is known, this can be used directly as weighting matrix W [5]. Alternatively if the errors are assumed independent, one may use a diagonal matrix with weightings reflecting uncertainty in the specific measures. For example if a quantity is unknown, the according weighting element is set to a very large value. This could be used to compute unknown components of the kinematics or ground reactions. For example when using a simplified skeletal model that represents the entire upper body (trunk, head, and arms) as one rigid segment, one may wish to determine effective lower back angles that account for the entire upper body dynamics. In this case, the low back angles can be computed considering the dynamics of the system simply by setting W large. Another example is experimental setups which only allow the measurement of a subset of the six components of the external reactions, such as when using insole pressure sensors instead of force plates. In cases like these, the least squares formulation serves as a practical tool to estimate the missing components of the external reactions [5, 6].

Our acceleration estimates were integrated forward in time to obtain dynamically consistent kinematic trajectories. Numerical optimization was employed to search for an initial model state that, when integrated, best replicated the recorded motion. It is important to note that using only the six overall equations of motion made this optimization problem computationally feasible to solve. While it would be possible to use the entire set of equations of motion[4, 5], this would also add additional unknowns (joint moments) to the estimation problem which could seriously degrade performance.

To achieve a truly optimal result, the initial positions and velocities of all generalized coordinates have to be adapted simultaneously. However, there are practical limitations. The computational costs of the optimization are substantial, as the entire motion has to be integrated for every evaluation of Eq. (20). It is therefore important that the optimization converges as fast as possible. For this reason, only a small subset of the initial generalized coordinates might become subject to optimization with the remaining ones being set to their experimentally recorded values. We showed in this study that simply adapting the motion of the pelvis reduced the computation time by over an order of magnitude, while only slightly degrading the accuracy of computed kinetics and kinematics. This favorable result is obtained because the ground reaction forces provide information about the motion of the COM which was greatly influenced by the motion of the pelvis, the base segment in our model. In contrast, ground reactions don’t provide much information about the acceleration of segments with small inertia. This is potentially problematic since small acceleration inaccuracies will lead to long term drift. To address this issue, feedback is used to compute a desired set of joint accelerations based on current errors in simulated generalized coordinates and speeds [12]. As a result, this ensures that these coordinates tend to track the desired trajectories.

The optimal determination of initial model states circumvented some of the inherent stability problems associated with simulating human walking. While feedback controllers are useful to ensure proper tracking of joint angles [12, 15], they cannot stabilize the overall ‘inverted pendulum’ characteristics of human gait. Precise identification of the initial conditions is not solving the underlying problem either, but can retard its effect. A short motion, like a cycle of gait, could be well simulated in spite of the inherent instability. This allowed the creation of forward dynamic simulations in an ‘open loop’ framework.

Some limitations arise from the approximations of the model used. These include the assumptions that the body segments are truly rigid and that all movements happen within the given degrees of freedom. We also assumed that all model parameters are precisely known. While other approaches used the redundancies of force and acceleration measurement to improve the assessment of model parameters [1921], our method attributed all differences to errors in the force and acceleration measurements. It is feasible to directly include physical model parameters as variables to be estimated within the optimization, though it would need to be demonstrated that such parameters can be robustly estimated from walking data alone. We also applied measured ground reactions directly to the foot segments rather than trying to model the complexities of the foot-floor interactions. However, it may be preferable to incorporate a contact model into the algorithm [15] when one desires to use the simulation to predict how movement changes in response to perturbations. Finally, we considered the ability of our algorithm to reject both systematic (bias) and white noise. Future studies should consider the use of continuous noise [20] or simple soft tissue models, which may better represent the marker displacements arising from soft tissue motion, which remains the major source of error in traditional motion analysis.

The proposed algorithm can be extended to accommodate a subset of prescribed motions. Examples might include experiments where the subject is attached to objects that are fixed in space or follow given trajectories (like handles, guardrails, or exercise equipment). If the prescribed motion is expressed directly in terms of generalized coordinates, the prescribed accelerations q¨^ can simply replace the measured ones q. Only the remaining generalized accelerations are adapted to achieve dynamic balance. However in many cases, the accelerations are probably defined for individual segments and therefore prescribed in Cartesian coordinates. In these case, the prescribed accelerations r¨^ and ω.^, don’t map directly into generalized accelerations q¨^, but instead form complex geometric constraints that can be incorporated into the optimization problem (see Appendix).

In summary, we have introduced and demonstrated a residual elimination algorithm (REA) that can be used to generate joint moment-actuated dynamic simulations of subject-specific gait dynamics. The algorithm was shown to eliminate the need for residual forces while also reducing the detrimental effects of measurement errors on joint kinematics and kinetics. REA-generated simulations can therefore be used as a basis for investigating the co-ordination of movement, and for predicting changes in movement due to injury, impairment or surgical intervention.

FIG. 4
Shown are the pelvic coordinates and joint angles estimated using inverse kinematics (solid lines) and the residual elimination algorithm (dashed lines) for one of the experimental data sets. Very good agreement is seen with average differences generally ...

V. ACKNOWLEDGEMENTS

This work was supported by Shriner’s Hospitals for Children, NIH AG24276, and NIH GM72970. The authors gratefully acknowledge Dr. Bryan Heiderscheit, Amy Silder and Ben Whittington for providing the experimental data used in this study.

VI. APPENDIX

REA can be extended to accommodate a subset of prescribed motions expressed in cartesian coordinates. To include these constraints into the formulation derived above, we refer back to Eq. (7). Here the Cartesian constraints can be introduced mathematically by expanding the matrix Acart. One row will be added for each constraint. Each of these rows contains a single 1 which maps a Cartesian acceleration to its prescribed value (denoted by a^ ) on the right hand side of the linear equation:

[[Acart(r)]0001000010]=[FMr¨ω˙]=[[fcart(r,ω)]r¨^iω.^i]
(21)

A^cart(r)[FMr¨ω˙]=f^cart(r,ω,r¨^,ω˙^)
(22)

This formulation has the advantage that it keeps a full vector of external reactions and Cartesian accelerations. It guarantees that the remaining problem formulation, as discussed earlier, can be kept as stated. It will inherently take care of the geometric constraints that result when the prescribed Cartesian accelerations are mapped into generalized coordinates. We only have to replace the matrix Acart and the function fcart by their extended counterparts and base the residual formulation on these new matrices. However, it is important to note that every constraint effectively removes one dimension of the least-squares solution space. If more constraints than generalized coordinates are introduced A will contain more rows than columns and Eq. (15) becomes overdetermined. However, even when dealing with fewer constraints, singularities in J map directly into A, and thus even a single constraint can cause (15) to become unsolvable for certain model configurations. For this reason, one would have to exercise caution when introducing motion constraints in a model that traverses near to singular configurations.

References

[1] Zajac FE, Neptune RR, Kautz SA. Biomechanics and muscle coordination of human walking part ii, lessons from dynamical simulations and clinical implications. Gait and Posture. 2003;vol. 14:1–17. [PubMed]
[2] Delp SL, Anderson FC, Arnold AS, Load P, Habib A, Chand J, Guendelman E, Thelen DG. Opensim: Open-source software to create and analyze dynamic simulations of movement. IEEE Transactions on Biomedical Engineering. 2007 in press. [PubMed]
[3] Winter D. Biomechanics and Motor Control of Human Movement. Wiley-Interscience; 1990. ch. 4 - Kinetics: Forces and Moments of Force; pp. 75–102.
[4] Cahout V, Luc M, David A. Static optimal estimation of joint accelerations for inverse dynamics problem solution. Journal of Biomechanics. 2002;vol. 35:1507–1513. [PubMed]
[5] Kuo A. A least-squares estimation approach to improving the precision of inverse dynamics computations. Transactions of the ASME. 1998;vol. 120:148–159. [PubMed]
[6] van den Bogert AJ, Su A. A weighted least squares method for inverse dynamic analysis. Computer Methods in Biomechanics and Biomedical Engineering. 2007 in press. [PubMed]
[7] Chao E, Rim K. Application of optimization principles in determining the applied moments in human leg joints during gait. Journal of Biomechanics. 1973;vol. 6:497–510. [PubMed]
[8] Anderson F, Pandy M. Dynamic optimization of human walking. Journal of Biomechanical Engineering. 2001;vol. 135(no. 5):381–390. [PubMed]
[9] Neptune R, Kautz S, Zajec F. Contributions of the individual ankle plantar flexors to support, forward progression and swing initiation during walking. Journal of Biomechanics. 2001;vol. 34:1387–1398. [PubMed]
[10] Greenwood D. Principles of Dynamics. Prentice Hall; 1988. ch. 8 - Dynamics of a Rigid Body; pp. 389–468.
[11] Yamaguchi . Dynamic Modeling of Musculoskeletal Motion. Kluwer Academic Publishers; 2002. ch. 7.5.3 - Properties of the Pseudoinverse; pp. 225–226.
[12] Thelen D, Anderson F. Using computed muscle control to generate forward dynamic simulations of human walking from experimental data. Journal of Biomechanics. 2006;vol. 39:1107–1115. [PubMed]
[13] McConville JT, Churchill TD, Kaleps I, Clauser CE, Cuzzi J. tech. rep.
[14] Lu T, J. J. O. Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints. Journal of Biomechanics. 1999;vol. 32:129–134. [PubMed]
[15] Seth A, Pandy M. A neuromusculoskeletal tracking method for estimating individual muscle forces in human movement. Journal of Biomechanics. 2007;vol. 40:356–366. [PubMed]
[16] Siegel KL, Kepple TM, Stanhope SJ. Using induced accelerations to understand knee stability during gait of individuals with muscle weakness. Gait and Posture. 2006;vol. 23:435–440. [PubMed]
[17] Cappozzo A, Leo T, Pedotti A. A general computation method for the analysis of human locomotion. Journal of Biomechanics. 1975;vol. 8:307–320. [PubMed]
[18] Challis JH, Kerwin DG. Quantification of the uncertainties in resultant joint moments computed in a dynamic activity. Journal of Sports Sciences. 1996;vol. 14:219–231. [PubMed]
[19] Vaughan C, Andrews J, Hay J. Selection of body segment parameters by optimization methods. ASME Journal of Biomechanical Engineering. 1982;vol. 104:38–44. [PubMed]
[20] Reinbolt JA, Schutte BJ, Fregly BJ, Koh BI, George AD, Mitchell KH. Determination of patient-specific multi-joint kinematic models through two-level optimization. Journal of Biomechanics. 2005;vol. 38:621–626. [PubMed]
[21] Reinbolt JA, Haftka RT, Chmielewski TL, Fregly BJ. Are patient-specific joint and inertial parameters necessary for accurate inverse dynamics analyses of gait? IEEE Transactions on Biomedical Engineering. 2007;vol. 54:782–793. [PMC free article] [PubMed]