Improving Inverse Dynamics Accuracy in a Planar Walking Model Based on Stable Reference Point

Physiologically and biomechanically, the human body represents a complicated system with an abundance of degrees of freedom (DOF). When developing mathematical representations of the body, a researcher has to decide on how many of those DOF to include in the model. Though accuracy can be enhanced at the cost of complexity by including more DOF, their necessity must be rigorously examined. In this study a planar seven-segment human body walking model with single DOF joints was developed. A reference point was added to the model to track the body’s global position while moving. Due to the kinematic instability of the pelvis, the top of the head was selected as the reference point, which also assimilates the vestibular sensor position. Inverse dynamics methods were used to formulate and solve the equations of motion based on Newton-Euler formulae. The torques and ground reaction forces generated by the planar model during a regular gait cycle were compared with similar results from a more complex three-dimensional OpenSim model with muscles, which resulted in correlation errors in the range of 0.9–0.98. The close comparison between the two torque outputs supports the use of planar models in gait studies.


Introduction
In inverse dynamics (ID) problems, the objective is to determine the joint torques and joint forces during movement [1,2].The torque calculations are based on availability of motion measurements, ground reaction forces, and the subject's anthropometrics.ID have been used in clinical applications to model the human body in attempts to better understand pathological gait.Based on the increasing use of mathematical models in human gait research, it has become more regularly utilized and a valuable tool.Generally, a bipedal locomotion system consists of several rigid bodies that are interconnected through articulated and actuated joints [3][4][5].The human body may be modeled as simply as three point masses connected with massless links [6,7], or it could be a complex construction with multiple muscles spanning joints with multiple DOF that bears closer resemblance to the actual human body [3,8,9].Further, a bipedal model can be characterized as an inverted pendulum system with constrained motion, which represents the effects of forward and backward impacts of the swinging limb with the ground [10][11][12].Knowledge about the internal structure of the locomotion system is essential in any study concerning gait dynamics [13].
Two different methodologies can be considered to model ambulatory motion of the human body.Models are based on either Newton-Euler method with the reference point at the hip [14,15] or employing Lagrange formulae [16][17][18].Blajer and Czaplicki utilized a multibody planar model with 9 DOF to study the human body performing a jump on a trampoline [14].Their model consisted of a one limb lower body having 3 DOF, an upper body of 3 DOF, and a reference point at the hip having 3 DOF.Kljuno and William II considered a 10-DOF planar walking model, distributed as follows: 6 degrees for the lower body and 4 degrees for the upper body [15].The latter used the hip as the reference point consisting of 2 translational and 1 rotational DOF for the trunk.Plestan et al. employed Lagrange formulae to build a 5 segment 7 DOF humanoid robot with 1 rotational DOF for each segment and 2 translations at the hip [16].Plestan et al. 's planar model consists of an upper body represented as one segment, two thighs, and two shanks with no feet.A similar model structure has been depicted by Imani et al. [17], who kept the model balanced by insuring an equal distance from the swing and the stance limbs to the body's center of mass.Hurmuzlu et al. [18] considered a 7 segment bipedal planar model with the reference point at the center of mass of the upper body.Hurmuzlu et al. also identified that, in modeling, the reference point can be considered at any place of the human body.However, considering the hip (and the pelvis in general) as a reference position generates error in modeling because in walking the pelvis undergoes rotations.
The pelvis in the human body is considered an unstable part of the body [20].Further, Van Leeuwen et al. [21] reported a pelvic rotation during human walk by approximately four degrees in the transverse plane and even more rotation in the frontal plane, that is, five degrees.Because the pelvis is unsupported, this rotation likely results from the reactive forces generated by moving the limbs and hitting the ground.Therefore, this places the two hip joints at different positions while walking.Hence, considering one of the hip joints as a reference point for modeling purposes will generate different results if the other joint had been considered.In addition, this consideration will affect measuring the joint angles trajectories.Despite instability of the pelvis and its rotation during walking, it has been considered as a reference point by several researchers.The errors resulting by the nonstationary nature of the pelvis will impact the calculations of the joint angles as well.These errors will be magnified especially if joint mechanics have been calculated from the external markers directly without considering the joint centers.Hence, essential elements in ID problems are the joint angle trajectories which require sufficient accuracy [22].
In this paper, a different reference point has been suggested in order to increase accuracy of a planar human body model.Modeling and validating a seven-segment structure had be undertaken based on the ID approach and considering a reference point at the top of the head.This particular reference position was chosen for two reasons: first, to avoid using an unstable part of human body as a reference point; second, this position physiologically signifies the location of human body acceleration sensors, that is, the vestibular system [23,24].In addition, mathematically, measuring this position from marker data is easier as it will mainly depend on the subject's height; further, its motion can be considered simply as sinusoidal during walking.In order to test validity of the model actual motion data of human walking was used as a comparison.Accurate measurements of the input trajectories were based on using the regression method developed by Davis III et al. [25], which has the advantage of simplicity and the results reported are approximate to a complex 3D human body model.The model developed in this study is cross-validated with outputs from the OpenSim software, which is able to model, simulate, and analyze 3D musculoskeletal movements [26].

