Modeling and computational issues in the inverse dynamics simulation of triple jump

The triple jump is a demanding field event consisting of an approach run, and then followed by a hop, a bound, and a jump. The three consecutive takeoffs are executed at high speed, during which a jumper must absorb extremely large impact forces. The purpose of this paper is to develop an effective formulation for the inverse dynamics simulation of all the jump phases separately. A planar model of the jumper is used, composed of 14 rigid segments connected by 13 hinge joints, and actuated by muscle forces in the lower limbs and resultant muscle torques in the upper body joints. The equations of motion of the model are obtained using a projective technique, allowing for effective assessment of the ground reactions as well as muscle forces and joint reaction forces in the lower limbs. Some numerical results of the inverse dynamics simulation of a triple jump are reported.

"cycling" movement-the athlete takes off from one leg, cycles it through, and completes landing on the same leg. The immediately followed step finishes with landing on the other leg, and the final jump from the non-take-off leg is very similar to the long jump. During the takeoffs and landings on the runway, all executed at high horizontal speed, the jumper must tolerate extremely large impact/impulsive forces in the leg contacting the ground, and all the distinct phases must be learned and practiced to combine them in one successful long distance performance. The triple jump is therefore recognized as one of the most technically and physically demanding events in track and field, also from the biomechanical point of view [1][2][3][4][5][6].
Earlier biomechanical studies of the triple jump focused mainly on qualitative comparisons of the triple jump techniques for individual athletes [2,6,7] or some more specific analyses like the optimum phase ratio [3,6,8], functions of arm swing motions [2,9], tradeoffs between velocity components during triple jumping [10], and biomechanical loading [4,11]. The present paper attempts to supplement these studies with the inverse dynamics analysis of the triple jump, based on musculoskeletal modeling and computer simulations, and aimed at quantitative evaluation of the lower-limb muscle forces and joint reaction forces during the movement. The insight into how the muscles interact to produce the motion, and the assessment of the involved internal loads may be of importance for a better understanding of both the triple jump techniques and possible injury mechanisms. While inverse dynamics simulation is a widespread and intensively used as a noninvasive method in biomechanics, especially in the analysis of lever walking (see, e.g., [12][13][14][15] and the literature contained therein), its applications to triple jump still remain to be designed.
Depending on the expected accuracy and scopes of the analysis, the human body can be modeled using different methodologies and ranges of approximation of its very complex mechanical structure and actuation; see, e.g., [16][17][18][19]. The present formulation is designed for the inverse dynamics simulation of the triple jump in all the movement phases, irrespective of the fact whether the jumper flies or contacts the ground with one of his legs. Under the assumption that the motion of all the body segments is performed with some accuracy in the sagittal plane, a planar human body model is developed, which extends and improves the model used previously in [20]. It is built as a simple kinematic structure consisting of rigid segments connected by ideal hinge joints, and branching from the head segment in the open chain linkages. Two models of control are used: the torque driven one that represents the resultant muscle action at the joints, and the other model with mixed set of muscle forces at the lower limb and the resultant muscle torques at the other body joints. The torque driven model results in the determinate inverse dynamics problem, from which, using the prescribed motion characteristics, the involved muscular torques in all the joints as well as the reactions from the ground can explicitly be determined. The force-control model allows then for a more detailed analysis of internal loads in the lower limbs, and results in the indeterminate inverse dynamics problem in which the previously obtained muscular joint torques in the lower limb joints are distributed into the individual muscle efforts by means of static optimization. A convenient decomposition-based formulation for solving of this muscle force sharing problem is proposed. Then, using the assessed muscle forces variations, an effective computational scheme is offered for the determination of reaction forces selectively in the lower limb joints. Some numerical results of the inverse dynamics simulation of a sample triple jump are finally reported.

