- Email: [email protected]
Modelling, simulation and optimisation of a human vertical jump
Journal of Biomechanics 32 (1999) 521—530
Modelling, simulation and optimisation of a human vertical jump T. Spa¨gele , A. Kistner , A. Gollhofer * Institute A of Mechanics, University of Stuttgart, 70550 Stuttgart, Germany Department of Sport Science, University of Stuttgart, 70569 Stuttgart, Germany Received in final form 21 September 1998
Abstract This paper describes an efficient biomechanical model of the human lower limb with the aim of simulating a real human jump movement consisting of an upword propulsion, a flying and a landing phase. A multiphase optimal control technique is used to solve the muscle force sharing problem. To understand how intermuscular control coordinates limb muscle excitations, the human body is reduced to a single lower limb consisting of three rigid bodies. The biomechanical system is activated by nine muscle—tendon actuators representing the basic properties of muscles during force generation. For the calculation of the minimal muscle excitations of the jump movement, the trajectory of the hip joint is given as a rheonomic constraint and the contact forces (ground reaction forces) are determined by force plates. Based on the designed musculoskeletal model and on the differential equations of the multibody system, muscle excitations and muscle forces necessary for a vertical jump movement are calculated. The validity of the system is assessed comparing the calculated muscle excitations with the registered surface electromyogramm (EMG) of the muscles. The achieved results indicate a close relationship between the predicted and the measured parameters. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Vertical jump movement; Human movement simulation; Optimal control; Dynamic optimisation
1. Introduction The determination of muscle forces during movement is not only essential for a complex analysis of internal loads acting on bones and joints, it also contributes to a deeper understanding of the underlying neural control. Due to the difficulty of measuring muscle forces directly within the living system by means of invasive techniques and due to the circumstance of the mechanically redundant arrangement of the actuators, static optimisation (Crowninshield, 1978; Hardt, 1987; Patriarco et al., 1981) or dynamic optimisation (Davy and Audu, 1987; Pandy et al., 1990; Spa¨gele, 1995) have been used to estimate muscle forces during movement. To produce motion of the body, the muscles must be activated by signals from the central nervous system (CNS) via electrical impulses. The electrophysiological excitations produce muscle forces, which are related to joint torques by the moment arms of those muscles. These torques cause angular accelerations, which are necessary for the calculation of the multibody movement that can be described by gene-
*Corresponding author.
ralised coordinates of the mechanical system. For complex movement simulations, the human system may be described by a model incorporating two distinct sections. The first part describes the behaviour of the actuators during force generation. The second part includes the musculoskeletal system with the mechanical properties of the limbs and the generation of joint torques as a function of the muscle forces. With the muscle excitations as results of an optimal control algorithm (Pandy et al., 1992, 1995) and an appropriate model of the musculoskeletal system, it is possible to simulate complex human movements and to resolve neural excitations that are responsible for the measured movement pattern. In order to use a multiphase optimal control approach for a satisfactory solution to the simulation problem, biomechanical models are necessary which are both accurate and highly efficient. Previous simulations (e.g. Pandy et al., 1990) investigated only the propulsion phase of the jump with the jump height as the performance criterion. The aim of this study is to simulate a measured complex jump movement consisting of an upward propulsion, a flying and a landing phase. The underlying performance criterion contains kinematic-based task criteria and neuromuscular penalty criteria.
0021-9290/99/$ — see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 1 - 9 2 9 0 ( 9 8 ) 0 0 1 4 5 - 6
522
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
2. Methods 2.1. Muscle model The muscles are the active elements controlled by the central nervous system. With their capacity of active shortening they can apply forces to the skeleton. One major problem in modelling large-scale motions is to find appropriate descriptions of the actuators with accurate mathematical models. In the literature, several very detailed and complex muscle models are described which can roughly be divided into two classes (Audu and Davy, 1985). The molecular models describe the processes of force generation on the sarcomer level; they range in complexity from single-molecular models (Huxley, 1957, 1974) to more complex approaches with many specific parameters (Hill et al., 1975; Zahalak and Ma, 1990). Despite their apparent beneficial effects on the representation of the biophysical processes, these models, however, do not appear to be practical for simulations of large-scale movements (Winters and Stark, 1987). More tractable approaches, single-element models, describe the dynamical properties of force generation of the muscles. They proved to be practical, especially in numerical applications (Spa¨gele, 1998). The muscle forces of the musculoskeletal system are represented by the vector f+. This vector is a function of the normalised muscle excitations u, the muscle lengths l+ and the muscle velocities v+ (Fig. 1). The muscle—tendon element (MTE) contains the parameters of the activation dynamics and the contraction mechanics of the muscle. The activation dynamics is considered to be related to the release of calcium ions from the sarcoplasmatic reticulum described by the vector c. Their subsequent function forming crossbridges to enable contraction, can be expressed by the active states a of the muscles (Audu and Davy, 1985; Winters, 1990). In parallel to the MTE, a damping element (DE) and a parallel elastic element (PEE)
describe the passive behaviour of the muscle. The PEE can be represented by a exponential spring, whereas the DE is modelled as a linear damper. The resultant forces f+ in the muscles comprise the forces of the three elements. They are the sum of the active muscle forces f+2#, the passive muscle forces f.## and f"#: f+ (u, l+ , v+)"f+2# (a(u), l+ , v+)#f.##(l+)#f"#(v+). (1) Based on the activation dynamics and the normalised muscle excitations u it is possible to calculate the vector a of the active states (Fig. 1). The muscle mechanics describing the dependencies of active muscle forces on the muscle lengths l+ and muscle velocities v+ contain the Hill force—velocity relationship and the tension—length relationship (Winters, 1990). The active muscle force F+2#3f+2# produced by the ith muscle group is calG culated as the product of the maximum isometric force F and the factors f (u ), f (l+ ) and f (v+ ), which
G "G G 2*G G $4G G depend on excitation and mechanical conditions of force generation (Spa¨gele, 1998) of the active muscle force. F+2#(u , l+, v+)"F f (u ) f (l+ ) f (v+ ). (2) G G G G
G "G G 2*G G $4G G The muscle model contains the basic active and passive characteristics of force generation (Fig. 1). It can be used in particular to calculate the muscle forces for the simulation and optimisation of human movements. 2.2. Multibody model Mechanical systems and objects such as vehicles, robots or even human bodies, can be modelled precisely as multibody systems. Their mechanical properties can be characterised by motion as a rigid body and they can be constrained by elements such as joints, bearings and supports. Force elements may be realised by springs, dampers or motors (Schiehlen, 1991). Mechanical principles (e.g. principle of d’Alembert in Lagrange’s version)
Fig. 1. Block diagram showing the calculation of the vector f+ of the muscle forces as a function of the normalised muscle excitations u, muscle lengths l+ and muscle velocities v+.
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
allow the equations of motion for a holonomic multibody system to be derived (Pfeiffer, 1989; Schiehlen, 1986) under these conditions. Subsequent equations can be arranged in a specific form and results for a planar multibody motion according to M(q)q¨ #gc(q, q , t)"gc(q, q , t),
(3)
where M(q) is the matrix of inertia L M(q)" (m JTTiJTi#HX JTRiJRi), (4) G 1G G gc(q, q , t) the vector of generalised gyroscopic forces L gc(q, q , t)" (m JTTia si#HX JTRiaN ), (5) G 1G G G and ge(q, q , t) the vector of generalised applied forces L ge(q, q , t)" (JTTif ei #JTRiTe,z (6) Si ). G For the calculation of Eqs. (3)—(6), the Jacobian matrices J (translation) and J (rotation) as well as the segment 0G 2G masses m , the moments of inertia HX , the local accelerG 1G ations a and angular accelerations aN , the applied forces G 1G e,z f and the torques T Si acting on the centre of masses are G necessary for all n modelled rigid bodies. The equations of motion (3) are constructive parts of the model of the musculoskeletal system (Fig. 2) to determine the generalised coordinates q, the muscle lengths l+ and the muscle velocities v+ as a function of the muscle forces f+. The anatomical arrangement of the muscles is essential for
523
the evaluation of the driving forces and the joint torques acting on the single-body segments. For the muscle force generation, there exists a subtle interaction between the required normalised excitations, the muscle lengths and the muscle velocities. With knowledge of the muscle forces and their moment arms, it is possible to calculate the active joint torques t contributing with the passive joint torques t to the joint torques t responsible for the active movement of the multibody system (Fig. 2). The functions of the moment arms d, the passive joint torques t, the muscle lengths l+ and muscle velocities v+ depend on the generalised coordinates q and their first time derivatives q . They are given in detail in Spa¨gele (1998). 2.3. Complete state space model The kinetics of the skeleton as well as the excitation—contraction dynamics of the muscles are determined by a set 2f#n differential equations of first order + x "f (x, u, t), x(t)"x0, x3IRD>L+, (7) where f is the degree of freedom of the mechanical system and n is the number of muscles modelled in the + musculoskeletal system. For state space notation, the vector q of the generalised coordinates is described by x , the vector q of the generalised velocities by x , and the normalised calcium ion concentrations c of the muscles by the vector x . The state vector x describing the comA plete biomechanical behaviour of the human body is
Fig. 2. Block diagram showing the calculation of the generalised coordinates q, the muscle lengths l+, and the muscle velocities v+ as a function of the muscle forces f+.
524
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
represented by first-order differential equations (ODE): x xv !1 e x " x " M (xp)(g (xp, xu, x , t)!gc(xp, xv, t)) A x C(u, x ) A A
(8)
where the vector C contains the excitation functions of the modelled muscles and describes the characteristic properties of all muscles during activation. With the state vector x of the musculoskeletal system and the dynamic optimisation approach, it is possible to compute the optimal controls u* responsible for a measured human movement with minimal cost functional. For a prescribed type of movement, the optimisation problem requires the control functions, the model design parameters and the phase separation times of the distinct movement phases to minimise a scalar cost functional I, which may be described by a combination of Mayer terms and Lagrange terms L (Ja¨nsch et al., 1989).
H K LH dt . (9) I" # H# RH\ H Additionally to differential equations (8) the solution is required to satisfy separable multipoint boundary conditions at the initial time of the movement as well as at the phase time points, the pure parameter constraints and the path constraints of the multiphase optimal control problem. The subdivision of the entire movement into m-independent phases requires certain conditions at their respective phase transients, which describe the connections of the single phases.
3. Application As an application, a vertical one-legged jump exercise is investigated to simulate human movement. This jump movement can be divided into an upward propulsion, a flying and a landing phase. In order to restrict the complexity of the musculoskeletal system, the active lower limb is described separately (Fig. 3). This approach requires both knowledge of the three-dimensional kinetics of the body (Fig. 4) and the ground reaction forces (see Appendix, Fig. 8). It allows, however, a detailed description of the reduced biomechanical system using more comprehensive musculoskeletal models. For simulation of the vertical jump, the reduced model consists of a single lower limb represented by three rigid bodies. The segments of the thigh, the shank and the foot are joined together by frictionless joints. They are allowed to be moved in the sagital plane by pure rotation around the hip, knee and ankle joint. The biomechanical system of the leg is activated by nine muscle—tendon actuators (Fig. 3). The force generation of each muscle group is described by the three-element muscle model.
With respect to the rheonomic constraints (measured hip point acceleration) of the hip joint, the mechanical system of the lower limb possesses f"3 degrees of freedom for the plane motion. The kinetics of the multibody system can be described by the independent absolute angles q , q and q of the thigh, shank and foot, repre sented by the vector q of the generalised coordinates q"[q q q ]2. (10) For the description of the entire biomechanical system of the jump movement, 15 first-order differentials are necessary incorporating six equations of motion and nine activation equations of the modelled muscle groups. In the approach described in this manuscript, the optimal control problem is defined to find the muscle excitations u minimising the scalar cost functional I . This cost functional contains not only the weighted differences between the measured and calculated kinematic data of the leg motion, but also the weighted sum of all neuromuscular values of muscle excitations I"I #I RH " wH ((u!u )#(t!t ) H\ R H RH u dt. #(m!m )) dt# wH (11) H\ G R H G The equations of motion for the jump movement are set up applying the principle of d’Alembert in Lagrange’s version. They can be specified by three non-linear differentials of second order (see Appendix A). The biological part of the system consists of the muscle model and the arrangement of the nine actuators within the musculoskeletal system. Together with the muscle excitations u these are inputs to the moments t acting at the joints, which contain the active t and the passive t joint moments responsible for the body motion.
¹ ¹ #¹ t"t #t" ¹ " ¹ #¹ . (12) ¹ ¹ #¹ The active joint moments t are the sum of the products resulting from the muscle forces f+"[F+, 2 , F+]2 and their respective moment arms d. For the jump motion, the active joint moments at the hip, knee and ankle joint can be as follows: u
F+du#F+d !F+du!F+du ta" F+dR#F+dR!F+dR!F+dR!F+dR , (13) F+dK !F+dK !F+dK 0 where d describes the moment arm of the ith muscle G group acting at the joint with the relative angle 03+u, t, m,. The moment arms d of the muscles, the passive moments t of the system and the muscle kinetics l+ and
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
525
Fig. 3. Schematic representation of the measured human vertical jump with one leg and the musculoskeletal model of the lower limb activated by nine muscle groups.
v+ are functions of the relative angles. The detailed functions and the anatomic parameters of the subjects investigated in this study are given in Spa¨gele (1998). The angles u, t and m (Fig. 4) describe the relative position of the segments and can be calculated from the absolute coordinates q of the mechanical system u"q ,
t"q !q ,
n m"q !q # . 2
(14)
4. Results and discussion The aim of the study was to simulate a complex jump movement consisting of an upward propulsion, a flying and a landing phase. The used performance criterion
contained kinematic based as well as neuromuscular penalty criteria (Zajac and Winters, 1990). Previous simulations (e.g. Pandy et al., 1990) of jumping have concentrated only on the propulsion phase with the jump height as the performance criterion. The basic result of the dynamic optimisation is shown in Fig. 5. The differences between calculated and measured trajectories are considerably small, according to the definition of the cost functional. This paper focuses on a description of how the biomechanical system controls the single muscle groups in order to perform the movement with minimal muscle excitations. The normalised muscle excitations of the modelled musculoskeletal system correspond to the elements of the control vector u of the optimal control algorithm. The results of the simulated vertical jump
526
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
Fig. 5. Comparison of measured (dotted line) and calculated (solid line) trajectory of the hip angle, the knee and the ankle angle of the vertical jump.
Fig. 4. Schematic representation of the relative angles u, t and m acting on the lower limb.
movement are summarised in Fig. 6. The calculated time trajectories of the minimal muscle excitation functions are presented for each muscle group separately. These computed optimal controls can be verified and validated by comparison with the actual muscle action potentials registered with surface electromyography (Fig. 7). The entire jump has been divided into three distinct phases and this is also reflected in the determined muscle excitations. In the first part of the jump, a plantar movement of the foot must be performed for the acceleration of the body in vertical direction. This motion is produced primarily by the gastrocnemius (muscle group 7) and the soleus (muscle group 8). Therefore, strong excitations in these two groups with maximal excitation rates at the end of the first phase are characteristic (Fig. 6). As tibialis anterior (muscle group 9) counteracts the plantar motion of the foot this muscle does not seem to be activated
during the first part of the jump. During vertical acceleration, extension of the shank is produced by simultaneous action of the rectus femoris (muscle group 2) and the vastus group (muscle group 3). This is indicated by rather similar trajectories characterised by increasing excitation rates. During the first phase, extension of the thigh can be produced by the gluteus group (muscle group 4) and the hamstring group (muscle group 5). The hamstrings and the biceps femoris c. b. (muscle group 6) counteract the required motion of the shank. They are therefore not controlled during the entire jump movement. The dynamics of the thigh is primarily generated by the gluteus group, as shown by an increasing excitation during the first jump phase. During flight phase only few muscles show considerable activation amplitudes. The largest torque acts in the knee joint at the beginning of the landing phase. Therefore, intense activation of the second (rectus femoris) and third (vastus) muscle group must be expected for deceleration of the body after landing. This is also achieved by the ability of the muscles to develop high forces during eccentric muscle actions. Due to the properties of the contractile apparatus in enhancing the force output during active lengthening the calculated
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
527
Fig. 6. Trajectories of the calculated control values u of the nine modelled muscle groups necessary for the vertical jump movement with minimal muscle excitations.
time histories of the muscle excitations are reduced in the deceleration phase as compared with the amplitudes calculated for the vertical acceleration of the body. Even preactivation can be observed in the plantar flexors in the short time interval before landing. Functionally preactivation is necessary to stiffen the joint complexes before mechanical loading (Gollhofer et al., 1992). The human musculoskeletal system is characterised by non-linear skeletal and excitation—contraction dynamics, high dimensionality of the state vector and a great number of control functions. Therefore, dynamic optimisation problems of human movement simulation are very
time expensive. The calculation time of the presented jump movement depends on the initial conditions as well as on the determination of the biomechanical parameters which must be frequently adapted to the problem during the simulation process. The optimisation of a single upward propulsion has already been published in the literature (Pandy et al., 1990). In the present study, however, a multi-phase motion consisting of an upward propulsion, a flying and a landing phase is simulated. To solve the described optimisation problem, it is necessary to register the contact forces and the kinetics of the jump movement. The
528
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
Fig. 7. Experimental EMG activities of the leg muscle groups measured during the jump movement.
presented biomechanical model has been validated with empirical data from the excitation of the entire muscular chain of the lower limb during the multi-phase jump performance, including the flying and landing phase. For the simulation of human movements, it is necessary to use efficient and attractive models of the human muscular system. The developed muscle model with its three-elements is especially attractive for the computational investigation of large-scale movements. The validity of the simulations can be assessed by comparing the calculated muscle excitations with originally registered surface EMG of the muscles of interest. The achieved results indicate a close relationship between the predicted and the measured parameters.
Appendix A A.1. Multibody dynamics The equations of motion describing the mechanical behaviour of the skeleton are given below. The differential equations of the multibody system can be generally written as M(q)q¨ #g(q, q , t)"g(q, q ,t). The inertia matrix M of the system M M M M" M M M , M M M
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
529
M "H #m z#(m #m )l, 1 M "(m z #m l )l cos(q !q ), M "m l (z cos(q !q )#z sin(q !q )), M "H #m z#m l, 1 M "m l (z cos(q !q )#z sin(q !q )), M "H #m (z#z). 1 The vector g of the generalised gyroscopic forces of the system g" (m z #m l )l sin(q !q )qR #m l (z sin(q !q )!z cos(q !q ))qR #(m z #(m #m )l )(cos(q )a (t)#sin(q )a (t)) W X (m z #m l )l sin(q !q )qR #m l (z sin(q !q )!z cos(q !q ))qR #(m z #m l )(cos(q )a (t)#sin(q )a (t)) W X
.
m l (z sin(q !q )#z cos(q !q ))qR #m l (z sin(q !q )#z cos(q !q ))qR #m ((z cos(q )!z sin(q ))a (t)#(z sin(q )#z cos(q ))a (t)) W X
The vector g of the generalised external forces of the system g" F (t)l cos(q )#((F (t)!(m #m )g)l !m gz ) sin(q )#¹ !¹ W X F (t)l cos(q )#((F (t)!m g)l !m gz ) sin(q )#¹ !¹ W X
Fig. 8. Measured trajectories of the ground reaction forces F (t) and W F (t) of the vertical jump movement and the distance ¸(t) from the ankle X joint to the ground reaction forces. .
F (t)(l cos(q )!(L(t) sin(q ))#F (t)(l sin(q )#L(t) cos(q )) W X !m g(z sin(q )#z cos(q ))#¹
The mechanical definitions are: a (t) measured horizontal acceleration of the hip W point a (t) measured vertical acceleration of the hip point X F (t) measured horizontal ground reaction force W F (t) measured vertical ground reaction force X g gravity constant H moment of inertia of segment i about its centre 1G of mass l length of segment i G ¸(t) measured distance from the ankle joint to the ground reaction forces m mass of the segment i G q generalised coordinates of the mechanical sysG tem t time ¹ torque acting at the joint i G z distance from distal point to centre of mass of G segment i
References Audu, M.L., Davy, D.T., 1985. The influence of muscle model complexity in musculoskeletal motion modelling. Journal of Biomedical Engineering 107, 147—157. Crowninshield R.D., 1978. Use of optimization techniques to predict muscle forces, Journal of Biomedical Engineering 100, 88—92. Davy, D.T., Audu, M.L., 1987. A dynamic optimization technique for predicting muscle forces in the swing phase of gait. Journal of Biomechanics 20, 187—201. Gollhofer, A., Strojnik, V., Rapp, W., Schweizer, L., 1992. Behavior of triceps surae muscle—tendon complex in different jump conditions. European Journal of Applied Physiology 64, 283—291. Hardt D.E., 1978. Determining muscle forces in the leg during normal human walking — an application and evaluation of optimization methods. Journal of Biomedical Engineering 100, 72—78. Hill, T.L., Eisenberg, E., Chen, Y., Podolsky, R.J., 1975. Some selfconsistent two-state sliding filament models of muscle contraction. Journal of Biophysiology 15, 335—372. Huxley, A.F., 1957. Muscle structure and theories of contraction. Progress in Biophysics and Biophysical Chemistry 7, 257—318. Huxley, A.F., 1974. Muscular contraction. Journal of Physiology 243, 1—43. Pandy, M.G., Anderson, F.C., Hull, D.G., 1992. A parameter optimization approach for the optimal control of large-scale musculoskeletal systems. Journal of Biomedical Engineering 114, 450—460.
530
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
Pandy, M.G., Garner, B.A., Anderson, F.C., 1995. Optimal control of non-ballistic muscular movements: a constraint-based performance criterion for rising from a chair. Journal of Biomedical Engineering 117, 15—26. Pandy, M.G., Zajac, F.E., Sim, E., Levine, W.S., 1990. An optimal control model for maximum-height human jumping. Journal of Biomechanics 23, 1185—1198. Patriarco, A.G., Mann R.W., Simon, S.R., Mansour, J.M., 1981. An evaluation of the approaches of optimization models in the prediction of muscle forces during human gait. Journal of Biomechanics 14, 513—525. Pfeiffer, F., 1989. Einfu¨hrung in die Dynamik. Teubner, B.G., Stuttgart. Schiehlen, W., 1986. Technische Dynamik. Teubner, B.G., Stuttgart. Schiehlen, W., 1991. Computational aspects in multibody system dynamics. Computer methods. Applied Mechanics and Engineering 90, 569—582. Spa¨gele, T., 1998. Modellierung, Simulation und Optimierung menschlicher Bewegungen. Ph.D.-Thesis, Institute A of Mechanics, University of Stuttgart.
Spa¨gele, T., Kistner, A., Gollhofer, A., 1995. An optimal control technique for the optimization of human movements. In: Ha¨kkinen, K., Keskinen, K.L., Komi P.V., Mero, A. (Eds.), Book of Abstracts XVth ISB Congress. Gummerus Printing, Jyva¨skyla¨, Finland, pp. 870—871. Winters, J.M., 1990. Hill-based muscle models: a systems engineering perspective. In: Winters, J.M., Woo, S.L.Y. (Eds.), Multiple Muscle Systems: Biomechanics and Movement Organization. Springer, New York, pp. 69—93. Winters, J.M., Stark, L., 1987. Muscle models: what is gained and what is lost by varying model complexity. Biological Cybernetics 55, 403—420. Zahalak, G.I., Ma, S.-P., 1990. Muscle activation and contraction: constitutive relations based directly on cross-bridge kinetics. Journal of Biomechanical Engineering 112, 52—62. Zajac, F.E., Winters, J.M., 1990. Modelling musculoskeletal movement systems. In: Winters, J.M., Woo, S.L.Y. (Eds.), Multiple Muscle Systems: Biomechanics and Movement Organization. Springer, New York, pp. 121—148.
Modelling, simulation and optimisation of a human vertical jump T. Spa¨gele , A. Kistner , A. Gollhofer * Institute A of Mechanics, University of Stuttgart, 70550 Stuttgart, Germany Department of Sport Science, University of Stuttgart, 70569 Stuttgart, Germany Received in final form 21 September 1998
Abstract This paper describes an efficient biomechanical model of the human lower limb with the aim of simulating a real human jump movement consisting of an upword propulsion, a flying and a landing phase. A multiphase optimal control technique is used to solve the muscle force sharing problem. To understand how intermuscular control coordinates limb muscle excitations, the human body is reduced to a single lower limb consisting of three rigid bodies. The biomechanical system is activated by nine muscle—tendon actuators representing the basic properties of muscles during force generation. For the calculation of the minimal muscle excitations of the jump movement, the trajectory of the hip joint is given as a rheonomic constraint and the contact forces (ground reaction forces) are determined by force plates. Based on the designed musculoskeletal model and on the differential equations of the multibody system, muscle excitations and muscle forces necessary for a vertical jump movement are calculated. The validity of the system is assessed comparing the calculated muscle excitations with the registered surface electromyogramm (EMG) of the muscles. The achieved results indicate a close relationship between the predicted and the measured parameters. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Vertical jump movement; Human movement simulation; Optimal control; Dynamic optimisation
1. Introduction The determination of muscle forces during movement is not only essential for a complex analysis of internal loads acting on bones and joints, it also contributes to a deeper understanding of the underlying neural control. Due to the difficulty of measuring muscle forces directly within the living system by means of invasive techniques and due to the circumstance of the mechanically redundant arrangement of the actuators, static optimisation (Crowninshield, 1978; Hardt, 1987; Patriarco et al., 1981) or dynamic optimisation (Davy and Audu, 1987; Pandy et al., 1990; Spa¨gele, 1995) have been used to estimate muscle forces during movement. To produce motion of the body, the muscles must be activated by signals from the central nervous system (CNS) via electrical impulses. The electrophysiological excitations produce muscle forces, which are related to joint torques by the moment arms of those muscles. These torques cause angular accelerations, which are necessary for the calculation of the multibody movement that can be described by gene-
*Corresponding author.
ralised coordinates of the mechanical system. For complex movement simulations, the human system may be described by a model incorporating two distinct sections. The first part describes the behaviour of the actuators during force generation. The second part includes the musculoskeletal system with the mechanical properties of the limbs and the generation of joint torques as a function of the muscle forces. With the muscle excitations as results of an optimal control algorithm (Pandy et al., 1992, 1995) and an appropriate model of the musculoskeletal system, it is possible to simulate complex human movements and to resolve neural excitations that are responsible for the measured movement pattern. In order to use a multiphase optimal control approach for a satisfactory solution to the simulation problem, biomechanical models are necessary which are both accurate and highly efficient. Previous simulations (e.g. Pandy et al., 1990) investigated only the propulsion phase of the jump with the jump height as the performance criterion. The aim of this study is to simulate a measured complex jump movement consisting of an upward propulsion, a flying and a landing phase. The underlying performance criterion contains kinematic-based task criteria and neuromuscular penalty criteria.
0021-9290/99/$ — see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 1 - 9 2 9 0 ( 9 8 ) 0 0 1 4 5 - 6
522
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
2. Methods 2.1. Muscle model The muscles are the active elements controlled by the central nervous system. With their capacity of active shortening they can apply forces to the skeleton. One major problem in modelling large-scale motions is to find appropriate descriptions of the actuators with accurate mathematical models. In the literature, several very detailed and complex muscle models are described which can roughly be divided into two classes (Audu and Davy, 1985). The molecular models describe the processes of force generation on the sarcomer level; they range in complexity from single-molecular models (Huxley, 1957, 1974) to more complex approaches with many specific parameters (Hill et al., 1975; Zahalak and Ma, 1990). Despite their apparent beneficial effects on the representation of the biophysical processes, these models, however, do not appear to be practical for simulations of large-scale movements (Winters and Stark, 1987). More tractable approaches, single-element models, describe the dynamical properties of force generation of the muscles. They proved to be practical, especially in numerical applications (Spa¨gele, 1998). The muscle forces of the musculoskeletal system are represented by the vector f+. This vector is a function of the normalised muscle excitations u, the muscle lengths l+ and the muscle velocities v+ (Fig. 1). The muscle—tendon element (MTE) contains the parameters of the activation dynamics and the contraction mechanics of the muscle. The activation dynamics is considered to be related to the release of calcium ions from the sarcoplasmatic reticulum described by the vector c. Their subsequent function forming crossbridges to enable contraction, can be expressed by the active states a of the muscles (Audu and Davy, 1985; Winters, 1990). In parallel to the MTE, a damping element (DE) and a parallel elastic element (PEE)
describe the passive behaviour of the muscle. The PEE can be represented by a exponential spring, whereas the DE is modelled as a linear damper. The resultant forces f+ in the muscles comprise the forces of the three elements. They are the sum of the active muscle forces f+2#, the passive muscle forces f.## and f"#: f+ (u, l+ , v+)"f+2# (a(u), l+ , v+)#f.##(l+)#f"#(v+). (1) Based on the activation dynamics and the normalised muscle excitations u it is possible to calculate the vector a of the active states (Fig. 1). The muscle mechanics describing the dependencies of active muscle forces on the muscle lengths l+ and muscle velocities v+ contain the Hill force—velocity relationship and the tension—length relationship (Winters, 1990). The active muscle force F+2#3f+2# produced by the ith muscle group is calG culated as the product of the maximum isometric force F and the factors f (u ), f (l+ ) and f (v+ ), which
G "G G 2*G G $4G G depend on excitation and mechanical conditions of force generation (Spa¨gele, 1998) of the active muscle force. F+2#(u , l+, v+)"F f (u ) f (l+ ) f (v+ ). (2) G G G G
G "G G 2*G G $4G G The muscle model contains the basic active and passive characteristics of force generation (Fig. 1). It can be used in particular to calculate the muscle forces for the simulation and optimisation of human movements. 2.2. Multibody model Mechanical systems and objects such as vehicles, robots or even human bodies, can be modelled precisely as multibody systems. Their mechanical properties can be characterised by motion as a rigid body and they can be constrained by elements such as joints, bearings and supports. Force elements may be realised by springs, dampers or motors (Schiehlen, 1991). Mechanical principles (e.g. principle of d’Alembert in Lagrange’s version)
Fig. 1. Block diagram showing the calculation of the vector f+ of the muscle forces as a function of the normalised muscle excitations u, muscle lengths l+ and muscle velocities v+.
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
allow the equations of motion for a holonomic multibody system to be derived (Pfeiffer, 1989; Schiehlen, 1986) under these conditions. Subsequent equations can be arranged in a specific form and results for a planar multibody motion according to M(q)q¨ #gc(q, q , t)"gc(q, q , t),
(3)
where M(q) is the matrix of inertia L M(q)" (m JTTiJTi#HX JTRiJRi), (4) G 1G G gc(q, q , t) the vector of generalised gyroscopic forces L gc(q, q , t)" (m JTTia si#HX JTRiaN ), (5) G 1G G G and ge(q, q , t) the vector of generalised applied forces L ge(q, q , t)" (JTTif ei #JTRiTe,z (6) Si ). G For the calculation of Eqs. (3)—(6), the Jacobian matrices J (translation) and J (rotation) as well as the segment 0G 2G masses m , the moments of inertia HX , the local accelerG 1G ations a and angular accelerations aN , the applied forces G 1G e,z f and the torques T Si acting on the centre of masses are G necessary for all n modelled rigid bodies. The equations of motion (3) are constructive parts of the model of the musculoskeletal system (Fig. 2) to determine the generalised coordinates q, the muscle lengths l+ and the muscle velocities v+ as a function of the muscle forces f+. The anatomical arrangement of the muscles is essential for
523
the evaluation of the driving forces and the joint torques acting on the single-body segments. For the muscle force generation, there exists a subtle interaction between the required normalised excitations, the muscle lengths and the muscle velocities. With knowledge of the muscle forces and their moment arms, it is possible to calculate the active joint torques t contributing with the passive joint torques t to the joint torques t responsible for the active movement of the multibody system (Fig. 2). The functions of the moment arms d, the passive joint torques t, the muscle lengths l+ and muscle velocities v+ depend on the generalised coordinates q and their first time derivatives q . They are given in detail in Spa¨gele (1998). 2.3. Complete state space model The kinetics of the skeleton as well as the excitation—contraction dynamics of the muscles are determined by a set 2f#n differential equations of first order + x "f (x, u, t), x(t)"x0, x3IRD>L+, (7) where f is the degree of freedom of the mechanical system and n is the number of muscles modelled in the + musculoskeletal system. For state space notation, the vector q of the generalised coordinates is described by x , the vector q of the generalised velocities by x , and the normalised calcium ion concentrations c of the muscles by the vector x . The state vector x describing the comA plete biomechanical behaviour of the human body is
Fig. 2. Block diagram showing the calculation of the generalised coordinates q, the muscle lengths l+, and the muscle velocities v+ as a function of the muscle forces f+.
524
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
represented by first-order differential equations (ODE): x xv !1 e x " x " M (xp)(g (xp, xu, x , t)!gc(xp, xv, t)) A x C(u, x ) A A
(8)
where the vector C contains the excitation functions of the modelled muscles and describes the characteristic properties of all muscles during activation. With the state vector x of the musculoskeletal system and the dynamic optimisation approach, it is possible to compute the optimal controls u* responsible for a measured human movement with minimal cost functional. For a prescribed type of movement, the optimisation problem requires the control functions, the model design parameters and the phase separation times of the distinct movement phases to minimise a scalar cost functional I, which may be described by a combination of Mayer terms and Lagrange terms L (Ja¨nsch et al., 1989).
H K LH dt . (9) I" # H# RH\ H Additionally to differential equations (8) the solution is required to satisfy separable multipoint boundary conditions at the initial time of the movement as well as at the phase time points, the pure parameter constraints and the path constraints of the multiphase optimal control problem. The subdivision of the entire movement into m-independent phases requires certain conditions at their respective phase transients, which describe the connections of the single phases.
3. Application As an application, a vertical one-legged jump exercise is investigated to simulate human movement. This jump movement can be divided into an upward propulsion, a flying and a landing phase. In order to restrict the complexity of the musculoskeletal system, the active lower limb is described separately (Fig. 3). This approach requires both knowledge of the three-dimensional kinetics of the body (Fig. 4) and the ground reaction forces (see Appendix, Fig. 8). It allows, however, a detailed description of the reduced biomechanical system using more comprehensive musculoskeletal models. For simulation of the vertical jump, the reduced model consists of a single lower limb represented by three rigid bodies. The segments of the thigh, the shank and the foot are joined together by frictionless joints. They are allowed to be moved in the sagital plane by pure rotation around the hip, knee and ankle joint. The biomechanical system of the leg is activated by nine muscle—tendon actuators (Fig. 3). The force generation of each muscle group is described by the three-element muscle model.
With respect to the rheonomic constraints (measured hip point acceleration) of the hip joint, the mechanical system of the lower limb possesses f"3 degrees of freedom for the plane motion. The kinetics of the multibody system can be described by the independent absolute angles q , q and q of the thigh, shank and foot, repre sented by the vector q of the generalised coordinates q"[q q q ]2. (10) For the description of the entire biomechanical system of the jump movement, 15 first-order differentials are necessary incorporating six equations of motion and nine activation equations of the modelled muscle groups. In the approach described in this manuscript, the optimal control problem is defined to find the muscle excitations u minimising the scalar cost functional I . This cost functional contains not only the weighted differences between the measured and calculated kinematic data of the leg motion, but also the weighted sum of all neuromuscular values of muscle excitations I"I #I RH " wH ((u!u )#(t!t ) H\ R H RH u dt. #(m!m )) dt# wH (11) H\ G R H G The equations of motion for the jump movement are set up applying the principle of d’Alembert in Lagrange’s version. They can be specified by three non-linear differentials of second order (see Appendix A). The biological part of the system consists of the muscle model and the arrangement of the nine actuators within the musculoskeletal system. Together with the muscle excitations u these are inputs to the moments t acting at the joints, which contain the active t and the passive t joint moments responsible for the body motion.
¹ ¹ #¹ t"t #t" ¹ " ¹ #¹ . (12) ¹ ¹ #¹ The active joint moments t are the sum of the products resulting from the muscle forces f+"[F+, 2 , F+]2 and their respective moment arms d. For the jump motion, the active joint moments at the hip, knee and ankle joint can be as follows: u
F+du#F+d !F+du!F+du ta" F+dR#F+dR!F+dR!F+dR!F+dR , (13) F+dK !F+dK !F+dK 0 where d describes the moment arm of the ith muscle G group acting at the joint with the relative angle 03+u, t, m,. The moment arms d of the muscles, the passive moments t of the system and the muscle kinetics l+ and
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
525
Fig. 3. Schematic representation of the measured human vertical jump with one leg and the musculoskeletal model of the lower limb activated by nine muscle groups.
v+ are functions of the relative angles. The detailed functions and the anatomic parameters of the subjects investigated in this study are given in Spa¨gele (1998). The angles u, t and m (Fig. 4) describe the relative position of the segments and can be calculated from the absolute coordinates q of the mechanical system u"q ,
t"q !q ,
n m"q !q # . 2
(14)
4. Results and discussion The aim of the study was to simulate a complex jump movement consisting of an upward propulsion, a flying and a landing phase. The used performance criterion
contained kinematic based as well as neuromuscular penalty criteria (Zajac and Winters, 1990). Previous simulations (e.g. Pandy et al., 1990) of jumping have concentrated only on the propulsion phase with the jump height as the performance criterion. The basic result of the dynamic optimisation is shown in Fig. 5. The differences between calculated and measured trajectories are considerably small, according to the definition of the cost functional. This paper focuses on a description of how the biomechanical system controls the single muscle groups in order to perform the movement with minimal muscle excitations. The normalised muscle excitations of the modelled musculoskeletal system correspond to the elements of the control vector u of the optimal control algorithm. The results of the simulated vertical jump
526
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
Fig. 5. Comparison of measured (dotted line) and calculated (solid line) trajectory of the hip angle, the knee and the ankle angle of the vertical jump.
Fig. 4. Schematic representation of the relative angles u, t and m acting on the lower limb.
movement are summarised in Fig. 6. The calculated time trajectories of the minimal muscle excitation functions are presented for each muscle group separately. These computed optimal controls can be verified and validated by comparison with the actual muscle action potentials registered with surface electromyography (Fig. 7). The entire jump has been divided into three distinct phases and this is also reflected in the determined muscle excitations. In the first part of the jump, a plantar movement of the foot must be performed for the acceleration of the body in vertical direction. This motion is produced primarily by the gastrocnemius (muscle group 7) and the soleus (muscle group 8). Therefore, strong excitations in these two groups with maximal excitation rates at the end of the first phase are characteristic (Fig. 6). As tibialis anterior (muscle group 9) counteracts the plantar motion of the foot this muscle does not seem to be activated
during the first part of the jump. During vertical acceleration, extension of the shank is produced by simultaneous action of the rectus femoris (muscle group 2) and the vastus group (muscle group 3). This is indicated by rather similar trajectories characterised by increasing excitation rates. During the first phase, extension of the thigh can be produced by the gluteus group (muscle group 4) and the hamstring group (muscle group 5). The hamstrings and the biceps femoris c. b. (muscle group 6) counteract the required motion of the shank. They are therefore not controlled during the entire jump movement. The dynamics of the thigh is primarily generated by the gluteus group, as shown by an increasing excitation during the first jump phase. During flight phase only few muscles show considerable activation amplitudes. The largest torque acts in the knee joint at the beginning of the landing phase. Therefore, intense activation of the second (rectus femoris) and third (vastus) muscle group must be expected for deceleration of the body after landing. This is also achieved by the ability of the muscles to develop high forces during eccentric muscle actions. Due to the properties of the contractile apparatus in enhancing the force output during active lengthening the calculated
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
527
Fig. 6. Trajectories of the calculated control values u of the nine modelled muscle groups necessary for the vertical jump movement with minimal muscle excitations.
time histories of the muscle excitations are reduced in the deceleration phase as compared with the amplitudes calculated for the vertical acceleration of the body. Even preactivation can be observed in the plantar flexors in the short time interval before landing. Functionally preactivation is necessary to stiffen the joint complexes before mechanical loading (Gollhofer et al., 1992). The human musculoskeletal system is characterised by non-linear skeletal and excitation—contraction dynamics, high dimensionality of the state vector and a great number of control functions. Therefore, dynamic optimisation problems of human movement simulation are very
time expensive. The calculation time of the presented jump movement depends on the initial conditions as well as on the determination of the biomechanical parameters which must be frequently adapted to the problem during the simulation process. The optimisation of a single upward propulsion has already been published in the literature (Pandy et al., 1990). In the present study, however, a multi-phase motion consisting of an upward propulsion, a flying and a landing phase is simulated. To solve the described optimisation problem, it is necessary to register the contact forces and the kinetics of the jump movement. The
528
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
Fig. 7. Experimental EMG activities of the leg muscle groups measured during the jump movement.
presented biomechanical model has been validated with empirical data from the excitation of the entire muscular chain of the lower limb during the multi-phase jump performance, including the flying and landing phase. For the simulation of human movements, it is necessary to use efficient and attractive models of the human muscular system. The developed muscle model with its three-elements is especially attractive for the computational investigation of large-scale movements. The validity of the simulations can be assessed by comparing the calculated muscle excitations with originally registered surface EMG of the muscles of interest. The achieved results indicate a close relationship between the predicted and the measured parameters.
Appendix A A.1. Multibody dynamics The equations of motion describing the mechanical behaviour of the skeleton are given below. The differential equations of the multibody system can be generally written as M(q)q¨ #g(q, q , t)"g(q, q ,t). The inertia matrix M of the system M M M M" M M M , M M M
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
529
M "H #m z#(m #m )l, 1 M "(m z #m l )l cos(q !q ), M "m l (z cos(q !q )#z sin(q !q )), M "H #m z#m l, 1 M "m l (z cos(q !q )#z sin(q !q )), M "H #m (z#z). 1 The vector g of the generalised gyroscopic forces of the system g" (m z #m l )l sin(q !q )qR #m l (z sin(q !q )!z cos(q !q ))qR #(m z #(m #m )l )(cos(q )a (t)#sin(q )a (t)) W X (m z #m l )l sin(q !q )qR #m l (z sin(q !q )!z cos(q !q ))qR #(m z #m l )(cos(q )a (t)#sin(q )a (t)) W X
.
m l (z sin(q !q )#z cos(q !q ))qR #m l (z sin(q !q )#z cos(q !q ))qR #m ((z cos(q )!z sin(q ))a (t)#(z sin(q )#z cos(q ))a (t)) W X
The vector g of the generalised external forces of the system g" F (t)l cos(q )#((F (t)!(m #m )g)l !m gz ) sin(q )#¹ !¹ W X F (t)l cos(q )#((F (t)!m g)l !m gz ) sin(q )#¹ !¹ W X
Fig. 8. Measured trajectories of the ground reaction forces F (t) and W F (t) of the vertical jump movement and the distance ¸(t) from the ankle X joint to the ground reaction forces. .
F (t)(l cos(q )!(L(t) sin(q ))#F (t)(l sin(q )#L(t) cos(q )) W X !m g(z sin(q )#z cos(q ))#¹
The mechanical definitions are: a (t) measured horizontal acceleration of the hip W point a (t) measured vertical acceleration of the hip point X F (t) measured horizontal ground reaction force W F (t) measured vertical ground reaction force X g gravity constant H moment of inertia of segment i about its centre 1G of mass l length of segment i G ¸(t) measured distance from the ankle joint to the ground reaction forces m mass of the segment i G q generalised coordinates of the mechanical sysG tem t time ¹ torque acting at the joint i G z distance from distal point to centre of mass of G segment i
References Audu, M.L., Davy, D.T., 1985. The influence of muscle model complexity in musculoskeletal motion modelling. Journal of Biomedical Engineering 107, 147—157. Crowninshield R.D., 1978. Use of optimization techniques to predict muscle forces, Journal of Biomedical Engineering 100, 88—92. Davy, D.T., Audu, M.L., 1987. A dynamic optimization technique for predicting muscle forces in the swing phase of gait. Journal of Biomechanics 20, 187—201. Gollhofer, A., Strojnik, V., Rapp, W., Schweizer, L., 1992. Behavior of triceps surae muscle—tendon complex in different jump conditions. European Journal of Applied Physiology 64, 283—291. Hardt D.E., 1978. Determining muscle forces in the leg during normal human walking — an application and evaluation of optimization methods. Journal of Biomedical Engineering 100, 72—78. Hill, T.L., Eisenberg, E., Chen, Y., Podolsky, R.J., 1975. Some selfconsistent two-state sliding filament models of muscle contraction. Journal of Biophysiology 15, 335—372. Huxley, A.F., 1957. Muscle structure and theories of contraction. Progress in Biophysics and Biophysical Chemistry 7, 257—318. Huxley, A.F., 1974. Muscular contraction. Journal of Physiology 243, 1—43. Pandy, M.G., Anderson, F.C., Hull, D.G., 1992. A parameter optimization approach for the optimal control of large-scale musculoskeletal systems. Journal of Biomedical Engineering 114, 450—460.
530
T. Spa¨ gele et al. / Journal of Biomechanics 32 (1999) 521—530
Pandy, M.G., Garner, B.A., Anderson, F.C., 1995. Optimal control of non-ballistic muscular movements: a constraint-based performance criterion for rising from a chair. Journal of Biomedical Engineering 117, 15—26. Pandy, M.G., Zajac, F.E., Sim, E., Levine, W.S., 1990. An optimal control model for maximum-height human jumping. Journal of Biomechanics 23, 1185—1198. Patriarco, A.G., Mann R.W., Simon, S.R., Mansour, J.M., 1981. An evaluation of the approaches of optimization models in the prediction of muscle forces during human gait. Journal of Biomechanics 14, 513—525. Pfeiffer, F., 1989. Einfu¨hrung in die Dynamik. Teubner, B.G., Stuttgart. Schiehlen, W., 1986. Technische Dynamik. Teubner, B.G., Stuttgart. Schiehlen, W., 1991. Computational aspects in multibody system dynamics. Computer methods. Applied Mechanics and Engineering 90, 569—582. Spa¨gele, T., 1998. Modellierung, Simulation und Optimierung menschlicher Bewegungen. Ph.D.-Thesis, Institute A of Mechanics, University of Stuttgart.
Spa¨gele, T., Kistner, A., Gollhofer, A., 1995. An optimal control technique for the optimization of human movements. In: Ha¨kkinen, K., Keskinen, K.L., Komi P.V., Mero, A. (Eds.), Book of Abstracts XVth ISB Congress. Gummerus Printing, Jyva¨skyla¨, Finland, pp. 870—871. Winters, J.M., 1990. Hill-based muscle models: a systems engineering perspective. In: Winters, J.M., Woo, S.L.Y. (Eds.), Multiple Muscle Systems: Biomechanics and Movement Organization. Springer, New York, pp. 69—93. Winters, J.M., Stark, L., 1987. Muscle models: what is gained and what is lost by varying model complexity. Biological Cybernetics 55, 403—420. Zahalak, G.I., Ma, S.-P., 1990. Muscle activation and contraction: constitutive relations based directly on cross-bridge kinetics. Journal of Biomechanical Engineering 112, 52—62. Zajac, F.E., Winters, J.M., 1990. Modelling musculoskeletal movement systems. In: Winters, J.M., Woo, S.L.Y. (Eds.), Multiple Muscle Systems: Biomechanics and Movement Organization. Springer, New York, pp. 121—148.