Walking dynamics from mechanism models to parameter optimization

A vailable online at www.sciencedirect.com Procedia IUTAM 2 (2011) 199–211 A vailable online at www.sciencedirect.com 2210-9838 © 2011 Published by Elsevier Ltd. doi:10.1016/j.piutam.2011.04.020 2011 Symposium on Human Body Dynamics Walking dynamics from mechanism models to parameter optimization Werner Schiehlena,* and Daniel García-Vallejob aInstitute of Engineering and Computational Mechanics, University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany bDepartment of Mechanical and Materials Engineering, University of Seville, Camino de los Descubrimientos S/N, 41092 Seville, Spain Abstract The paper deals with the historical development of human body dynamics, and it presents results received for simple models by parameter optimization. The scientific research on human walking dynamics started already in the 19th century and was later promoted by the physicist and physiologist Otto Fischer who published in 1906 his fundamental book. Fischer used the mechanism theory for modeling and analysis of human walking. His research was based on the barycenters of the corresponding reduced mechanisms. By the end of the 20th century computational multibody dynamics provided more complex models which were applied to human body dynamics, too. More recently parameter optimization has been used to deal with the overactuation of biomechanical systems still a very active research topic in biomechanics. As an example a gait disorder simulation is presented showing that even today mechanism models, muscle group selection, inverse dynamics approaches and parameter optimization techniques using energy and aesthetics criteria are essential tools. © 2011 Published by Elsevier Ltd. Peer-review under responsibility of John McPhee and József Kövecses Keywords: Walking dynamics; mechanism models; parameter optimization 1. Introduction The scientific research on human walking is an interdisciplinary science. It started in physiology, it was investigated in biomechanics and considered later in multibody dynamics. Due to the overactuation * Corresponding author. Tel.: +49-711-685-66388; fax: +49-711-685-66400. E-mail address: [email protected]. Open access under CC BY-NC-ND license. Open access under CC BY-NC-ND license. 200 Werner Schiehlen and Daniel García-Vallejo / Procedia IUTAM 2 (2011) 199–211 by the many muscles inverse dynamics and parameter optimization are used. In this introduction some early books and a few recent contributions representing the historical development are mentioned. Physiology research is devoted to experiments, data collection and model design used for simulations. An early book was published in 1836 by Weber and Weber [1] entitled “Mechanics of the human walking tools – an anatomical-physiological investigation”. In this book, Fig. 1(a), the authors deal with the science of walking and running, position and joints of the parts of the skeletal, and experiments on walking and running. Forty-five years later in 1881 Vierordt [2] published a book entitled “Walking of healthy and sick humans: using auto-recording methods”. In this book the author considered the differences between normal gait and gait disorders. An engineering approach was used in 1906 by Fischer [3] in his book “Theoretical principles for mechanics of the living body – with special applications on the human body and some motions of machines”, Fig. 1(b). Fischer’s approach will be reviewed in more detail in the second chapter. Many research results of the 20th century are summarized in the books of Vaughan et al. [4], Winter [5] and the two review articles by Zajac, Neptune and Kautz [6]. In multibody dynamics, the focus of the research is related to the modelling, the equations of motion of walkers and efficient simulation tools. Two more recent papers are due to Silva and Ambrosio [7] and Chenut, Fisette et al. [8]. For the overactuation problem to be solved by parameter optimization the human muscle energy expenditure or metabolical cost, respectively, is a most valuable criterion presented by Umberger et al. [9]. Also Kim et al. [10] are using optimization techniques considered in the third chapter. The metabolical cost of transport depend on the muscle forces, the muscle activation and the neural excitations. While the muscle forces can be found by inverse dynamics, for the muscle activation and the neural excitation an extended inverse dynamics approach is useful as proposed by Ackermann and Schiehlen [11] and Ackermann [12] also used by Garcia Vallejo and Schiehlen [13] together with parameter optimization. The extended inverse dynamics is based on the model’s equations of motion, and the scalar differential equations of contraction and activation dynamics. More details are presented in fourth chapter. Fig. 1. (a) Title page of Ref. [1]; (b) Title page of Ref. [3] Werner Schiehlen and Daniel García-Vallejo / Procedia IUTAM 2 (2011) 199–211 201 2. Mechanism models Fischer used link systems for the modeling of the human body as shown in Fig. 2. Step by step the local centers of mass of the body parts are assembled to the center of mass of the whole walking human. Fig. 2. Table 1 from Fischer [3] Thus, Fischer identified the center of mass of the human body by means of static principles. Then, he used the approach of barycenters for the evaluation of the equations of motion. A general theory of this approach was published by Wittenburg [14]. Fischer dealt with the swing phase of a leg in detail, too, and he presents the equations of motion for the thighs, the shanks and the feet. By inverse dynamics Fischer provided a formula for the torques generated by the muscles. He found that the swing motion of the leg is not only a pendulum oscillation as suggested by Weber and Weber [1] but also driven by muscles. Fischer stated also that the computed muscle torques do not identify the muscles directly involved in the different phases of the leg motion. He suggested only that most probably Iliopsoas, Rectus Femoris and Tibialis Anterior are initiating the swing phase and their contraction disappear after one third of the swing period. Mechanism models have been applied for controlled biped walking by Beletzky [15], see Fig. 3. His mechanism models were analysed with respect to the limit cycles of the highly nonlinear equations of motion. However, Beletzky used only mechanical principles and did not consider the overactuation problem of the many human muscles. Center of mass of the human body 202 Werner Schiehlen and Daniel García-Vallejo / Procedia IUTAM 2 (2011) 199–211 Fig. 3. Biped walker from Beletzky [15] Today very detailed musculoskeletal models can be obtained by the commercial software AnyBody [16] providing a broad spectrum of applications. There is also available the free community AnyScript [17] to discuss models of the human body. 3. Parameter Optimization The problem of muscle coordination of human walking already stated by Fischer [3] is a central research topic in biomechanics. The review articles by Zajac et al. [6] present a comprehensive introduction to concepts, dynamics, simulations and clinical implications. With respect to optimization and control the contributions by Anderson and Pandy [18, 19] are very essential using global parameter optimization for jumping and human gait. Furthermore, Thelen et al. [20] and Thelen and Anderson [21] proposed the computed muscle control method with local optimizations to produce empirically reasonable muscle activations. Seth and Pandy [22] use the neuromusculoskeletal tracking method which shows a high computational speed up compared to parameter optimization. There are several choices for the optimization criteria of human walkers. Healthy walker usually will walk with minimal energy per distance traveled what is useful for walkers with orthoses or prostheses, too. Elderly or injured walkers may prefer a minimal pain criterium which is related to the joint reaction forces. More general criteria are the least action principle by De Sapio et al. [23] and the min/max criterion by Rasmussen et al. [24]. 4. Gait Disorder Simulation Mechanism models and parameter optimization are important tools for studies of human walking. These approaches are applied to a gait disorder simulation. First of all a spatial mechanism model is introduced, then the inversion of contraction and activation dynamics are added, the constraints of the parameter optimization are considered and a numerical example is presented. The gait disorder is modeled in the simulation by adding a weight of two kilograms to the right foot resulting in an unsymmetrical walking. 4.1. Spatial model The three-dimensional model of the human body used by Garcia-Vallejo and Schiehlen [13] for the parameter optimization approach is composed of 7 moving bodies: a body called HAT representing the Werner Schiehlen and Daniel García-Vallejo / Procedia IUTAM 2 (2011) 199–211 203 pelvis and the trunk, two thighs, two shanks and two feet, see Figure 4. The thighs are connected at the hips to the HAT by spherical joints, the shanks and thighs are connected by revolute joints representing the knees and the foot and shanks are connected by revolute joints representing the ankles. This is a simplification compared to other three-dimensional models that can be found in the literature. However, this simplification allows the derivation of the symbolic equations of motion by software Neweul-M2 [25] for the free flying tree structure composed by the mentioned 7 bodies without any external constraint. Fig. 4. Model of the human body The kinematic tree in Figure 4 is described by the following vector of nc = 16 generalized coordinates: T IIIIII zyx ][ 78671616164534131313111111 EEJEDEEJEDJED q (1) where the subscript I refers to the inertial frame, subscript 1 refers to body HAT, which is in composed of the pelvis and the trunk, subscripts 3 and 6 refer to right and left thighs, respectively, subscripts 4 and 7 refer to right and left shanks, respectively, and subscripts 5 and 8 refer to right and left feet, respectively. When a subscript is written as ij it means a relative motion of body j with respect to body i. Once the kinematic tree representing the skeleton is described, the contact of this tree with the ground is added. The contact conditions in the different walking phases must be included by means of unilateral constraints. However, due to the use of an optimization framework in which it is possible to constrain the normal contact forces to be positive only, the contact with the ground can be modelled using the simple bilateral constraints associated to revolute joints. Therefore, the contact forces can be easily added to the model by using Lagrange multipliers as:      8...,2,1,,   ph ph T ph m rȜCfABqqqqqkqqM  (2) where is the vector of Lagrange multipliers at phase ph of the motion. Note that the previous equation is used together with constraint equations forcing the normal contact forces to be always positive. Moreover, hard impacts will be avoided. ph Ȝ 204 Werner Schiehlen and Daniel García-Vallejo / Procedia IUTAM 2 (2011) 199–211 Furthermore,  qM is the (nc×nc)-mass matrix of the system, q, q  and q  are the (nc×1)-position, velocity and acceleration vectors, respectively,  qqk  , is a (nc×1)-vector describing the generalized Coriolis forces,  qqq  , r is a (nc×1)-vector including generalized gravitational forces, passive generalized moments at the joints due to tissues interacting with the joints and generalized viscous damping torques at the knees and hips. The (nc×1)-vector m fAB includes the generalized forces exerted by the muscles actuating the model and matrix A is the constant (nb×Nm)-matrix of moment arms and is used to calculate the torques generated by all muscles at the actuated joints and matrix B is a (nc×nb)- distribution matrix used to get the generalized torques from the torques at the actuated joints. Finally, the generalized coordinates (1) are parameterized by spline polynomials within the parameter optimization approach where the equations of motion (2) serve as constraints. 4.2. Muscle forces, inversion of contraction and activation dynamics By parameterizing the muscle forces using spline polynomials, the calculation of the time derivative of the muscle force by means of numerical differentiation using finite divided differences can be avoided by the implementation of the analytical derivative of the spline polynomial in the interpolation subroutine. Furthermore, the time derivative of the activation, what is required when inverting the activation dynamics, can be also obtained due to the implementation of the first and second derivatives of the spline polynomials. The inversion of the contraction dynamics is needed to obtain the values of the muscle activation a since they are required to evaluate the energy expenditure according to the procedure proposed by Umberger et al. [9]. Once the activations are obtained, using their time derivative a  it is possible to invert also the activation dynamics so that the neural excitations u are also obtained. The neural excitations are required for two reasons: one is that they are also involved in the calculation of the muscle energy expenditure and the other is that they are involved in some of the nonlinear constraints of the optimization procedure since their values are restricted to the interval [0, 1]. The calculation of the time derivative of the muscle force is carried out by implementing the analytical derivative of the spline polynomial used in the interpolation subroutine. Furthermore, the time derivative of the activation, what is required when inverting the activation dynamics, can also be obtained by the implementation of the first and second derivatives of the corresponding spline polynomials. Figure 5 shows a flow diagram summarizing the inversion process of the contraction and activation dynamics. 4.3. Parameter optimization and constraints The complete set of design variables are summarized in vector Ȥ. This vector is built from four vectors: i) a vector containing all nodal values of the generalized coordinates, ii) a vector containing all nodal values of the different muscle forces, iii) a vector with eight components representing the duration of the eight phases of a walking cycle, and iv) a vector collecting all geometrical parameters describing the kinematic constraints of the feet on the ground. According to the previous explanation, the vector of design variables can be written as: >@ T T g T ph mT j T iptfqȤ (3) where indices i, j and ph are running from 1 to nc, Nm and 8, respectively, and g is just a subscript meaning that parameters in are geometrical. In the three-dimensional model presented before the number of coordinates nc is equal to 16 while the number of muscles Nm is equal to 28. g p Werner Schiehlen and Daniel García-Vallejo / Procedia IUTAM 2 (2011) 199–211 205 Fig. 5. Flow diagram of the inversion of the activation and contraction dynamics It is generally accepted that the central nervous system tries to minimize the metabolical energy consumption during normal walking when dealing with muscles recruitment, specially for walking long distances. For this reason, it makes sense to obtain muscle forces and generalized coordinates by minimizing the metabolical cost of walking. Umberger et al. [9] provided a measure of the metabolical expenditure including thermal and mechanical energy liberation rates during simulated muscle contractions of mammalians at normal body temperature. According to their model, the total energy rate of a single muscle is written as follows:  p,,,,, uafvlEE cecece  (4) where ce l is the contractile element length, ce v is the contractile element velocity, ce f is the contractile element force, a is the muscle activation, u is the neural excitation and p is a vector summarizing all constant muscle parameters required to evaluate the energy rate, see Umberger et al. [9]. A more meaningful measure of energy expenditure for steady-state walking in normal conditions is the energy expended per unit of length. This is called the total energy of transportation and reads as LR t LL E E (5) where E is the time integral of the total energy rate E  over one cycle, and R L and L L are the right and the left steps walked in the simulated cycle. Inverse Dynamics E q uations of motion 206 Werner Schiehlen and Daniel García-Vallejo / Procedia IUTAM 2 (2011) 199–211 Since the time histories of the muscle forces and of the generalized coordinates are obtained by optimization techniques trying to minimize the energy consumption there is a need to follow a certain motion pattern. This fact has some advantages in the case of prosthesis design. For example, the simulated motion of an individual wearing a prosthesis will be close to normal walking patterns because that is desirable for aesthetical reasons. Also, the simulated contact forces will be close to those of a normal walking cycle to avoid significant modification of the contact forces at the non-damaged foot. The deviation with respect to normal walking patterns is evaluated as follows:    ³¦  fx t t n ii m ii dev dt txtx J 01 2 2 V (6) where i x is a time dependent variable of the model and m i x refers to the experimentally measured value of the same variable. These variables i x with x ni ...,2,1 , include the generalized coordinates and ground reaction forces. In (6), i V is a characteristic measure of the time variability of i x. Dividing by 2 i V , one can scale the differences between measured and simulated values of all i x. In this investigation, the mean square deviation with respect to the averaged mean is used as i V . The measured motion used in this research was obtained by Ackermann and Gros [26] by measuring the walking motion of a subject wearing sport shoes and walking at his preferred velocity. The final value of the cost function is calculated using the metabolical cost of transportation, t E, and the measure of the deviation from normal walking patterns, dev J, as follows: devJ t EJEf ZZ  100/ (7) where t E is divided by the factor 100 for balancing of the two terms of the cost function (7) to get comparable numbers, and E Z and J Z are two weighting factors. The procedure suggested by Ackermann [12] and used in this research, avoids the expensive forward time integration by parameterization of the time histories of the generalized coordinates. The parameterized solution is found by searching for optimum parameters at certain node positions. Since walking is a periodic motion, other authors have also used Fourier series to parameterize the motion, i.e. Peasgood et al. [27]. In fact, the time derivatives of the parameterized function can be easily implemented, avoiding the numerical differentiation used by Ackermann [12]. In addition, the interpolation can be split into two parts: a more computationally expensive one that can be done in a preprocessing stage and the other that is done during the optimization. The solution of the optimization algorithm must fulfill a set of constraints as stated in the introduction of this section. In principle, the set of constraints is similar to the one used in [12] and is summarized as follows: 1. Neural excitations must be bounded in the interval [0, 1]. This kind of constraint ensures that muscle forces are consistent with the activation and contraction dynamics of the muscles. 2. Ground clearance must be positive or equal to zero to avoid penetration of the feet into the ground. 3. Positive normal contact forces to model unilateral constraints between the feet and the ground. 4. The feet tangent contact forces must be consistent with Coulomb’s friction model to avoid foot sliding. 5. The averaged velocity is fixed. 6. Design variables are bounded. These bounds may be due to some physiological reasons like for example the amplitude of the relative motion allowed by a certain joint. Werner Schiehlen and Daniel García-Vallejo / Procedia IUTAM 2 (2011) 199–211 207 7. Other physiological constraints that may help to the convergence of the optimization algorithm like for instance constraining the maximal achieved knee flexion during the swing phase or the maximal achieved hip extension during the stance phase. 8. Equations of motion must be fulfilled according to a certain tolerance. 9. Kinematic constraints must be fulfilled within a certain tolerance. 4.4. Numerical example of a gait disorder simulation In this section, the three-dimensional non-symmetrical (3DNS) model is studied, see Figure 4. In order to keep the computation time as small as possible 18 nodes are used to parameterize each time dependent variable. However, the optimization algorithm used took up to 65 hours in a 64-bit quad core processor with a clock speed of 2.40 GHz to find a solution within the tolerances of TolX = 10í4 for the solution Ȥ, TolFun = 10í2 for the objective function, and TolCon = 10í2 for the constraint violation, what requires the scaling of constraints and design variables to ensure the enough accuracy of the constraint satisfaction. For comparison, the corresponding two-dimensional symmetrical (2DS) required only 33 minutes to find the solution. In this example, the measured motion as well as the measured ground reaction forces has been used as a reference motion (or normal walking pattern). Two 18 nodes 3DNS models are considered with one of them having a weight of 2 kg attached to the right foot. The attachment of a weight to the ankle results in a higher energy expenditure during walking. The metabolical cost of transportation of the model without the attached weight is 569.56 J/m while 631.14 J/m are obtained when the weight is attached. On the other hand, the deviation from normal walking patterns increased from 3.06 to 3.40 when the weight is attached. Figure 6 shows the time histories of the three components of the ground reaction forces at both feet. The similarity between the solutions with and without the attached weight is remarkable. The solutions are almost coinciding because of the ability of the three-dimensional non-symmetrical model to accept minor variations of its variables to absorb the effect of one-sided disorders. Enforcing the model to follow the measured symmetrical motion without disorder makes the simulated motion to maintain the same level for the normal contact force by adapting the time histories of the generalized coordinates accordingly. Changing the weighting factors in Equation (7) to give less importance to the motion to be followed this effect may be reduced. Some neural excitation time histories of selected right leg muscles of the 18 nodes 3DNS model are shown in Figure 7 and denoted by Table 1. It can be observed a higher neural excitation level of certain muscles between 60% and 70% of the normalized cycle time coinciding with the increase in the maximum dorsiflexion at the right ankle. As shown in this chapter the 3DNS model provides rich information of the dynamics of human gait with one-sided disorders. The model could also be improved by adding extra degrees of freedom to better represent the upper limb motion offering new research directions.