|Home | About | Journals | Submit | Contact Us | Français|
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.
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.
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.
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 . 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 . 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 . 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 [4–6]. 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  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.
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 .
Similarly, the equations for rotational motion about the origin of the coordinate system can be stated as:
and are external forces and moments respectively, is the location of a point along the external reaction force , points to the center of gravity, is the angular velocity, and mj and Ij are the mass and the inertia dyadic of segment j. Gravity is denoted by the vector . 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.
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 and . These generalized coordinates (that mostly correspond to joint angles) can be mapped to the Cartesian coordinates by a nonlinear transformation t.
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 , it would be more appropriate to state )
Differentiating this expression leads to:
Where Acart is a matrix of dimension 6 × (n · 6 + N · 6) and 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 . 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.
which can be used to state the equations of motion as a linear function of the generalized coordinates:
where the Cartesian coordinates and speeds have been calculated from the generalized co-ordinates using
The matrix A is of the dimension 6×(N · 6 + nq), where nq denotes the number of generalized coordinates q.
As a result of model and measurement errors, the experimental recordings of ground reactions and accelerations will generally not satisfy Eq. (9). Small variations (δ, δ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:
where , and 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 () and moments (). The desired accelerations include adjustments to the experimentally derived accelerations () that account for current errors between the experimental generalized speeds () and coordinates () and the corresponding simulated values:
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.
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:
that will minimize the Euclidean norm . 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 , where only the six pelvic degrees of freedom were altered.
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).
The agreement of recorded and calculated motion can be assessed by the squared and weighted distances between calculated and experimentally derived generalized coordinates:
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:
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:
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).
For all analyses, a 3-dimensional full body model consisting of 12 independent segments was used. The generalized coordinates 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:
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).
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.
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).
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).
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).
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.
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).
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 , or to predict how movement will change in response to interventions that alter the musculoskeletal or nervous system .
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 , 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 . 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 . 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 [19–21], 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  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  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 can simply replace the measured ones ’. 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 and , don’t map directly into generalized accelerations , 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.
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.
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 ) on the right hand side of the linear equation:
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.