Materials and Methods
An average North American male subject performed selfselected normal speed walking for 8 trials.Tracks of 16 passive optical markers were used to capture the body motion configured in Helen-Hayes marker set.Ground reaction forces were measured using 2 Bertec (Columbus, OH) force platforms.

Biped Model.
The planar bipedal model undertaken in this study has 7 rotational DOF represented as revolute joints, as in Figure 1.These include two ankle angles, two knee angles, two hip angles, and one rotational degree for the headarms-trunk (HAT) segment.Two translational degrees are also included, which represent the reference position (, ) of the head.The upper body segments have been accumulated into one segment because they all act together to keep the body balanced [27].Furthermore, no significant bending occurs during gait at the upper part of the body [15].Human hip and ankle have the ability to rotate in three dimensions, which justifies representing them as a ball-socket joint [2,28].However, for simplicity, these joints will be considered as revolute joints [8,29]; hence, analysis of the model will only consider motion in the sagittal plane.

Human Walking.
Due to the complexity of modeling the human body, specifically for walking, a brief description of the terms and the sequence forms of this movement is described here.Human walking involves locomotion using both legs alternately for balance and propulsion [27].Two different events separate single support from double support during the normal (i.e., healthy) gait cycle, as in Figure 2: heel strike (HS) is characterized as the instant when the swing leg contacts the ground (i.e., with the heel), while toe-off (TO) is the instant when the ipsilateral toe leaves the ground [30] at the end of the support phase which also signifies the initiation of the swing phase [30].The pathway of the body's center of mass and any other point on the upper body during a gait cycle takes a periodic form due to limbs swing, which makes the torso drop slightly as shown for C7 and SACR markers, as in Figure 3.It is obvious that the upper body has a periodic cycle which will be easy to generate through fitting equations.Although, the waveform for the lower body repeats itself for each step (half stride), it still shows nonlinearities and randomness, as in Figure 4. Altitudes of the heel and the toe of right limb are shown in Figure 4, and HS (dashed blue line) and TO (dashed red line) events are clear when their altitude is at minimum.The gap between HS and TO events is the period of the double support and between TO and HS is the single support period and is repeated for each limb alternatively [20].