Basic assumptions
The developed model of the jumper is designed as a planar kinematic structure composed of b = 14 rigid segments connected by k = 13 ideal hinge joints, and branching from the head segment (point T ) into the open chain linkages (Fig. 1). Considering the model as a flier, its number of degrees of freedom is r = 3 + k = 16. The flier can then contact the ground with one of his feet, which yields external (ground) reaction forces on this particular foot.
The dynamic equations of the jumper can conveniently be introduced in n = 3b = 42 absolute coordinates p = [x C 1 y C 1 θ 1 . . . x C b y C b θ b ] T that specify locations of the segment mass centers and the segment angular orientations with respect to the inertial frame XY , where Mp = f g are the combined dynamic equations of unconstrained segments, with the generalized mass matrix M = diag(m 1 , m 1 , J C 1 , . . . , m b , m b , J C b ) and the generalized force vector due to the gravitational forces f g = [0 −m 1 g 0 . . . 0 −m b g 0] T , and where m i and J C i are the masses and the central moments of inertia of the segments, and g is the gravity acceleration. The generalized force vectors f r and f u are due to the external ground reactions and muscle force enforcement (actuation), which will be discussed later. Finally, the generalized force vector f c = −C T λ results from l = 2k = 26 constraints (p) = 0 on the segments due to their connections in the joints, C = ∂ /∂p is the l × n constraint matrix, and λ = [λ 1 . . . λ l ] T are the corresponding constraint reactions (Lagrange multipliers). As it will be seen further, neither the constraint equations (p) = 0 nor constraint Jacobian C need to be introduced explicitly in the followed formulation. The modeling effort is therefore limited to the formulation of f r and f u .

Reactions from the ground
The triple jump is composed of consecutive flying phases, when there is no contact with the ground, and single-support phases, when the jumper touches the runway with one of his feet (the final, not analyzed in this study landing in the sand pit is usually on two feet). The external reaction exerted on a supporting foot can be reduced to point P (Fig. 2), and modeled by means of l r = 3 components λ r = [R x R y M P ] T irrespective of the possible contact scenarios, where R x and R y are the X and Y components of the reaction force, and M P is the reaction force moment about point P of the foot segment. During the flying phases λ r is expected to vanish, which may be a criterion of validity of the developed dynamical model and accuracy of the recorded kinematic characteristics.
The generalized force vector f r introduced in Eq. (1) can symbolically be represented as where the n × l r (42 × 3) matrix A r of distribution of λ r in p directions has nonzero entries only in the rows corresponding to the supporting foot absolute coordinates. The explicit formulation of A r is a standard task, however, two versions of A r must be formulated since the supporting foot can either be the right or left foot, and one must switch between the formulations accordingly.