Equations of Motion.
Segmental inertia was obtained by regression formulae that are based on a number of anthropometric measurements, including body mass, inertia, and segment dimensions [31], as in Figure 5. Equations of motion can be generated by using either method: the Newton-Euler method or the Lagrange form [2].
To obtain equations of motion using the Newton-Euler method, it is required to determine the segments center of mass (COM) and the joint positions from top down.Hence, there will be a total of 21 equations, p, concerning translation and rotation of each segment using nine generalized variables, q: Substituting q into p yields the system in (2) and can be differentiated with time to find velocity (3) and accelerations (4): where g is an arbitrary function that contains the kinematic constraints ( 5)-( 7); D = g/q is a 21×9 dimensional matrix; and E = D q can be determined using Implicit function theory [32].Using the first segment as an example to clarify the way to find the matrices in (3)-( 4), one has Forces and moments of inertia of the lower body segments will be used to determine the moments at the joints to generate an equation for each joint, as shown in (10); calculations will be from bottom up, for example, measuring ankle joint kinetics from ground reaction forces: A total of 6 kinetic equations (for each of the forces and moments) will be developed to be distributed among the seven segments, as shown in (11).For example, the HAT segment is affected by both hips' torques so this segment equation will be generated from the addition of the two hip torque equations.Therefore, depending on these six torque equations, seven equations of motion have been formed which are functions of the 9 generalized variables, as shown in (11).This conversion is not required using the Lagrangian method because the equations depend on the energy of each segment which will generate the seven equations of motions [2].Consider Writing equations of motion in a compact form to be a function of generalized variables, q, will reduce the matrices size and consequently will reduce the time required to compute inverse matrices.A process of converting the system of equations to matrices will depend on the number of segments and the generalized variables [33].Compact matrix dynamic equations will have the form as described.
For the current model, the resultant system will have a mass matrix M of dimension 7 × 9, and the same dimension for the Coriolis matrix C. The gravitational vector G is 7 × 1, while the distribution torques matrix Q is of dimension 7 × 6 dimension.The constraints matrix  contains 2 equations ( 12)-(13) when the system is on single support (SS), and contains 4 equations ( 12)-( 15) when it is on double support (DS), as shown in Figure 5. Lagrange coefficients  correspond to ground reaction forces, as follows:  1 ,  1 for SS and  1 ,  1 ,  2 ,  2 for DS.Consider 2.4.The Choice of Reference Point.Global position, or reference point, of the humanoid model can be located at any point on the bipedal model [18].The interconnection of upper and lower body parts, that is, the hip joint, has been previously used by several researchers [14,34].However, during walking, the contralateral hip joints have slightly different sagittal plane locations throughout the gait cycle; thus, considering either one of them as a reference develops different results than considering the other joint.Furthermore, hip joint location will affect other joint angle trajectories which contribute to the system stability [18].Therefore, this section involves the developed error at the hip joints when considering either of the two hip joints as a reference point.Hip joint center (HJC) is located at the center of the acetabulum hemisphere at which the femur head meets the pelvis [35], as shown in Figure 6.Locating the HJC precisely is important [36], and for accuracy purposes the hip joint center will be considered.Estimating the HJC has been developed previously based on the anatomical frame of the pelvis and its origin [25].Formulae used to calculate the HJC can be found in [25] and will not be repeated here.Further, only two-dimensional results will be considered by neglecting the mediolateral location because only sagittal plane of the human body is interested here.
HJC of an average North American male has been estimated doing 8 trials of walking.Anatomical positions of the joint centers of both hips have been estimated using the proposed formulae by Davis III et al. [25].After determining the anatomical positions of the HJCs, it is required to convert them to global position, as described by Winter [31].

Simulation Results
To show the drawback of picking either joint as a reference point, the difference between both HJCs was taken for all trails for one walking cycle, as shown in Figure 7. Horizontally, in the sagittal plane, the maximum difference obtained was 9.5 mm, which is a high figure concerning center of joint [35], as in Figure 7(a).This error will increase if the marker on the tissue has been considered as a reference point instead, because hip rotation is higher on the pelvic tissue rather than the HJC based on pelvic circumference.In the transverse plane, in Figure 7(b), the difference was higher and was approximately 11.5 mm.In addition, the vertical position difference of the HJCs expose randomness in measurements that will be an additional drawback.As a result, picking either joint will generate a Euclidean distance to the other joint of approximately 14.9 mm.Also, this error will affect calculating the joint angles because precise joint angles depend on the joint centers [31], hence affecting accuracy of measuring the joint angle trajectories, which is not recommended [22].
According to Newton's third law, when the heel strikes the ground a reaction is rebuilt at the hip joint which pushes it posteriorly from its position.Therefore, this increases the difference between the two HJCs in the duration from HS to TO, as in Figure 7.Moreover, largest differences occur in the DS phase, for both the sagittal and the transverse planes.To further illustrate the effect of pelvic rotation while walking, Table 1 was prepared which contains average and standard deviation values of minima and maxima of all the frames of all trials.Largest difference has been recorded in the DS phase, which is approximately, (18.4 mm) horizontally and (11.7 mm) vertically.It is worth mentioning that these differences will be higher if minima and maxima of the differences between the HJCs have been considered instead of the average values, as in Table 1.Therefore, selecting the reference position at a more stable position, the top of the head, for example, will reduce these errors.In addition, there will be just one point to select instead of picking from two different joints which was the main cause of these differences especially in the horizontal axis.Besides, generating a formula for the top-head reference point is simple due to its sinusoidal periodic behavior.Furthermore, it has similarity to human body physiological position of the balance sensor.For these reasons, the top of the head is a preferred reference point when considering modeling a walking human body.
Before applying the trajectories to the model, we would like to mention the validation method that has been depicted in this paper.Lately, an open source software called Open-Sim was developed by the National Center for Simulation in Rehabilitation Research (NCSRR), Stanford, CA, USA.This software is able to model, simulate, and analyze 3D musculoskeletal movements [26], which uses 12 DOF for the lower body.OpenSim software has been used to model and simulate walking models by several researchers [37][38][39][40].OpenSim can do inverse dynamics to find joint kinetics based on providing a marker set and the ground reaction forces.Therefore, its outputs will be considered by comparing it with the proposed human body model outputs in order to validate the model.Ground reaction forces (GRF) of both limbs are sketched in Figure 8, where first letter in the legend indicates the limb and  and  represent horizontal and vertical components, respectively.Number of the trajectories used is 9 and consists of two translational coordinates representing the body reference point and seven segments rotations, as in  Applying the joints trajectories and the GRFs on the model developed in Section 2.1 using ID equations proposed in Section 2.3 will generate the net joint torques, as in Figure 10.The figure shows the joint torques outputs of both the proposed model (blue) and the OpenSim software (red) for both limbs (R, L).Close results have been obtained by comparing both outputs which reveals effectiveness of the developed humanoid model.The improved accuracy of simulation outputs results from several factors.First, picking the reference point at a stable position increased the accuracy of the starting point of equations of motion.Second, accurate joint trajectories were computed, which were based on determining the joint angles from center of joints.Third, picking two different joints for the hips gave the model the close look of a 3D model.Finally, differentiating the signals requires accuracy and they have to be filtered [31].Therefore, to accomplish both requirements (differentiation and filtering), the method of digital low-pass differentiation has been used in order to find first and second differentiations [41].