Models of actuation
Two models of actuation are considered in this study. The deterministic model of actuation involves k = 13 torques u τ ≡ τ = [τ 1 . . . τ 13 ] T that represent the resultant muscle action at the respective joints ( Fig. 1). In the nondeterministic model of actuation, the k = 3 torques in the supporting leg are replaced with the action of m = 9 muscle forces seen in Fig. 3, and the respective control parameters are muscle stresses σ = [σ 1 . . . σ 9 ] T , σ j = F j /A j , where F j is the j th muscle force and A j is its physiological cross sectional area. Denoting symbolically u τ = [τ T τ T ] T , where τ are the k torques in the supporting leg, and τ are the k = k − k = 10 torques in the other joints (including the nonsupporting leg joints), the nondeterministic model of actuation is related to u τ σ = [τ T σ T ] T , which is a mixed set of resultant muscle torques and the muscle stresses. Applying this hybrid model of actuation a more detailed analysis of internal loads in the supporting leg is granted, while retrieving the dynamic interaction between the upper body and the locomotion apparatus; see also [20] where a similar methodology was used.
The generalized force vector f u in Eq. (1), respectively for the deterministic and nondeterministic models of actuation, is modeled as where the control distribution matrices B τ and B τ σ are of dimensions n × k and n × (k + m), respectively. The n × k matrix B τ is the same in both formulations, and the dimensions of B τ and B σ are, respectively, n × k and n × m. The deterministic and nondeterministic models of actuation are not equivalent, i.e., B τ u τ = B τ σ u τ σ , and more strictly B τ τ = B σ σ . Evidently, the muscle forces in the nondeterministic model must result in the same control torques at the respective joints as those (1) by f c = −C T (p)λ, and this contribution is not taken into account when the deterministic model of actuation is used. The determination of joint reactions should therefore include the nondeterministic model of actuation, which will be discussed further.
The formulation of B τ in the deterministic model is obvious. It is a constant sparse matrix with only two nonzero entries in each column, one of them equal to 1 and the other equal to −1; see [20] for details. The formulation of is much more challenging. Each muscle in the lower limbs must be individually modeled as concerns its action on the skeleton segments. Of importance is an accurate identification of the anatomical origin and insertion points of the muscles, and then the muscle paths relative to the skeleton to determine the actual muscle moment arms about the joints [16,[21][22][23]. Implicit in the model is that, in the neighborhood of the joints, the tendons of some muscles may wrap on joint capsules, bones, and other muscle bellies (a certain moment arm with respect to the joint is reached irrespective of the relative configuration of the adjacent segments), or may be constrained by retinacula. Consequently, the effective origin and insertion points can be estimated, in which the muscle forces are applied to the segments, and the effective musculotendon paths are then defined as lines connecting the effective attachment points. All these issues allow for an explicit formulation of B σ (p) for the supporting leg in a way similar to that described in [23] for the upper limb.

Independent and open-constraint coordinates
The initial dynamic equations of the jumper model, introduced in Eq. (1), are formulated in For the aims stated below, it is desirable to project the equations into the directions of r = 16 independent coordinates defined here as q = [x T y T θ 1 . . . θ b ] T , where x T and y T are the coordinates of point T at the top of the head segment, and θ i are the segment orientation angles used in p. The relationship p = g(q) expresses then the joint constraint equations given explicitly [24][25][26], so that the constraint equations given implicitly, z = (p) = 0, are satisfied identically when using p = g(q), i.e., (g(q)) ≡ 0. As shown in [27], the explicit form of joint constraint equations can also be generalized to the augmented explicit form, which involves the open-constraint coordinates z = [z 1 . . . z l ] T that describe the prohibited relative motions in the system An advantageous feature of using open-constraint coordinates is that they can be introduced only in those joints in which the reaction forces are to be determined; see [27] for the theoretical background and [20,23,28] where the approach was applied to biomechanical modeling.
Since, by principle, z = 0, the augmented and traditional formulations of explicit constraint equations, p = g(q, z ) and p = g(q), are virtually equivalent, and the dependence on z is introduced here only to grasp the directions ofż , which are also the directions of constraint reactions λ (Fig. 4). Time differentiation of Eq. (4) leads tȯ where the matrices D and E are of dimensions n × r and n × l , respectively, and E is constant (and simple) for the case of hinge joints in planar systems; see [20,23,28] for more details. Evidently, the matrix D introduced in Eq. (5) can also be derived from the traditional formulation p = g(q) asṗ = (∂g/∂q)q ≡ D(q)q, and it is an orthogonal complement matrix to the constraint matrix C [24][25][26], i.e., Then, as shown in [20,23], related to the symbolic denotation where 0 is the l × l null matrix, and I is the l × l identity matrix, l = l + l .

Projection of the initial dynamic equations
The initial dynamic equations introduced in Eq. (1), after applying the modeling issues of Sects. 2.2 and 2.3, become where the vector of generalized actuating force is either f u = B τ u τ or f u = B τ σ u τ σ , respectively, for the deterministic and nondeterministic models of actuation. With the use of D (related toq) and E (related toż and λ ), the initial dynamic equations can be projected into the directions of q and z /λ , respectively. The projection formula is [27] D Since, according to Eq. (6), D T C T = 0, the dynamic equations projected in the directions of q become It should be emphasized that in the above projections one uses the initial dynamic equations in absolute coordinates p, i.e., M, f g , A r , and f u = Bu (either B τ u τ or B τ σ u τ σ ). The involvement of independent coordinates q and open-constraint coordinates z is required only to formulate p = g(q, z ), and then D and E . Also, while in general D depends on q, for the case at hand D depends only on the segment orientation angles θ i 's (sinuses and cosines of the angles), which are included in both p and q. That is why one can also write D(p), and, as stated above, E is a constant matrix. Therefore, both Eq. (10) and Eq. (11) are dependent on p andp only (ṗ are not involved).

Determinate inverse dynamics problem using the initial dynamic equations
In case of the deterministic model of actuation, f u = B τ u τ = B τ τ , in the initial dynamic equations, Eq. (8), one can augment an n × (3 + k + l) = n × n invertible matrix Then, applying the kinematic characteristics of the observed movement, p d (t) andp d (t), the inverse dynamics solution can be stated as ⎡ from which λ rd (t), τ d (t) and λ * d (t) in all the joints can explicitly be determined. While the obtained variations of external reactions λ rd (t) and resultant muscle torques τ d (t) are exact in the sense of the modeling assumptions, the variations of joint reactions λ * d (t) are erroneous as the influence of the tensile muscle forces on the internal loads is omitted in the deterministic formulation. The determination of joint reactions should therefore follow the nondeterministic model of actuation. We will come to this problem in Sect. 3.4.

Determinate inverse dynamics problem using the projected dynamic equations
Using the deterministic model of actuation, f u = B τ u τ = B τ τ , and the projected dynamic equations obtained in Eq. (10), first, one can augment an r × (3 + k) = r × r coefficient matrixW which corresponds to W defined in Eq. (12). Then, as previously, after applying the kinematic characteristics of the observed movement, p d (t) andp d (t), the determinate inverse dynamics problem can be resolved to from which λ rd (t) and τ d (t) in the observed motion can explicitly be determined. The present formulation of determinate inverse dynamics problem corresponds to that introduced in Eq. (13), and the solutions λ rd (t) and τ d (t) are in both cases the same. The solution to Eq. (13) includes also λ * d (t), which are the joint reactions with the lack of influence of the tensile muscle forces. In the present formulation, by applying the projected dynamic equations, the joint reactions are not computed, and they must be determined separately. As it has already been mentioned in order to incorporate the contribution of muscle forces to the internal loads, their determination must follow the indeterminate inverse dynamics problem.

Indeterminate inverse dynamics problem
The indeterminate inverse dynamics problem, also widely known as the redundant problem in biomechanics or the muscle force distribution problem (see, e.g., [29][30][31][32]), relates to the task of sharing the muscular joint torques from the determinate inverse dynamics analysis into individual muscle forces, and is usually solved using optimization techniques. For the present hybrid model of actuation, u τ σ = [τ T σ T ] T , the control redundancy is limited to the supporting-leg joints, and the muscular load sharing problem relates to τ d (t) → σ d (t), where k = 3 resultant muscle torques in the supporting-leg are distributed into the effort of m = 9 muscle forces F j = σ j A j . The redundancy of muscular load sharing is usually addressed by minimizing a cost/objective function appropriately selected for the movement under investigation. Prior to this process, one must know an explicit mathematical relationship between the considered muscular forces/stresses and the torques they exert on the respective joints.
As said in Sect. 2.3, the deterministic (f u = B τ u τ ) and nondeterministic (f u = B τ σ u τ σ ) models of control are not equivalent, i.e., B τ u τ = B τ σ u τ σ , and more strictly B τ τ = B σ σ . The generalized force vector f c = −C T (p)λ in Eqs. (1)/(8) will then be different when using f u = B τ u τ or f u = B τ σ u τ σ . Therefore, the muscular load sharing problem cannot be based on these initial dynamic equations, and one must use the projected dynamic equations introduced in Eq.  [24,26]. With this projection, one receives The distributed formulation of these k = 13 relations, following the denotation of Eq. (3), is from which one can state the following explicit relationship between the k = 3 resultant muscle torques τ in the supporting-leg joints and the respective m = 9 muscle stresses σ G(p)τ = H(p)σ (18) where G = B T τ D(D T MD) −1 D T B τ and H = B T τ D(D T MD) −1 D T B σ are of dimensions k × k and k × m, respectively. Using this explicit relationship between the muscular forces/stresses in the supporting-leg and the torques they exert on the respective leg joints, the considered muscular load sharing problem can be stated as the following optimization scheme: where J is a chosen cost (objective) function, and σ min and σ max are the physiologically allowable minimal and maximal values of muscle stresses. In this way, the resultant muscle torques τ d (t) in the supporting-leg joints, known from the determinate inverse dynamics solution, are distributed into the individual muscle efforts σ d (t) in the lower limb. The cost function J used in this study is that proposed in [33] that minimizes the sum of cubed muscle stresses involved in Eq. (18). This cost function, which maximizes endurance, is commonly accepted for lower extremity analysis, with many applications to human gait [19]. For more dynamic movements like the triple jump, a more justified cost function should possibly reflect other criteria than the maximum endurance, which was not addressed in this study.

Determination of joint reactions in the supporting-leg joints
Using τ d (t) and λ rd (t) from the solution to the determinate inverse dynamics problem, and σ d (t) assessed from the indeterminate inverse dynamics problem, the joint reaction λ d (t) in the supporting lower limb joints can be determined from Eq. (11) as By contrast to the internal reactions λ * d (t) obtained from Eq. (13), the present solution λ d (t) is limited to the supporting-leg joints, and takes into account the contribution of the tensile muscle forces to the internal loads.

Anthropometric data and musculoskeletal geometry
For the developed jumper model, the anthropometric data used in the inverse dynamics simulation of triple jump are: the lengths l i of the segments, locations ξ C i and η C i of their mass centers in the local coordinate frames, and their masses m i and mass moments of inertia J C i with respect to C i , i = 1, . . . , b. The locations of the shoulder and hip joints in the local reference frames of segments 2 and 4 ( Fig. 1), and the distance from the ankle joint A to point P (Fig. 2), are also required. The lengths and the body mass can be measured directly from the subject. The segment masses and the mass center locations must then be estimated using the regression equations reported in, e.g., [32,[34][35][36], which is concerned with a series of additional measurements of characteristic circumferences and segment lengths of the subject body. In addition, the lower limb musculoskeletal model (Fig. 3) requires estimating the cross-sectional areas A j of the specified muscles, and then the actual origin and insertion point locations in the local reference frames of appropriate segments. Of special importance is the identification of the muscle paths relative the skeleton [21][22][23]. This involves, in particular, the estimation of the muscle moment arms with respect to appropriate joints, resulting in the effective origin and insertion points of the muscle forces as described in Sect. 2.3; see also [23] where a similar model was developed for the upper limb. In our studies, we estimated the data following, among others, [32,35,36]. The relevant anthropometric data used in this study for the inverse simulation of triple jump (Sect. 4.4), estimated for a 20-year-old male subject of mass 76 kg and height 179 cm, are provided in Tables 1, 2, and 3 of the Appendix.

Kinematic data
The triple jump performance was recorded using four synchronized digital cameras (100 Hz), and the three-dimensional coordinates of a set of p = 19 base points on the athlete body were converted (digitized) from the photographic images using the direct linear transformation method [37]. The base points were marked on the jumper skin so that to coincide with k = 13 model joints, and additional p − k = 6 base points were defined at the external segment tips, seen as black dots in Fig. 1. For the purpose of the present analysis, only the X and Y coordinates of the base points were used, r j = [x j y j ] T , j = 1, . . . , p, and r = [r T 1 . . . r T p ] T , so that each segment position in the sagittal plane was defined by appropriate two base points. Then, for the ith segment, using the coordinates of its two assigned base points and the segment anthropometric data, the segment absolute coordinates p i = [x C i y C i θ i ] T could easily be determined. If the length of ith segment, arising from the measurements as a distance between the two assigned base points, was different from l i introduced in the anthropometric data, the ξ C i location of mass center C i in the local reference The one-leg jump: landing and take-off from the force platform system was appropriately scaled, and the scaled value was then used in calculations of x C i and y C i (values of η C i are usually small and their scaling is meaningless). Repeating this for all segments, the transformation r(t) → p(t) can effectively be achieved. The raw kinematic data were finally smoothed to p d (t) using a second order Butterworth filter [34] with cutoff frequency of 10 Hz. The requiredp d (t) characteristics were then computed from p d (t), sampled with fixed time intervals t = 0.01 s. The acceleration at the kth sample was found numerically asp

Partial verification of the simulation model
The validity of the developed simulation model was verified through numerical tests for a one-leg jump with the right leg on the force platform (Fig. 5), similar to a transition from the step to jump in a triple jump. The movement, performed by a 22-year-old athlete of mass 67.4 kg and height 170 cm, starts with a short prejump from the left to the right leg. Then, after landing with the right leg on the force platform, the athlete takes-off with the same leg in the upward-forward direction. Vanishing values of the computed external reactions during the flight phases and the compatibility of the calculated and measured external reactions during the contact phase are taken as criteria for the validity of the simulation model. As seen from Fig. 6, the calculated ground reactions fit well the measured values, being appropriately high during the contact phase while diminishing during the flying phases. The approximate validity of the (deterministic) simulation model is thus proven, which involves also an acceptable accuracy of the measured/assessed anthropometric data of the jumper as well as the accuracy/adequacy of kinematic characteristics of the analyzed movement. On the other hand, the processed (smoothed) kinematic characteristics used, with attenuated peaks and valleys, preclude the precise identification of the impact, and the "calculated contact" with the platform begins before and ends after the actual contact known from measurements. The analysis was focused on three contact (landing and take-off) phases, together with some short flying periods before and after the contacts with the runway. The contact phases were chosen for the expected high external reaction forces, reflected in the internal loads in the lower limbs. Phase 1 starts after the approach run, just before the jumper touches the takeoff board with the takeoff leg, covers the whole subsequent contact, and, after the take-off for the hop, includes a short period of the hop flight. Likewise, Phase 2 corresponds to the transition from the hop to the step (from the first to the second jump), and Phase 3 relates to the transition from the step to the final jump. All the three phases are illustrated in  As seen in the graphs, the most dynamic stage is the third phase when the competitor takes-off for the final jump. The assessed maximal vertical reactions R y from the ground were found somewhere in the middle of the contact phases, and were equal to 3.54, 4.62, and 5.14 times bodyweight (G = 745.3 N), respectively for phase 1, 2, and 3. The maximum resultant reactions from the ground, R = R 2 x + R 2 y , were almost the same, respectively: 3.55, 4.68, and 5.23 times bodyweight, and the maximal horizontal ground reactions were found just after landing and before take-off actions, and vary between 1 and 1.5 times bodyweight. The significant fact is that the peak values of ground reactions during the triple jump reported in [4] are much bigger, c.a. 15 and 7 times bodyweight for the vertical and horizontal components in the braking part of the step, respectively. These were, however, the values measured from the force platform, and as such included the impact forces, which are not represented in our calculations due to the modeling simplifications and smoothed kinematic data (attenuated picks and valleys). There were also differences in the values of contact times. Respectively for the hop, step, and jump, the contact times reported in [4] were 0.129 s, 0.157 s, and 0.177 s, while in our study these were c.a. 0.18 s, 0.21 s, and 0.18 s, which were estimated from the film pictures. Finally, the run-up velocity reported in [4] was 8.65 ms −1 , while we measured around 7.5 ms −1 . The less velocity and longer contact times may explain the smaller reaction values, also. The increasing ground reactions, from Phase 1 (Fig. 8), through Phase 2 ( Fig. 9), and to Phase 3 (Fig. 10), are reflected in the internal loads in the supporting leg-here the right leg in Phase 1 and Phase 2, and the left leg in Phase 3. This is represented, first, in the courses of resultant muscle torques, which are qualitatively and quantitatively different in the three contact phases, with the maximum values in Phase 3. According to the definitions in Fig. 1, during the contact phases, the torques in the knee and ankle joints are positive in the extension and plantar flexion directions, respectively, while the torque in the hip joint changes its direction from acting as an extensor during the landing and as an flexor at the end of take-off. The torques are then distributed in the muscle forces following the indeterminate inverse dynamics scheme provided in Sect. 3.3. In Figs. 8-10, we demonstrate the assessed forces of only three muscles of the supporting leg: vastus, soleus, and gluteus, denoted by m2, m5, and m9 in Fig. 3, the extensors, respectively, in the knee, ankle, and hip joints. As seen in the graphs, gluteus is active mainly during the landing stages, while vastus and soleus are responsible for the take-offs. The external loads and the muscle forces contribute finally to the joint reactions. For example, in the ankle and knee joints of the supporting leg, the maximum joint reaction was estimated in the knee joint in the middle of Phase 3, whose magnitude λ = λ 2 x + λ 2 y is almost 18 times bodyweight of the jumper, and 3.6 times higher than the maximum value of ground reaction. This increase in the internal loads is due to the contribution of tensile muscle forces in the joints.
The assessed internal loads, both muscle forces and joint reactions in the supporting leg, seem high. This might be due to many reasons: underestimation of the muscle moment arms resulted from the via points, wrapping surfaces and the patellar mechanism in our rough musculoskeletal model, inaccurately estimated anthropometric data, and the simplifying assumptions such as the direct contribution of the tensile muscle forces to the joint reactions. Therefore, the reported simulation results should rather be treated as qualitative, illustrating the potential of the developed formulation.

Discussion and conclusion
A compact and effective formulation for the inverse dynamics simulation of the triple jump was presented which can also be used for the analysis of other sagittal plane movements like the long jump, walking, and running. The essential advantages of the present formulation can be enumerated as follows.
• The formulation is valid for all movement phases, irrespective of the fact whether the modeled human flies or contacts the ground with one of his legs. The double support case can possibly be analyzed after determining the way the external reactions are distributed between the two feet in contact with the ground. • The human dynamics is introduced in absolute coordinates, which results in an easy formulation of M, f g , A r , and B τ (determinate model of actuation). Only the determination of B τ σ , and particularly B σ for the lower extremities for nondeterministic model of actuation involves some modeling complexity.  10) and (11), respectively. In this way, the inverse dynamics simulation problem is split into the determination of actuation (including the ground reactions) and the determination of joint reactions. • With the hybrid model of actuation, using the muscle forces in the supporting leg and the resultant muscle torques in the other body joints, the attention can be focused on a more detailed analysis of internal loads in the lower limb while preserving the influence of the whole body motion. • The determinate inverse dynamics formulations provided in Sects. 3.1 and 3.2 allow for the explicit determination of external reactions and resultant muscle torques in all the model joints during the observed movements. • The indeterminate inverse dynamics formulation described in Sect. 3.3 is then focused on the direct distribution of the previously assessed muscle torques in the supporting leg into the respective muscle stresses/forces, resulting in the decomposed optimization scheme introduced in Eq. (19). • Following Eq. (20), the joint reactions in the supporting leg can finally be determined, involving the influence of tensile muscle forces, possible reactions from the ground, and the whole system dynamics.
The human motion apparatus is extremely complex and, consequently, very difficult to model. The simplified models used for the inverse dynamics simulation always involve major inaccuracies. The critical modeling parameters are associated with the assumptions related to the musculotendon paths and the effective attachment points of the tendons [20,21,32]. The physiological cross-sectional areas of the muscles and other anthropometric data, assessed for a particular subject with a rather large margin of uncertainty, can also significantly affect the estimated magnitudes of muscle forces. Of great importance is then the way the raw kinematic data are processed before they are used in the inverse dynamics simulation [34]. Finally, the muscle force estimates are influenced by muscle decomposition and recruitment criteria used in the force sharing optimization process. Nonetheless, though the inverse dynamics simulation is always burdened with possible high inaccuracy, it remains the only available noninvasive method for the assessment of the internal loads during human movements.