Conclusions
In this study, a simplified bipedal human locomotion model was developed and employed in an inverse dynamic approach to obtain estimates of joint torques.At the same time, accuracy of joint torques was improved by avoiding known sources of error and providing accurate input trajectories.Simplification was exacted by including only those joints that have main effect on the human gait and approximating the upper body with a single segment.Precise measurement of joint trajectories becomes essential in human body inverse dynamics because they are the main input to the model.Known sources of error in modeling were identified and methods of avoiding them were adopted.These included avoiding the unstable pelvis as a reference point and considering the center of joints to determine the joint angles.In addition, filtering the piece-wise function form developed via numerical differentiations of the joint angles was undertaken before utilizing those signals in the model.
Studies have shown existence of pelvic rotation during walking in the frontal and the transverse planes for approximately 5 ∘ and 4 ∘ , respectively [21].Because the frontal plane was not considered in this study, the modeling error is expected to increase if a 3D model had been considered.This study shows that pelvic rotations can develop errors between the two hip joint centers of approximately 9.5 mm and 11.5 mm in the frontal and transverse planes, respectively.In addition, maximum differences obtained between both limbs' HJCs for a single walking stride were 18.7 mm and 27.5 mm for the vertical and the horizontal planes, respectively.These errors were sufficient to avoid using the hip as a global position indicator for the humanoid model and simulations results support that hypothesis.

Figure 5 :
Figure 5: Forming constraints on single support (a) and double support (b).

Figure 6 :
Figure 6: Anatomical structure of the pelvis showing HJC location.

Figure 7 :Figure 8 :
Figure 7: HJC difference between left and right limbs for horizontal axis (a) and vertical (b) axis.

Figure 9 :
Figure 9: Trajectories used as input to the humanoid model.

Figure 10 :
Figure 10: Joint torques output of the developed model (blue) and OpenSim outputs (red): (a) for right limb, and (b) for left limb.

Figure 9 .
Figure 9. Trajectories considered are for two strides and the joint angles are of absolute horizontal type.Applying the joints trajectories and the GRFs on the model developed in Section 2.1 using ID equations proposed in Section 2.3 will generate the net joint torques, as in Figure10.The figure shows the joint torques outputs of both the proposed model (blue) and the OpenSim software (red) for both limbs (R, L).Close results have been obtained by comparing both outputs which reveals effectiveness of the developed humanoid model.The improved accuracy of simulation outputs results from several factors.First, picking the reference point at a stable position increased the accuracy of the starting point of equations of motion.Second, accurate joint trajectories were computed, which were based on determining the joint angles from center of joints.Third, picking two different joints for the hips gave the model the close look of a 3D model.Finally, differentiating the signals requires accuracy and they have to be filtered[31].Therefore, to accomplish both requirements (differentiation and filtering), the method of digital low-pass differentiation has been used in order to find first and second differentiations[41].

Table 1 :
Average (avg) and standard deviation (std) values of minima and maxima of the HJC differences for both single support (SS) and double support (DS) phases.