Home

Modelling Martial Arts Techniques

1. L 1 1 5 2 25 3 3 5 4 4 5 5 55 6 migration a Best objective per migration 0 16 T T T 0 14 0 12 Objective o 2 oe 0 06 0 02 i i i i t 1 1 5 2 25 3 3 5 4 45 5 5 5 Migration b Figure 4 6 Best result of experiment 1 continued a The evolution of the best solutions of each regions Here the approximated objective vs migration is plotted The migration method we used is the universal migration method With this method each region will choose another region at random and replace the worst solution of the chosen region with its own best solution during migra tion b The evolution of best objective of the whole population all regions Because of elitism best solutions are saved for next generations we are guar anteed of a monotonous decrease for the objective The best found objective is around 0 0224 62 Activation of best solution l S SS ee SSS Ss SS ogf d o8 h J 07 H 0 6 J B05 H 0 4 H 0 3 H 02 H 0 1 H Hip angle vs time 12 I T T T q deg a T f 0 AR I aa I L L 1 1 L L 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 time s b Figure 4 7 Best result of experiment 2 a The found activation function which is a bit mysterious b The 0 vs time belonging to the best solution It looks
2. Algorithm 1 Algorithm for Hill s muscle model A stabilization mass m is used to keep track of the length of CE The position of the mass is given with x and tells us what the length of CE is In the algorithm f ce and fv ce functions are the normalized force length and force velocity relation ships with respect to Fmaz Furthermore faf is normalized to 0 1 If we want to know the real Fog we should scale fi cr foce af With Fmaz We used Euler method to solve x and ay for easier understanding 3 5 2 The Virtual Muscle model Hill s model only has good prediction for maximal activation Most human motions occur at sub mazimal activation levels Virtual Muscle is developed by the Alfred E Mann Institute California and its purpose is to provide an accurate muscle model Virtual muscle is a Hill s based muscle model where the force characteristics of the different elements are obtained through many animal experiments The results are different characteristics for the elements Especially the active state a has great differences with the ay of Hill s model Observed phenomena are rise time fall time sag yielding used to describe ay rather than chemical reactions Virtual muscle model is also supposed to give good force prediction at sub maximal activation levels To do this the muscle group is divided into multiple units In algorithm 2 we see the algorithm for Virtual Muscle For a precise formulation of the
3. 227 2 Simulation see 2 ence dodo ake e Gee SO Se be 2 8 Skeletal data 2 44 4 4 4 eA Ais ee ae See Be 2 8 1 Free falling thigh with realistic hip 20 2 9 Conclusion ey wis aod S BA EA ee RA eS 3 Building skeletal muscle models 3 1 Structure of skeletal muscles 000 3 2 Activation of the skeletal muscles 2 3 2 1 The cross bridge theory 00 0 3 3 Brief overview of current models 0 3 4 Muscle experiments 2 00000 ee eee 3 0 Virtual Muscle 22 25 4 a e ae ee ew ed ee SS ee eee 3 5 1 Hill s muscle model 00484 3 5 2 The Virtual Muscle model 3 5 3 The model s variables 0 00 0 3 6 Using Virtual Muscle 3 1 5 000 3 6 1 Building the fibers 3 6 2 Building the muscles 0 00 0 3 7 Coupling the skeleton and mtc models iO Eperme n fe Ae ce gs Bes Ales nese ee ee tao ae ena el as a 3 9 Project Mae Geri 2 2 0 0 0 00 0 02 0 ee 3 LO COnClUSIONs 0s lt b shoe aerate te test ake dP dh Sele ele hese 4 Optimizations with Evolutionary Algorithms 4 1 Optimizing our motions 2 0008 4 2 Evolutionary Algorithms EAs 0 4 3 Terminology ai oid Be ee eg A a ee d 4 4 The basic Evolutionary Algorithm 4 5 Extending the basic EA
4. 4 2 Evolutionary Algorithms EAs Evolutionary Algorithms look at a number of possible solutions for a given problem The collection of these solutions is called a population The population dX dt f X U X Y gt F Y Figure 4 1 Schematic overview of optimizing our system For a given input act in our case we get an output Y which we will rate with F Our optimization problems can be seen as finding act such that h act is optimal 52 will evolve over time until the problem is solved or the maximum processing time is met So at each time step we have a different generation supposedly to fit better to the problem than its predecessors An evolution step is done by applying several evolutionary operators on it These operators mimic the behavior of nature s evolution processes The main operators are Selection Select the better solutions from the current generation Whether a given solution is better than another is determined by a user input function The objective function Recombination Solutions are created by recombining the selected members from previous generations Mutation Complex problems tend to have many local optima With mutation we want to escape from a local optimum Mutation is done be changing some solutions slightly These operators make that evolutionary algorithms differ from a simple ran dom search Good information are passed through next generations For each operator we have diffe
5. 2 3 Simplifications In modelling it is always a tradeoff between complexity and realism The more realistic we want the model the more complex it will be Below are the simpli fications for our skeletal models 2D motions Only planar motions are considered A main reason for this is the absence of data needed for 3D simulation Rigid bodies We divide the human body into several segments and treat them as rigid bodies The segments are connected with each other via ro tational joints Skeletal models are therefore known as articulated models Rigid bodies are idealized physical objects and their mechanics are well understood Part of the body we will only model and simulate parts of the body that undergo motions A limb for example There reason is that we need to build muscle models as well The fewer bodies the fewer muscle models are needed Fixed base and loose end Our skeletal model can now be seen as a chain of rigid bodies It has a base point and an end point The base point is closer to the whole body than the end point We fix the base and let the end loose By doing this our skeletal models have only rotational degrees of freedom Simplified joints the modelled joints allow pure rotational motions within a certain range With real life joints we have gliding and sliding 5 With the above simplifications we can build a leg or arm model A martial arts technique which can be simulated with a 2D leg model is the front
6. It becomes clear that some types of experiments cannot be done on human mus cles Therefore animals are used When performing such experiments scientists try to keep the number of animals used as low as possible 32 3 5 Virtual Muscle Virtual Muscle 3 1 5 is a program that builds muscle tendon complexes and describes their behavior It is written in Matlab and the muscle models can be exported as Simulink blocks In the Simulink environment these block can be fully integrated with other Simulink blocks The equations that Virtual Muscle uses to describe the behavior of skeletal muscles are determined by best fit procedures on data that come from many experiments These experiments studied the muscle behavior under a wide range of physiological conditions For people who are interested in how this is precisely done we refer to 16 Hopefully the behavior can be extrapolated and we get realistic results 3 5 1 Hill s muscle model The Virtual Muscle Model is a Hill s based muscle model and is categorized as a lumped muscle model It treats the fascicle as a scaled version of a sarcom ere It assumes that all sarcomeres behave uniformly within a muscle fiber and that all fibers within a fascicle also behave the same System of engineering alike elements are then used to represent the muscle tendon complex mic To understand how Virtual Muscle works it is better to look at the more simple Hill muscle model first In Hill s muscle model
7. Unlike in mechanics we do not have theories that describes the complex muscle behavior well We have to rely on muscle experiments Once these experiments are done we can formulate equations and specify coefficients that describe the muscle s behavior The most well known muscle experiments are listed below 31 isometric The length of the muscle is kept constant and the muscle is maxi mally activated The maximal produced force is that length is called the tetanic force isotonic A known constant force is put on the muscle and quantities of interest are measured isokinetic The contraction speed is kept constant and relevant quantities are measured To get accurate data from the muscle of interest it should be isolated from other muscles Dependent on which state the experiment is done the muscle is cut away from its neighboring muscles The experimented muscle can be in different states in vivo Latin for within in the living The experiment is done on the living tissue of a whole living organism ex vivo Latin for out of the living The tissue is still alive but outside in the organism in vitro Latin for with in the glass These experiments are done in a con trolled environment where the biological structure is separated from its owner in situ Latin for in the place The experiment is done on the normal oper ation location of the tissue The owner could be dead or alive
8. except for the case where we use act 0 as input Furthermore we do not know what activation function we should feed the model It is amazing how the human mind can do all kind of amazing motions with these complex force producers We should however not forget that we needed years of training to do simple movements like walking for example From a mathematical point of view the behavior of our musculoskeletal system is very nonlinear Optimizing the output is therefore a real challenge In the next chapter we will see how we can use evolutionary algorithms to perform optimizations 51 Chapter 4 Optimizations with Evolutionary Algorithms 4 1 Optimizing our motions We have built a musculoskeletal model and can perform all kinds of experiments with it The input is an activation function act for each mtc model The output vector Y contains variables we are interested in Many times we want to find an optimal output Examples fastest kick economic kick etc To do this we should first define an objective function F for Y It is obvious that depending on one s interest many objective functions can be defined Then we introduce a function h that maps act to F Y h is a very complex function and with each simulation we can only know one point of h Our optimization problems become as follow find act such that h is optimal Figure 4 1 makes it all more clear In this chapter we will use Evolutionary Algorithms to find the optimal act
9. 1215 84 bad Muscle 4 Vastii VU 4185 93 Muscle 5 Gatrocnemius VU 733 55 391 44 Muscle 6 Soeusvu 14667 s55 i 2454 30 Muscle 7 SoleusvU test 14667 55 30 Muscle 8 lf 0 Muscle 9 0 Muscle 10 0 Comments Next 10 Muscles Figure 3 8 With the program BuildMuscles we can build and export mtc s as Simulink block after we have filled in necessary data Furthermore it can divide the whole muscle into different units We should fill in the fiber distribution is the whole muscle and the number of units we want The higher the number of units the more computation intensive the simulation will be The fiber type distribution are taken from the Yamaguchi 12 table Once the number of units are filled in we can export the mtc models as Simulink blocks In these blocks the muscle equations are implemented In A we see how these equations are implemented The equations are given in 10 and we verified the Simulink block with them During our verification we discovered a small error in the Simulink blocks Two force signals were switched see figure 3 9 Although it was very easy to fix it manually we fixed the source code of BuildMuscles The reason is that we need to build more mtc models and fixing them one by one costs time 3 7 Coupling the skeleton and mtc models The skeleton model outputs joint coordinates and needs net
10. 2 relTol 12 3 relTol 1e 6 relTol 1e 9 10 16 10 14 10 12 10 1 a degrees 10 08 10 06 10 04 10 02 0 9985 0 999 0 9995 1 time t s Figure 4 10 Plot of 0 belonging to the best solution vs time t solved with different relative tolerances relTol 1e 1e 3 1e le ODE45 solver Here we plotted the calculated 0 instead of the interpolated 6 We have zoomed in at t 1 and the usage of variable time steps is clearly visible We assume the higher the relative tolerance the higher the accuracy From the plot we can see that solving with a relative tolerance of 1076 is pretty accurate At the endpoint t 1 the error between the found 0 with relTol le 6 and le 9 equals 1 4e 9 67 pulses 010 1 act 0 1 t s Figure 4 11 Triangular shaped pulses instead of constant pulses A triangular pulse has less activation than a constant pulse the model Things change if we use other type of activation function that are finer Solutions with higher frequencies or recorded EMG data for example If the taken time steps are too large pieces of the activation can be missed by the model By specifying the maximal time step manually we can prevent this The price for this is more computation time however 4 7 4 Experiment 3 another type of activation function For our activation functions we used piecewise constant functions that we call pulses What if the pulses are not constant but of t
11. Chapter 2 Building skeletal models 2 1 Introduction In this chapter we will build skeletal models First we take a look how we model the skeletal system as a mechanical system Then we will study the dynamics of such a system To get familiar with some basic mechanics we derive the equations of motions for a simple thigh model After that we will use SimMechanics to build skeletal models 2 2 The skeletal model The skeletal model is a sub system of our musculoskeletal system It will receive input from the muscle models and produce motion Moving objects on earth obey certain laws These laws are known as Newton s laws of motion and have already been formulated in 1687 by Isaac Newton 3 They are 1 Every object continues in its state of rest or of uniform motion in a straight line unless it is compelled to change that state by forces impressed upon it 2 The change of motion is proportional to the motive force impressed and is made in the direction of the line in which that force is impressed 3 To every action there is always imposed an equal reaction or the mutual actions of two bodies upon each other are always equal and directed to contrary parts Newton s laws are only valid in an inertial frame An inertial frame is a coordi nate frame that does not undergo acceleration nor rotation Real inertial frames do not exist but for many applications we can use a frame that is attached to the ground as an inertial frame
12. Each simulation produces a certain motion A martial arts technique is a optimal motion for a given goal We thus have to optimize the model Muscles pull Skeleton motion Articulated Skeleton model MTC model I Multi body system SE PE Lumped models Figure 1 1 Schematic overview of how we try to find martial arts techniques We model the human locomotor system that has neural activation as input and will produce motion A martial arts technique is an optimal motion with respect to a certain goal If we want to find an optimal motion we should find the neural activation that produces this optimal motion The locomotor system consists of 2 very different subsystems the muscular and skeletal system The muscular system is made up with a number of skeletal muscles which are actually muscle tendon complexes mtcs They are activated by the neural activation As result of the neural activation the skeletal muscles will contract and produce a pulling force This will cause the skeletal system to change and we get motion We will model the skeletal muscles from a mechanical point of view and these type of models are known as lumped models The skeletal system will be modelled with a system of rigid bodies that can only rotate with respect to each other articulated models Because the skeletal model is a system of multiple bodies it can be considered as a multi body system 1 3 Modelling dynamic systems With modelling we mean math
13. J Cheng Ian E Brown Gerald E Loeb Virtual muscle a com putational approach to understanding the effects of muscle properties on motor control Department of Biomedical Engineering Alfred E Mann In stitute for Biomedical Engineering Southern California USA June 2000 Ernest Cheng Ian Brown Jerry Loeb Virtual Muscle 3 1 5 Muscle Model for MATLAB User s Manual Feb 21 2001 Jack M Winters Savio L Y Woo Multiple Muscle Systems Biomechanics and Movement Organization Springer Verslag New York 1990 Stuart I Fox Pierce College Human Physiology Sth edition McGraw Hill Science Engineering Math 2003 85 14 15 16 Hartmut Pohlheim GEATbx Introduction Evolutionary Algo rithms Overview Methods and Operators version 8 7 2005 http www geatbx com Agamemnon Despopoulos Stefan Silbernagl Color Atlas of Physiology 5th edition Thieme Medical Publishers Georg Thieme Verlag New York 2003 http www thieme com Ian Edward Brown Measured and Modeled Properties of Mammalian Skeletal Muscle Queen s University Kingston Ontario Canada 1999 86
14. by the muscles These netjoint moments depends on the stimulation but also on the state of the leg Hip fixed Machine Environment Input net Moments for the joint generated by the muscles 100 Figure 3 15 The top level view of our leg model b The skeletal sub system which is an extended version of the thigh skeletal model 48 El sim_maegeri_VU musc_maegeri File Edit View Simulation Format Tools Help DSHS teels Po S ps Nom Ranen Rame musc_maegeri muso_maegeri isthe muscular part of the maegeri model Given the stimulation and state variable q as input this subsystem will return the net joint moments k cale_joint_moments mt caloMTClength mte_mom_arm jadeg q deg Embedded MATLAB Function2 Embedded MATLAB Function U UE ial Selector Scopet 1100 Figure 3 16 The musculo system of the leg model The mtc lengths and the moment arms should be calculated to allow interactions with the skeletal model 49 a sim_maegeri_VU musc_maegeri participating MICs Eile Edit View Simulation Format Tools Help Dees Bles rl gt m 05 Normal JRH Deal RE The participating muscle tendon complexes MTCs in the maegeri model Input activation and length of the MTCs Output force of the MTCs Activation MT path length cm Biceps Femoris VU Force N Activation MT path length cm Glut
15. cm can be expressed in known terms by differentiating equation 2 11 twice l tom p cos 06 sin 00 2 14 ticm p cos 66 sin 06 14 Figure 2 4 A new angle definition The thigh can be seen as a pendulum driven by gravity When the thigh is held and released at angle a it will swing back and forth from a to a F thus depends on 6 6 and 6 Our system can be seen as a simple pendulum driven by gravity The behavior of the simple pendulum is well known Before we describe the motion we use a different angle definition We name a new angle coordinate a and is defined to be the angle between the negative y axis and the thigh see figure 2 4 When we release the thigh at angle a amaz it will swing back and forth towards the negative direction of the y axis Because of conservation of energy the motion will take place between amar and Qmaz Special situations occur when we release from a 0 degrees or a 180 degrees The system is then in static equilibrium and will not move 2 5 1 Implementation with Matlab The free falling thigh can be very easily implemented in Matlab Our state vector X 0 0 and its time derivative 0 f 6 pce 5 aa 2 15 Ihip To start simulation we should specify the simulation time T and the initial conditions Xo Implementing the free falling thigh with Matlab is fairly easy Below is the Matlab code options odeset RelTol 1e 6 relative tolerance T X ode45
16. de coordinate of that axis Euler s equation can be used to describe the rotational motion with respect to CM bp gt M lome 2 4 Here X M is the contribution of all force moments and torques with respect to CM The force moment is given with M rFrx F 2 5 where F is the distance vector between CM and F s point of exertion The perpendicular component of 7 on F s line of action is also known as the moment arm of force F The larger the moment arm the larger the effect of F We speak of a torque T if two equally but opposite forces F and F act on the body with distance If d is the distance vector between both forces we have Ta dxF dx F 2 6 The perpendicular distance between those forces is known as the moment arm The magnitude of T is given with T moment arm x F In musculoskeletal models we should speak about torques rather than moments When a muscle that originates at body A exerts a force Fm on an adjoined body B body B will react with an equal but opposite force Fm Newton s third law The joint rotation can be calculated with Euler s equation The knee cap and the heel bone are thought to enhance the torque contribution Figure 2 2 makes it more clear The equations of motion for a free moving rigid in 2D can now be written as SO Fy m 0 0 x SF 0 m_ 0 2 7 SM 6 20 in 6 If we name the coordinate vector x y 0 we can write 2 7 in a more compact form we have M
17. depends on the shape of the body If we have a full symmetric body and the body is made of uniform material then its CM is located in the geometric center for example this happens for an iron ball The inertia moment I is usually determined with respect to the CM of the rigid body and is given by Tow fram 2 1 Here r is the perpendicular distance between an element of mass and point CM Icm tells us about the amount of resistance against rotational motion around CM In many applications the rotation is around another point If we name the point O and the distance between O and CM equals r we can calculate Io with Steiner s parallel axis theorem Io Icm mr 2 2 Note that Euler s equation see next section is only valid if O is fixed This restriction does not hold for CM Because of the quadratic contribution of r the shape of the body plays a very important role In practice we can experience this For example it is much easier to raise a chambered leg than a fully extended leg 11 2 4 1 Equations of motion An unconstrained rigid body in 2D has 3 degrees of freedom 2 translational and 1 rotational General motions of the rigid body can be seen as combinations of translational and rotational motions The dynamics of the translational motions can be written with Newton s equation Newton s equation of motion along an axis is x STen 2 3 With gt gt F the sum of all forces that acts along that axis and x is
18. kick With this technique we can kick a target in front of ourselves 2 4 Dynamics of a rigid body in 2D Before we start to build a skeletal model we will look at the dynamics of a rigid body Dynamics studies the relation between motions and forces A rigid body is an idealized physical object It can be seen as a system of particles where their relative distances are fixed Real rigid bodies do not exist but one might think of a brick or iron bar The skeletal bones can be considered rigid as well If we want to study the planar motion of a rigid body we need to know the following properties of the rigid body first e mass m e inertia moment Jo with respect to a certain point O e length l of the body e center of mass CM by specifying p or d where p length between the CM and the proximal end d length between the CM and the distal end We have p d l 1 The rigid bodies that models body segments includes the bone s and other organic ma terials around it such as muscles and skin 10 Figure 2 1 The thigh in an inertial frame Because the thigh is heavier on the proximal side p is smaller than d The configuration of the thigh is determined by 0 In figure 2 1 we have a sketch of a rigid body in an inertial frame The CM is a very special point The gravitational force F mg is modelled to act on this point 3 The position of CM is denoted with rom xom yom where xdm TCM aap ydm YCM Te The CM
19. know that X ty At Xngi Xn X At O At 1 3 Euler method uses the first order Taylor series to approximate X tn 1 Ke AX ING 1 4 The error is thus of order O At and this method is not very accurate A more accurate method is the midpoint method or Euler Heun and goes as follow ky A X n tn 1 5 gt gt gt 3 1 ko f X t 5 At Unptitn 54 1 6 Xay Xn k AER OE 1 7 With this we have an error term of O At and should be more accurate Euler s method and the midpoint method belong to the family of Runge Kutta methods They are also known as the first order Runge Kutta and second order Runge Kutta methods respectively Runge Kutta methods use different approximations for X and combine them in such a way that a higher order error term is obtained Higher order Runge Kutta methods use more function evaluations than lower order but should be more accurate A common way to bring down the computation costs is to use variable time steps Numerical integrators can be categorized into two main methods implicit and explicit Explicit methods use known state variables from the current state to calculate state variables for the next state while implicit methods do not We will use Euler s method to show the difference The explicit Euler method is known as Euler forward method and calculates the next state as follows Xe SN ee NG 1 8 The implicit form of Euler method is known as Euler Backward metho
20. produce output at fixed time points that are equally distributed along the time axis with distance dt The objective function is then approximated with N 1 Fw dt XO 0 t 10 4 2 i 0 with to 0 75 and ty 1 2With variable time step solvers linear interpolation is used to produce results at specified time points This will result in extra approximation errors 56 Volts Oo N 0 0 o s 1 0 1 5 2 0 2 5 3 0 Time Seconds Figure 4 2 A typical EMG signal Voltage vs time The electrical activity in the muscles The amplitude tells us how active the muscle is and can be related to the neural activity If we want to use this signal we should clip the signal to 0 1 Encoding the activation function Although the objective function is very important formulating it was fairly easy The next step is encode the solution which is the activation function This is a difficult step because the activation function is a very complex function It stands for the normalized neural activation function If we want to do it right we should model the neural activation and parameterize this function A solution can then be seen as an array of the parameters that determine the activation function According to the manual of Virtual Muscle 3 1 5 11 we can use recorded electromyography EMG as activation EMG measures the sum of all synchronous currents in each muscle fiber In figure 4 2 we see
21. ranking The sp should be in the range of 1 0 2 0 Roulette nmates sp size size is the tournament size Linear ranking With tournament and roulette selection fitter members have more chance to be selected These selection chances are proportional to their fitness evaluation Sometimes we have too dominant solutions Their proportion is way bigger and it is very likely that weaker members have practically no chances to be selected Think of 1 percent or less Linear ranking gives weaker members relatively fair chances to be chosen too This will guarantee the diversity The sp parameter regulates the rate of diversity Elitism With elitism we keep a number of best solutions of the current generation for the next generation These are called the elite members By doing this best solution of each generation will be better or equal to the best solution of the previous generation Typical numbers of elite members are 1 or 2 C 2 3 Recombination The implemented recombination operators are discrete Two children are created from two parents with the use of 2 random index arrays 83 uniform Same as with discrete but with one random index array crossover m m is the number of cross over points If m is small eg 1 or 2 the children will inherit much of their parent s look This of course depends on how the solution is encoded C 2 4 Mutation We have only implemented one type of mutation and that is mutation wi
22. techniques With 8 we can calculate some average skeletal parameters This source only provide skeletal data however For our musculoskeletal model we need muscle data as well These muscle data should fit with the skeletal data and we thus have a problem Luckily the faculty of Human Movement Sciences of de Vrij Univer siteit of Amsterdam has provided us some realistic musculoskeletal data The skeletal data are listed in 2 1 and table 2 2 2 8 1 Free falling thigh with realistic hip In real life the body joints and thus also the hip allow a limited rotation range Ligaments and anatomical structure of the skeleton are responsible for this We model this by adding stiffness and damping to the joint if the joint angle exceeds one of the two extreme angles The damping and stiffness will produce a counter moment that is proportional to Oextreme The proportion is given with damping constant k and stiffness constant k respectively We have damping il if 2 16 k 0 Tii Omax if 0 gt Omax k2 0 Omin Lf OX Gr n stif fness 2 if 2 17 k2 0 Omax if 0 gt Omax For the hip we have 23 El ushi ro_skel_customjoint File Edit View Simulation Format Tools Help Discus IQS gt mfo Normal Fane iJ Free falling thigh with realistic joint Models the free falling thigh with a realistic hip joint i h pe Em Machine Environment 100 Figure 2 8 The new free falling thigh model Here
23. the tendon is modelled with a passive elastic force element whilst the muscle group is modelled with a passive and an active force element The force of the active element CE depends on its length Dog and neural activation act The force of the passive elements PE and SE depends on their lengths Where pg lcg and lsg lor lce and lor is the mtc length In figure 3 4 we see a scheme of the Hill s model The characteristics of the force producing elements are Passive element SE Fsg lor lce kgp L Lestack if L ee L stack half parabola ksp is the stiffness constant of the tendon Experiments have shown that the tendon is very adaptive We have Isz stack 0 96L 5 opt where LsE opt is the length of the tendon when maximal muscle force is applied to the tendon Passive element PE Fpg lce kpg L mi Lsiack if L 2 Lstack half parabola kpg is the stiffness constant and is determined experimentally Active element CE this the most challenging element and 3 relations are used to model the behavior we have Isometric Force length Fi icz Fmax 1 Z225 Eort 2 This is a parabolic curve and W is a measure for the width of the parabola This relation agrees with the cross bridge theory Isokinetic Force velocity F icg Enmggb av This relation is noticed by Hill in 1938 and the constants a and b are fiber type dependent Active State as ay stim t The active state can be seen as the concentration of boun
24. uniformly distributed and therefore the force scales linearly with the amount of overlap Change in length of the sarcomere the rate of forming cross bridges de pends on wether the sarcomere lengthens or shortens Active state of the sarcomere binding sites for the cross bridges are only freed if the concentration Cat exceeds a certain threshold During activation calcium ions are released When activation stops the calcium ions are actively pumped back 3 3 Brief overview of current models Current muscle models can be divided into 3 categories Namely Huxley This model keeps track of the state of a finite number of cross bridges Furthermore the transition of these states over time is described by a set of equations It needs lot of computation power and lots of equation coefficients Theoretically this type of modelling has the highest accuracy Distribution Model DM In order to decrease the computational effort the DM model groups the individual cross bridges into several groups The behavior of the groups is described Lumped Element Model These models are composed of force producing el ements from the world of engineering The elements can be connected to each other either parallel or serial Hill s muscle model is the most well known lumped element model Most muscle models currently used are of this category The muscle model that we will use falls in the category Lumped Element Model 3 4 Muscle experiments
25. used to SimMechanics uses blocks that represent idealized physical objects like rigid bodies revolute joints etc With these blocks we can build fairly complex phys ical systems 6 SimMechanics can automatically formulate the equations of motions SimMechanics requires Simulink and Matlab to be installed 2 7 Free falling thigh with SimMechanics 2 7 1 Building the model In section 2 5 1 we have modelled a free falling thigh with the Matlab code We will use SimMechanics to model the free falling thigh again with same data The steps are explained one by one 1 Specify the machine environment With the environment block we model the environment for our system Gravity the motion s dimension 2D and analysis mode forward dynamics are specified here Error tolerances for several cases can also be set We use the default values 2 Choose an inertial frame The ground is the ever present body in our models We attach a frame to the ground or some object that is fixed to the ground This frame is our inertial frame We as observers are particulary interested in the output motion and configuration expressed in this frame The origin of the frame is set to 0 0 0 and it coincides with the hip Because the hip is fixed to the trunk we name our ground block trunk 3 Add a revolute joint If we want to add another body to an existing body we should use joints If we add an airborne body in space we can use a joint with 6 DoF We want t
26. various relations we refer to 10 35 Effective length Rise and fall time Sag Yield Parallel elastic Thick flament Foroe tength ae er eo element compression Fpe Fpg1 Af FpE2 Eeer Total parallel pe contractile elastic force force FTota F petFce mass D Total muscle force Series elastic Total contractile element force element Activation force relationship K Figure 3 5 Diagram of relations that Virtual Muscle uses to model a muscle unit Only the diagram of one unit is given A muscle unit has multiple state variables that should be solved The more units are modelled the more state variables are needed The sum of the forces of all units forms the total contractile force This force will be transferred to the skeleton via the SE tendon element A mass is used to keep track of the length of the muscle group Input Neural activation function act t and mtc length loz t Output Tendon force Fsg for t 0 dt T do Fse fse lor 2 stim freeuit act t Emus D Funit il stimi where Funiti fpealx forlz stimi wa Fse Fmuscle m t e dz x x dtt end Algorithm 2 Algorithm that Virtual Muscle uses if more units are mod elled and sub maximal activation is considered During sub maximal acti vation not all muscle units or fascicles are recruited The function frecruit turns the activation into different stimulation frequencies for the different units fc
27. xprime 0 10 0 0 options numerical integrator with simulation time and initial conditions function dx xprime t x g 9 81 gravitational constant m s 2 m 16 944 mass kg 15 p 0 2100 f m g mom f p cos x 1 I_cm 0 4180 I_hip I_cm m p p dx x 2 mom I_hip proximal distance m gravitational force N moment around hip Nm x 1 theta inertia moment wrt CM k m 2 Inertia wrt hip and not Hip_cm k m 2 time derivative of state vector In figure 2 5 we see the results of a simulation for Xo 0 0 6 0 0 0 7 and T 10 seconds The results agree with our expectations 16 2 5 2 Multi body dynamics We have successfully simulated the free falling thigh with Matlab For more realistic motions we need more bodies however and we enter the field of multi body dynamics Multi body dynamics is well studied and can be divided into two areas forward and inverse dynamics In inverse dynamics the motions are known and we want to know the forces that cause these motions In forward dynamics the forces are known and we want to study the motions that are caused by these forces Our problem can thus be considered as a forward dynamics problem Currently there are several software packages specialized in solving multi body dynamics problems There is even a chip PhysX that is specialized in physical computations www ageia com This chip is mainly used to make game
28. 000 4 5 1 Multipleregions 000 4 5 2 Different strategies 0 4 4 6 Project Static Thigh iuao p e DEDERE RAA AA A AOI WS gp EAS aporia a dee a e E R E AA a 4 7 The experiments aoaaa eee eee 4 7 1 Experiment 1 f 400 ose ee ae ee pee 4 7 2 Experiment 2 f 80 scene ea ee 4 7 3 How reliable is our best solution 002 4 7 4 Experiment 3 another type of activation function ArS Conclusion fea a a Pe ee oe A e Appendices A How does Virtual Muscle work B Grieve s method B 1 Calculating the mtc length with Grieve constants B 2 Calculating the moment arm with Grieve constants The GA operators C 1 Population Parameters o oo 0000000 0s C 1 1 Region Parameters 2 500004 C 2 Initializing the population 0 0 C 2 1 Terminal constraints 000008 C22 Selection ae ao nt Soden oS de kk ee ee A C 2 3 Recombination 0000008 C2 Mutation 62 5 a A Ae Set a ee att Chapter 1 Introduction We want to let the computer find optimal movement as from martial art tech niques We should model human motions and optimize these motions It in volves mathematical modelling programming and optimization In this chapter we look how we can formulate and solve these problems in the most general sense 1 1 Martial arts Since the beginning of mankind people have to protect themselves against other
29. 2 8 F and M are known as the generalized force vector and generalized mass matrix respectively In our case F is known and our problem becomes q M F 2 9 This type of problem is known as forward dynamics We study the motion of a system as a result of known forces acting on it The other way around is known 12 Figure 2 2 A torque example When force F acts on a body A that is con nected to another body via another body B body B will react with equal but opposite force Fn We then speak about a torque arm is the perpendicular distance between those forces and is known as the moment arm The magnitude of the torque is given with T arm x Fm If the torque is not zero the joint will rotate as inverse dynamics The aim of inverse dynamics is to calculate the forces and torques for known motions In mechanical engineering inverse dynamics is used to solve control problems The state of the rigid body is determined by and and they are our state variables Our state vector is X q q T and we have uh aw 2 5 Example free falling thigh With the above knowledge we start to model and simulate simple cases A free falling thigh for example The goal is to get familiar with some terms as center of mass inertia moment etc and use them In this case we assume the proximal end of the thigh to be connected with the hip and free at the distal end Our hip is fixed and only allows planar rotations Studying the sy
30. Modelling Martial Arts Techniques J C K Wong activation Muscles pulls Skeleton motion Articulated Skeleton V model SE PE Muscle Model A Thesis submitted for the degree of Master Scientific Computing Supervisor Dr Gerard L G Sleijpen Department of Mathematics University Utrecht October 2007 Contents 1 Introduction EI Martial arts sc 445 hai eh ee OUR BS SHA eae oa 1 2 Building a musculoskeletal model 00 1 3 Modelling dynamic systems 000 3L Simulation a 6 Bish es dik et ee Ae 1 3 2 Numerical solvers aoaaa 1 4 Modelling and simulation with Simulink 0 002 1 4 1 Modelling with Simulink 2 2 0 0000 1 4 2 Simulation with Simulink 0 1 5 Optimization problems 0 2 0004 2 Building skeletal models 241 IntroductiOtis lt lt haha of Bae ae Re ee He ee a 2 2 The skeletal model 02000000000 e 2 3 Simplificdtions ec ah osteo we oo ee ee we e 2 4 Dynamics of a rigid bodyin2D 0 2 4 1 Equations of motion 004 2 5 Example free falling thigh 000 2 5 1 Implementation with Matlab 2 5 2 Multi tbody dynamics 0 0 2 6 Modelling with SimMechanics 000 2 6 1 What is SimMechanics and why 2 7 Free falling thigh with SimMechanics 0 0 2 7 1 Building the model
31. a realistic hip is added Omin 5 deg maxz 110 deg e damping constant 100 ven e stiffness constant 40 a In figure 2 8 we see the block diagram Things become more interesting if we open the hip block in figure 2 9 In figure 2 10 we see the motion of the our free falling thigh With 49 90 degrees T 10 seconds This time the motion does not go from 90 to 90 degrees 2 9 Conclusion With SimMechanics we can very easily build a skeletal model and study its motion with known forces acting on it This is known as forward dynamics With forces we mean both rotational torques and translational forces We have built a skeletal model that models the thigh with a realistic hip Extending the model with more bodies is very easy The only forces that act on this model are gravity and force produced by ligaments In mechanical systems joint actuators are placed at the joints to drive the systems These actuators can produce any force at certain time if it is within their limit In the musculoskeletal system skeletal muscles are the joint actuators Their force production depends on many factors Controlling the musculoskeletal system is therefore more complex than controlling mechanical systems Luckily we have a very powerful controller 24 Elus hiro_skel_customjoint Hip Joint File Edit Yiew Simulation Format Tools Help Olea Sl st BOeleotrt 22 Joint Initial Condition joint_constants 110 5 100 40 hip jo
32. a smaller relTol needs more time But there are no big differences Comparing with the ODE32 table we see that sometimes ODE23 needs less time but sometimes more We are more interested in pulses of type rand There we see that ODE45 needs more time about 5 approximation errors and we should keep this in mind We directly use the different strategies EA because we think that these will do better 4 7 1 Experiment 1 f 40 We used a pulse array with frequency f 40 The values of the pulses are 0 or 1 Below are the basic settings for the EA e number of regions 4 e region size 40 solutions e migration time 5 generations e a unique EA strategy for each region The search is stopped if the the best solution did not improve after 3 migrations Results of the search are saved and the search restarted with random solutions as initial solutions and so on We let the search run for one night and in figures 4 5 and 4 6 we see the best result after 9 restarts The best found objective is 0 0224 which is quite small Wether this is a local or global optimum we do not know We can know it of course if we evaluate all possible solutions that equals 24 If we assume that it takes around 1 8 seconds for 1 evaluation for our computer then we should wait around 6 x 104 years If we use lots of computers however we can have the answer within our lifetime Before we try to do this or think of trying to do this we will do another experiment fi
33. a typical recorded EMG signal Currently our knowledge of the neural activation function is too low and therefore we will use a simplified function to model the activation function We have chosen our activation function to be made up with piecewise constant functions that we call pulses A pulse can only be zero or one All pulses have the same width and that depends on the number of pulses that we use The frequency f is a quantity for the number of pulses and is defined as the number of pulses per second The width w of a pulse is then w Now a solution is an array of binary valued variables Because we want to simulate the thigh for one second the size of the array equals f In figure 4 3 we see the case for f 4 Suppose we have an array of pulse of size f Then we can calculate act with f 1 Y pliait if O lt t lt 1 i 0 pf 1 if t 1 act t 4 3 where t 1 iw lt t lt 1 w er 0 elsewhere and p i is the i th variable of the pulse array 57 pulses 010 1 act 0 1 t s Figure 4 3 A binary valued pulse function with frequency f 4 Tmax 1 will get the value of the last pulse The extended raising thigh model The EA will find solutions in the form of pulse arrays These pulse arrays should be converted into activation functions In the extended raising thigh model we let the thigh model do the pulses to act conversion Suppose we have solution p which is an array of length f then the activation fu
34. actically impossible experiments This time we vary the mass and thus inertia The activation function is set to u t 1 In figure 3 13 we see the results for the normal thigh a thigh with twice the mass and a thigh with halve the mass Experiment 3 different activation functions In reality the activation function is very complex with varying amplitudes and frequencies In this experiment we used some standard functions Constant u t 1 Linear u t 0 8 t Sinusoid u t 0 9 0 1 sin zt 43 i sim_raising_thigh_00 DER File Edit View Simulation Format Tools Help Dee amp amp amp gt ar eS mp Normal o sim_raising_thigh_00 Models a raising thigh The hip is implemented with damping stiffness The rectus fermoris is responsible for the raising Figure 3 11 The raising thigh model Many blocks are masked to get a clear view The results are given in figure 3 14 Conclusion After all these experiments we are wondering if we can control the thigh Can we keep it still for example This can be seen as a control problem but also as an optimization problem We want to minimize the difference between the output angle and a reference angle In the next chapter will will use a self implemented genetic algorithm to find an activation function that can do this 3 9 Project Mae Geri It was very tempting to extend the thigh model Therefore we built a leg model that can simulate the mae geri or front ki
35. ation Below are listed 2 extensions to the basic EA 4 5 1 Multiple regions With this method we divide the population into several sub populations Each of them will evolve independently for some generations The number of gener ations where the sub populations evolve independently is called isolation time After each isolated evolution the subpopulation will exchange their best solu tion with each other This is called migration and there are several ways to do this After migration the sub populations will evolve independently and so on In algorithm 4 we see the algorithm for EA with different regions 54 Input Fitness function f Output Best solution found regions Initialize f random while terminal constraint is not met do foreach region do region evolve region n g end regions migrate regions end return Best solution population Algorithm 4 The different regions EA which is an extension of the basic EA 4 5 2 Different strategies The EA with multiple regions can be seen as an extension of the basic EA It can be extended further In 4 5 we have mentioned that we can improve the basic EA by choosing the right operators and their parameters A given set of choices for EA operators and their parameters is called an evolution strategy In most cases we do not know which are the better ones One solution to solve this problem is to give each subpopulation a different evolution strategy By doing thi
36. ave implemented an EA from scratch to solve optimization problems If we want to use the EA for our model we should formulate the goal in terms of the output and encode the activation function into an array of variables We used the EA to solve a problem that can be seen as a control problem The goal was to hold the thigh static for a moment of time For this goal we for mulated an objective function in terms of 0 and used the EA to minimize this objective function Encoding the activation function was a bigger problem In the ideal situation we can parameterize the activation put these parameters into an array and use EA to find the optimal array of parameters Because we do not how we can parameterize the activation function properly we used a simplified type of acti vation function These functions are made up with piecewise constant functions that we called pulses Real neural activation functions do not look like these functions With the simplified activation functions we get good results The EA we imple mented can also be used to solve other optimization problems After the experiments we realized that the approximated objectives can con tain errors if the used step sizes were too big With very big step size we not only have large interpolation errors but also risk that other activation is used instead of the found activation After a thorough examination with the best found solution we concluded that our approximated objectives are rel
37. ck It is a 3 segmented skeleton model with 6 mtc models The most difficult part was to couple the skeleton model and the mtc models Once this was done we can start simulation A major problem is that we did not know what activation we should feed the model besides zero and maximal activation In figure 3 15 3 16 and 3 17 we see some details of our leg model The most difficult part was to couple the skeleton and the mtcs During simulation with maximal activation for each mtc we encountered integration errors The precise reason is unknown but we suspect a bad initial state to be the cause of this In the initial state the system should be in rest Because there are 6 mtcs and 3 bodies finding a proper initial state is not easy In this thesis we did not solve this problem 44 Hip angle vs time various constant activation T T T T T theta deg time s Muscle moment vs time various constant activation netM Nm time s Figure 3 12 We used different constant activation functions to activate the rectus femoris Top the hip angle Bottom the generated moment The moment arm is a constant see table and Grieve constant Therefore the mo ment graph can be seen as a scaled force graph When the activation is 0 the thigh will not move and the angle is thus zero all the time Our mtc models consist of 1 fascicle with fast twitch fibers The activation threshold is set to ur 0 8 A
38. d and goes as follows Xati Xn AG at 1 9 Implicit methods that are used in practice have better stability properties than explicit ones per time step but they usually are more expensive 1 1 4 Modelling and simulation with Simulink In order to get numerical results we should implement the mathematical model on the computer and let the computer do the numerical integration There are several ways to do this Simulink is a very popular and user friendly platform to model and simulate multi domain dynamic systems For more information we refer to www mathworks com products simulink A Simulink model is a block diagram that represents the set of equations that describe the system It also has a wide variety of built in numerical integrators to do the simulation Throughout this thesis we will use Simulink as the modelling environment for our muscu loskeletal models and the ease of simulation with Simulink will become very clear 1 4 1 Modelling with Simulink In Simulink we use a block diagram to represent the equations of one model The block diagram is made up of blocks and signals A block can be an operator constant output etc Signals send time dependent values between the blocks The values can be of different data type double integer and multi dimensional Relevant blocks can be grouped into a larger block that represents a subsystem for example It is very common to build subsystems in a way that they have an one on one corresp
39. ded Catt ions This equation models the re lation between the stimulation and active state in a very simplified 4People with a weak stomach are not recommended to read this Many cats were sent to heaven 50 and I stand for origin point and insertion point of the mtc The origin point is closer to the body than the insertion point 33 Tendon PE ee Figure 3 4 Schematic view of the Hill s muscle model With Hill s model we want to predict the force behavior of the muscle tendon complex mtc with inputs neural activation and mtc length lor The muscular part is modelled with two elements a parallel elastic element PE and a contractile element CE The tendon part is modelled with a serial elastic element SE The force of the SE is the output of the muscle model way The active state is very difficult to model because it depends on many factors bio chemical bio electrical etc Hill based muscle models differ the most in how they model ay Virtual Muscle for example uses a set of coupled ODEs to model this relation 34 The algorithm for Hill s model In algorithm 1 we see how we can simulate Hill s muscle model Input Neural activation function act t and mtc length loz t Output Tendon force Fsg for t 0 dt T do Fsg fse lor 2 Fpp fpx 2 For Pea fcelz fo cE 4 af where f fap act ar af af a Fse Freet Fce 44 dtz z ax2 dtz end
40. e various activation function 80 T a const y Y lineair i oO n I f 3 j 70H j 1 theta deg time s b Figure 3 14 Different types of activation functions a The activation func tions b The joint angle 0 We can see some relations between activation and joint angle 47 ci sim_maegeri_VU File Edit Yiew Simulation Format Tools Help DSHS taele le muse_maegeri skel_maegeri c gt m jos Normal J Repe Maegeri model Maegeri is the Japanese term for a frontkick With this model we can model and simulate the maegeri The model is divided into a muscular and a skeletal part The muscular part contains 6 muscle models while the skeletal model contains 3 segments representing the thigh shigh and foot Simulation is done by feeding 6 activation functions into the model which will determine the activation of the 6 muscles The output is a leg motion Optimalisation is finding the proper activation function for a given goal An obvious goal eg is to let the ball of the foot make contact with a target in front of the leg the fastest way File Edit View Simulation Format Tools Help DISES reles A ga 2 o D gt fis Normal X Banga skel_maegeri This is the skeletal part of the maegeri model The human leg is modelled as 3 rigid bodies interconnected via joints The model is activated by netjoint moments around the joints caused
41. e muscle force The work done is given with 61 Ff They should be equal to each other thus a 0 F loiF Rearrangement leads to Or 1i B 3 With infinitesimal small steps we have a dlo 0 With B 2 we have an equation for loi Differentiation of this equation with respect to 0 gives us the equation for the moment arm a A Ao 0 B 4 If Ag is zero we have a constant moment arm Grieve constants are experimen tally determined The length of the mtc is measured at various joint angles and regression analysis is used to yield the Grieve constants see B 81 Appendix C The GA operators For our EA we have implemented several operators and methods We only list their names and their parameters For more details we refer to 1 The mutation and recombination operators depend on how we have encoded the solutions The implemented operators are suited for multi valued variables not for real valued variables C 1 Population Parameters Population size N The number of solutions we look at every time If N is too small we risk convergence to local optima If N is too big we need too much computer resources Evolution time G Once a generation is created the evolution time increase with one If G is small we risk we risk premature termination Streak s If G is very large but there is no improvement for a number of gen erations it is likely that our EA has found the optimal solution or got stuck in a local optim
42. e outcome of the force length force velocity and effective activation relation The FL block models the force length relationship The shortening and the lengthening block models the force velocity for concentric and eccentric contraction respectively The effective activation re lation is fairly complex and it takes diverse phenomena into account Some of them are fiber type dependent 77 D 5 g P 2 2 2 EI g B Z g 7 21145 14 1 1Tt ef lt 0 UCZ WO 09406 740 050493 u 3 LTE eff gt 0 UJU YIK FS Se 005 Figure A 6 The effective activation block It outputs a value between 0 and 1 that represent the effective activation This block takes all complex observed phenomena into account like yielding sag rise and fall times etc 78 Appendix B Grieve s method B 1 Calculating the mtc length with Grieve con stants In order to compute the muscle force we need to compute the length of the muscle tendon complex mtc I 1 loi depends on the joint angles of the joints they cover Mtcs that cover one joint are called mono articular mtc and mtcs that cover 2 joints are called bi articular mtc In figure B 1 we see a sketch of how we can compute the loi of a mono articular mtc We call this the radian approach and we have loil lo rx0 B 1 where r is the radius of the joint and l is the length of the mtc for 0 0 The way we computed loi however is with Grieve s method Here the joint lo
43. e quite limited In figure 4 10 we see a plot of the calculated 0 vs time where we zoomed in at t 1 There we see that solving with a relTol 1e is accurate enough Then there is the question is the activation function really used The width of our pulses equals 5 12 5e7 If we look at the maximal time step for relTol 107 in table 4 3 we see that it equals 8e 3 This means that every pulse is taken Furthermore we are talking about the maximal time step and therefore we conclude that the activation function is really used by 3This is only true if the used step sizes are very small compared with 0 001 the distance of the output points 65 Force vs time 250 T T 200 150 force N 50 0 f f f 1 L 1 f fi 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 time s a Force vs time 200 T T 180 F force N 2 3 R E a Oo o o o T T T T T T EN ts T L L i I L 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 time s 0 fi 1 b Figure 4 9 The tendon force vs time a The force plot belonging to the best found solution of experiment 1 Here we used pulses with f 40 We honestly did not expect this shape We expected a less or more constant force at the end b The force plot that belongs to best found solution of experiment 2 f 80 The produced force is more constant at the end 66 Found 6 vs time various relative tolerances 10 16 relTol 12
44. e triangular instead of constant pulses b The used activation vs time variable time steps Here triangular pulsus are used We can see that when large step sizes t 0 3 were used pulses are missed The simulation is done with a different activation 70 Found 6 vs time 0 12 T T 0 1 F 4 d Fi f a f 0 08 F i r f f as f w x 8 0 06 F F 4 Oo KeA Fd 5 Fd Fd 0 04 H f 4 af 7 F 0 02 we 4 0 il L eh ere sill f f 0 3 0 4 0 5 0 6 0 7 0 8 0 9 4 time t s a Figure 4 13 c Results of experiment 3 continued The found 0 vs time Omax 0 1 degree With triangular pulses we have less maximal activation and thus muscle force and we expected a smaller maz but not this small When the state of the model is less dynamic larger time steps are used Non zero pulses that cause the muscle to produce force risk being missed by the model 71 Conclusion We wondered if we could use the computer to find martial arts techniques Mar tial arts techniques are optimal motions for a certain goal We therefore have to model the locomotor or musculoskeletal system of the human body and perform optimizations on the model The musculoskeletal system is very complex and can be split into the skeletal and the muscular system Both are highly dynamic systems To get results on screen we should write software We have chosen Simulink a platform for modelling and simulation of multi domain dynamic s
45. ematical modelling A mathematical model is a simplification of a real world system It uses a set of equations to describe the system The complexity of the equations depends on the system and how real we want it For dynamic systems we have 3 main types of time dependent variables Input variables U t User given input The variables are known State variables X t The minimum number of variables that describe the state of the system at a certain time Output variables Y t These are the variables we are interested in They depend on the input and state variables We have Y t g U t X t The state of a dynamic system changes over time and the rate of change de pends on the input In physical systems the current state typically depends on previous states The type of equations that describe the relations are the well known ordinary differential equations dX gt 3 Z Xl O TO 1 1 1 3 1 Simulation If we want to study the behavior of a dynamic system over a certain time T we should simulate the model This is done by specifying the state at the beginning X 0 and solving equation 1 1 X 0 is also known as initial conditions vector X 0 Xo ra F 1 2 a X t dt X t 1 3 2 Numerical solvers Equation 1 2 can only be solved in some very simple cases In general we should use numerical methods to solve it The easiest way to understand how this can be done is with the Euler method With Taylor series we
46. ents lie parallel to each other with partial overlapping resulting in a striated pattern see figure 3 2 Skeletal muscles are also called striated muscles In other type of muscles the sarcomeres are not aligned in series These muscles are known as smooth muscles Examples of smooth muscles are heart and kidneys Furthermore a sarcomere is bordered with Z lines 3 2 Activation of the skeletal muscles Skeletal muscles are voluntary muscles We can activate them at will The small est unit that we can activate individually is the fascicle A fascicle contains only one type of fiber There are two main types of muscle fibers slow twitch and fast twitch They differ in the way they perform their metabolism Slow twitch fibers have relatively large endurance whilst fast twitch fibers are not Each fascicle is controlled by a a motoneuron that is located in the spinal cord A fascicle with its motoneuron is called a motor unit The motoneuron innervates all the muscle fibers of its fascicle with an action potential The larger the fascicle the larger the action potential 2 As a result each muscle fiber will fire a new action potential that will travel along their length see figure 3 3 This will result into the release of Ca ions With the release of Catt ions 1Muscle strain is something that we are all familiar with It indicates broken Z lines Muscle strain is usually the result of eccentric contraction In eccentric contractions the musc
47. eus Maximus VU Force N Force N MT path length cm Rectus Femoris Vu Activation MT path length cm Vastii VU Force N Activation MT path length cm MTClenath om atrocnemius VU From10 Soleus VU Goto11 From11 Figure 3 17 The participating mtcs for our leg model There are 6 of them 50 3 10 Conclusion We have seen how we can build models for skeletal muscles or muscle tendon complexes mtcs with Virtual Muscle The force producing behavior is very complex and currently there is no theory that describes the behavior Virtual Muscle models the mtcs empirically We hope that the behavior can be extrap olated and thus get realistic results for our musculoskeletal system A very hard problem is to get realistic data for the mtcs If we want to build musculoskeletal models the skeletal model and the mtcs should fit together To obtain these data lots of practical measuring experiments should be done For us this is almost impossible Luckily we have data from a jumping program supplied by the VU Amsterdam Virtual Muscle exports its mtc models as Simulink blocks and we can couple them relatively easily with skeletal models that we build with SimMechanics the Simulink environment The simulation is then very easy We have added a hip flexor to our thigh model from the previous chapter and done some experiments with it The outputs of our thigh model were very unexpected
48. good 63 Fitness evolution of all regions 0 25 T T T r region 1 region 2 region 3 region 4 0 2 F 4 Fitness 1 2 3 4 5 6 7 8 9 migration a Best objective per migration 0 16 T T T 0 14 4 O12 4 O1 4 o 2 8 0 08 4 T O 0 06 F 4 0 04 H 4 0 02 M M E 0 L L L L L L L 1 2 3 4 5 6 7 8 9 Migration b Figure 4 8 Best result of experiment 2 continued a The best objective evolutions of the 4 different regions during the best search The time is given in migrations A migration counts for 5 generations b Evolution of the best objective of the whole populations The best found objective 0 0180 64 relTol stepmin StePmax to Fa le de 70 0 014 0 75 0 0196 le 3e719 0 012 0 75 0 0180 le Te 18 0 008 0 75 0 0176 le 4e718 0 008 0 75 0 0175 Table 4 3 The results of simulating with the best found solution for different relative tolerances relTol stepmin and stepmaz is the minimal and the maximal step size respectively that the model used to is the start time that is used to calculate the objective Surprisingly t 0 75 for all relTol This indicates that very small time steps are used around t 0 75 F is the new objective approximation where we also take the variable time steps into account functions of both solutions we can not really tell the difference An interesting wa
49. he thigh to rotate in the XY plane This is done via a revolute joint with axis of action 0 0 1 The revolute joint models the hip 4 Add the thigh We will model the thigh with a rigid body SimMechanics s body block represents a rigid body In the body block we should fill in its mass and inertia tensor with respect to the center of mass Then the The inertia tensor is a 3x3 matrix that specifies resistance of the body against rotation In 2D there is only rotation around one axis in de Z direction Therefore we only need to fill in I33 19 cI ushiro_skel x File Edit View Simulation Format Tools Help DISHE sBeBleot 22 gt flo Noma l Be ushiro_skel This model models a free falling thigh The hip is a revolute joint and is fixed This system can be seen as a swinging pendulum Initial joint angles can be specified in IC block Joint angle q and joint reaction forces Fr are the modefs output The scope shows the joint angle over time directly Figure 2 6 The thigh model in Simulink built with SimMechanics It is a block diagram A SimMechanics block models a physical object The joint sensor block for example can sense different quantities like joint angle joint angular velocity joint reaction force etc 20 configuration of the body should be specified With 3 coordinate systems we can specify the main points of the body They are the center of mass the proximal end and the distal end The
50. http www fotosearch com 42 Name Rectus Femoris Mass kg 1 215 Fax N 4500 limar cm 44 lee opt cm 8 1 lse opt cm 35 36 Grieve const hip Ao A1 A2 0 110 0 035 0 Grieve const knee Ao A41 A2 0 308 0 042 0 nr units 1 Table 3 2 The data we used for our rectus femoris The grieve constants are needed to couple the rectus femoris with the thigh We only use one unit to model the whole muscle group the knee Because we do not have a knee we assume the knee to be at zero degrees for the rectus femoris At this position the virtual shank is thus aligned with the thigh The rectus femoris flexes the hip If it is activated and the generated moment is larger than the gravitational moment the thigh will rise In table 3 2 we see the data we used to model and couple the rectus femoris Therefore we named the model raising_thigh 3 8 Experiments Now that we have coupled the skeletal model with a mtc model we can start experimenting In all the experiments the initial angle is set to zero that is trunk and thigh aligned The solver is set to ODE45 with reltol 1076 Experiment 1 varying the activation In these experiments we use a constant activation act t c as input and simulate the model for 2 seconds In figure 3 12 we see the result for c 0 0 0 8 and 0 9 Experiment 2 varying mass and inertia One of the greatest advantages of musculoskeletal models is the ability to do pr
51. i relationship is measured via experiments Then a second order polynomial is fitted The found coefficients are called the Grieve constants J loi 5 Ajo Ain x 6 Aiz x 6 B 2 i 0 For mono articular muscles we have j 0 For bi articular muscles we have j 1 Grieve constants depend on the way 0 is defined Usually the null state 0 0 is when both segments are aligned Notice that when A2 0 we have a similar equation like with the radian approach B 2 Calculating the moment arm with Grieve constants Once the geometry of the joints and segments is known we can determine the moment arm This can be rather complex because each muscle has its own geometry A more generic way to solve this problem is with virtual work If we rotate the joint over an infinitesimal angle 6 an amount of virtual work is done lo stands for origin and is closer to the trunk i for insertion and is more away from the trunk 79 Figure B 1 An intuitive way to determine loi Seg i and seg i 1 represent the midline of 2 segments The circle represent the covered joint On top we have the zero state that is where 0 is defined as zero We assume that the MTC will follow the joint closely which leads to a radian increase This is however not always the case A clear example are the flexors of the knee There we have a straight line crossing 80 The work done to cause the rotation is then 6M 6a 0 F We also know that loi will change by th
52. iable In another small experiment we used triangular pulses instead of constant pulses Here we see that big step sizes will produce unreliable objective approximations We thus should not only focus on the relative tolerance but also pay attention to the maximal step size 4We used the default st pmaz which is set to 1 50 of the simulation time In our case it equals 1 50 0 02 the default stepmaz is thus used 69 Activation of best solution SS 3 S SS SSS Ss Se l 0 9 F 4 0 8 F 0 7 H 4 0 6 lt Bos H 0 4 H H 03 4 0 2 4 0 1F H 0 m y 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 time s a activation vs time 1 P I plat ron uR he Ki T F Eo 0 9 it Lt eH T HaT jel of 7 4 H l Il l 0 8 sak PEE aE al ht Ht J 0 7 jE ii ial r k l ik os Wil I Wad ik ji K K ji 05 mir a 0 4 0 5 0 6 0 7 0 8 0 9 3 time t s b Figure 4 12 Result of experiment 3 a The activation function belonging to the best solution where the pulses are constant This is a reference for the case that we us
53. int Models a realisti hip joint CounterMoment is added when extrema of the joint are exceeded Joint constant are specified in join_constant block Ready 100 lode45 Figure 2 9 The more realistic hip joint The joint sensor senses the joint angle In the extra moment block we calculate the contribution of the hip ligaments With the joint actuator we can actuate the hip with this contribution 25 Hip angle vs time realistic hip 100 T T T 6 degrees kkk k k 20 f f f f f 0 time s Figure 2 10 The motion of the free falling thigh with realistic hip started from a horizontal position The minimum angle is set to 5 From the plot we can see that it is not a hard limit When the minimum angle is exceeded counter moment is added This can be seen as an elastic collision in real life and that is the human mind In the next chapter we will build skeletal muscles that produces torques for the skeletal model 26 Chapter 3 Building skeletal muscle models The skeletal muscle is a complex biological structure Its main function is to generate force through contraction The force is transferred to the skeleton via tendon and aponeurosis that connect the muscle with the skeleton This force allows us to move but also protects our bones against impact forces Bones are relatively weak against bending forces compared with compressing forces By contracting the right muscles the bending force can be par
54. ion problem where we want to minimize the difference of the hip angle 6 with a reference angle Ores for a certain amount of time Can the EA solve it 4 6 1 Using EAs If we want to use EAs to solve this problems we should do the following 1 Formulate the objective function F this function tells us how good the simulation is It is defined in terms of the model s output 2 Encode the solution our solution is an activation function act This function should be encoded into an array of variables The evolutionary operators depend on what type of variables we use 3 Choose Evolutionary Strategy s if we use the basic EA we only have to choose one evolutionary strategy If we use the different strategies EA we should of course choose more evolutionary strategies 4 Perform evolutionary search we set some stopping criteria goal reached maximal search time etc and start the search Formulating the objective function The objective function tells us how good a solution is The definition of objective functions is problem dependent The EA only gives solutions for what we asked for and therefore a well defined fitness function is of paramount importance This time we want then the thigh to move up and stand stationary from t 0 75s to t 1 0s We choose the objective function F to be the r2 norm of A t 10 for 0 75 lt t lt 1 We have Fe Vl 10 2dt 4 1 We can only approximate the objective function We let the model
55. joint moments as input whilst the mtc models outputs muscle force and needs mtc lengths as input If we want to couple both models we have to calculate the mtc length and the moment arm of the mtc s force 8 The mtc lengths and moment arm both can be calculated with Grieve method see appendix B In figure 3 11 we see our previous thigh model coupled with the rectus femoris model The rectus femoris is a bi articular mtc see figure 3 10 It flexes the thigh and extends 8moment a F 9In mechanical systems torque drivers are used to activate the joints They only cover one joint 41 El controlled_ushiro_ga mtcs Rectus Femoris Vu_ss File Edit View Simulation Format Tools Help Die S BBE OP Qe r gt wh pn Sealan e Activation Out Activation Frequencies f0 5 Activation Frequencies f0 5 Recruitment Length L0 Velocity L0 s fascicles CE amp PE Convert to FO Force N series elastic element muscle mass Lee L0 path length for init Cem Vee ae Force F0 MT path length crn 100 Figure 3 9 A small error that comes with Virtual Muscle 3 1 5 In this pic ture the error is already fixed In the original blocks the force signals were interchanged These forces are used to calculate the length of the CE elements racer ine IUDYRK Crheom Figure 3 10 A lateral side view of rectus femoris The rectus femoris flexes the hip and extends the knee Source
56. kel The solver is set to ODE 45 with relative tolerance reltol 10 We have chosen 0 to be a joint coordinate The initial condition is 0 90 which means that the thigh is parallel to the ground The simulation time is set to 10 seconds Simulating the model with Matlab command line goes with sim ushiro_skel In figure 2 7 we see the output of the model This time with joint reaction forces Assigning the model s output to a variable y is done with t x y sim ushiro_skel Here t is the time vector x the state matrix and y the output matrix The dimension of y equals nt x n where nt is the number of time steps and n is the number of output variables So each column i of y contains the value of variable i over time When using variable time step solvers the time steps will vary Sometimes we want to have values at specific time steps We can specify an output time vector The output values are interpolated T t_start dt t_end t x y sim ushiro_skel T 21 Hip angle vs time 100 80 1 o L o Ll a o s u p e 204 40b 60H 80 100 10 time s Joint reaction forces vs time N eq04 10 time s se forces ensure Figure 2 7 a The hip angle It agrees with what we have already found but this time the angle definition is different b The joint reaction forces that
57. ld fit together This table however does not contain those needed data Therefore we used muscle data from a cycle program Soest The data for different mtc s of the leg are listed in table 3 1 As we see we are missing some parameters Most of the missing parameters can be derived however The only parameters that we do not have are the U and Lmaz We choose U 0 8 for all mtc s and we guessed Linaz mtc for each mtc such that Lmao 1 3 Building units Real skeletal muscles consist of many units They can easily sum up to more than 100 Each unit behaves differently within the muscles The way they are recruited depends on the size and fiber type of their fibers The complexity of the model scales with the number of units we want to model Remember that we need 4 state variables for a unit Is is very common to group similar units into big units BuildMuscles can automatically divide the whole muscles into units TIn cartoons we associate thick muscles with lots of power 40 BuildMuscles musc_maegeri_vu_new mat File Edit Muscles Model Help Fiber Type Database human_fiber_types mat Clipboard contents Empty Fiber type distribution PCSA _Units Muscle name Muscle FascicleLo Muscle MuscleFo Tendon Whole Fascicle S mass g cm PCSA N LoT cm muscle LMax Lo LMax cm Muscle 1 Biceps Femoris vU 1386 7 104 I Muscle 2 luteus Maximus VU 3333 3 x2 n i Muscle 3 Rectus Femoris Vu
58. le 73 Appendix A How does Virtual Muscle work The equations that Virtual Muscle uses to describe the behavior of a mtc is given in chapter 3 While some are relatively easy to understand others are not With BuildMuscles we can export mtc models as Simulink block diagrams From figure A 1 to figure A 6 we will explain how the output force of a mtc block is calculated We have built a mtc model that models the rectus femoris We only used one unit to model the whole muscle group In3 Activation Force N aa MT path length em f rectus1 Rectus Femoris Wu_ss In4 Figure A 1 The top level view of a mtc block The input and output are clearly given 74 Force N MT path length cm O Activation Frequencies 10 5 Activation Ng Frequencies f0 5 Recruitm ent Length L0 Force N r gt velocity LOS m fascicles CE amp PE series elastic elem ent perene muscle mass Lee L0 Path length for init om Figure A 2 View under the mtc block We have a recruitment block a fascicle block a serial elastic block and a mass block The recruitment block turns the normalized activation into stimulation for each fascicle The fascicle block cal culates the fascicle force The serial elastic block calculates the pulling force on the skeleton The mass block keeps track of the fascicle length and contraction speed Frequency O 5 Reached recruitmen
59. le force is overcome by a larger force and will lengthen instead of shorten If the difference is too big microscopic tears will take place in de Z lines 2 An action potential is a local voltage change that travels along the membrane of a cell 29 Charge mV 50 4o 0 as 50 bis Time 0 1 2 ms t t 4 4 4 4 4 4 Figure 3 3 The action potential of a muscle fiber top and how it is propagated through the fibers bottom Notice that the action potential starts from the middle of the muscle fiber where the junction of the fiber and motoneuron is formed Image taken from 9 cross bridges can be formed and the fibers will contract The force result from an action potential is called the twitch force Slow twitch fibers produces slower twitch force than fast twitch fibers The fire frequency is the number of action potentials per second The force production depends on the fire frequency The sum of electrical activation of the fascicle can be measured with electromyog raphy EMG The output is called a myogram and depends on the number of activated fibers fiber type motor unit composition and many more The am plitude of the signal corresponds with the activation level A fascicle has an activation threshold It will only be activated if the activation level exceeds this threshold The threshold depends on the fiber type and is proportionate with the size of the fascicle number of muscle fibers 3 The process of how fa
60. les Note in our original publication cited above the parameters listed for the slow twitch fibers had an error which was corrected in a subsequent corrigendum the correct values are used in this database Next 5 Fibers Figure 3 7 With BuildFiberTypes program we can build muscle fibers by spec ifying different coefficients The behavior of the muscle fibers are normalized in optimal force and optimal state variables A fascicle is treated as a scaled fiber the optimal fascicle length corresponds with its resting length So in a relaxed state the muscle has the greatest force production potential 3 6 2 Building the muscles With the BuildMuscles program we can build mtc models by specifying a num ber of parameters The parameters that scale the normalized force state rela tions to real values are Muscle mass m g the mass of the muscle Optimal fascicle length L cm average length of the muscle belly to pro duce optimal isometric tetanic force Physical cross sectional area PCSA cm can be seen as a measure for the amount of muscle fibers It is automatically calculated via PCSA Venyacte oe Pmuscle is the muscle density and is assumed to be 1 06 cm3 6The starting positions of sprinters and weightlifters make use of this knowledge 39 id Name Function Lo em lor cm F N 1 Biceps femoris hip extensor knee flexor 10 4 37 0 4000 2 Gluteus maximus hip e
61. m is a very complex function and a system of coupled ODEs should be solved to obtain for The more units we use the more time the simulation will cost In figure 3 5 we see a schematic view of the relationships between the different elements In figure 3 6 we see the way Virtual Muscle models a mtc 3 5 3 The model s variables Input variables the model has 2 input variables 36 Tendon amp Muscle Tendon amp Aponeurosis Aponeurosis SE1 SE2 100 E Muscle s5 Figure 3 6 Schematic overview of how muscle tendon complexes are modelled with Virtual Muscle Here the muscular part is divided into multiple units Each unit is modelled with a CE and PE element During sub maximal activation not all muscle units are active Smaller units have lower activation thresholds than bigger units 37 e Neural activation act it should have a values between 0 and 1 0 means no activation all the fibers are at rest 1 means full activa tion all the fibers fires action potentials at their maximum frequen cies With this value the effective firing frequency for each fascicle is calculated This models the recruitment The result is scaled to fo 5 fo s is the firing frequency where half of the isometric force is reached Each fascicle has its own fo 5 Recorded data from EMG could be used as input The data should be scaled to the level of maximal voluntary contraction Another source is data from a simu lated a mot
62. nction can be calculated with equation 4 3 By doing this we guarantee ourselves that the activation function will not be interpolated In figure 4 4 we see a screen shot of the extended model that we named sim act Benchmarks With the choice of EAs to solve the problem a large number of simulations should be done Usually we speak in terms of 10 thousands Therefore we did a little benchmark with our computer AMD 2600 WindowsXP for different activation solvers and relative tolerances The results are given in table 4 1 and 4 2 We are not happy with these results Each simulation takes around 1 8 sec onds Because there are no very big time differences we used the ODE45 solver with RelTol le for our simulations In this way we hope to get simulations with a relatively high accuracy 4 7 The experiments We let the model produce 0 at specified time points The distance dt between these time points is fixed and equals 1 1000 107 The used solver is the ODE45 variable time step solver with relative tolerance 1076 which is quite accurate The approximation of the objective is calculated with equation 4 2 Output 0 and thus interpolated values of 0 are used This will result in extra 58 E sim_act 7 lox File Edit View Simulation Format Tools Help D seSl st Beleot l2 gt wf Normal gt pulse2sct converter fi rasing thigh model sim_act This model is an extension of the raising thigh model The input of thi
63. nts and should be normalized to optimal state variables Virtual Muscle treats the muscle fiber as scaled sarcomere Therefore it uses normalized force state relation to describe the fibers behavior If we know the force state properties of a muscle fiber then the force state properties of the fascicle is the product of the number of muscle fibers A measure for the number of muscle fibers is the physical cross sectional area PCSA The PCSA is given with PCS A fasc Volumefase When we know length fase the optimal length of a fascicle and its PCSA we can calculate its optimal force with a constant called specific tension N cm The specific tension for human muscle fibers is around 31 8 ay and is experimentally determined For humans 38 BuildFiber Types human_fiber_types mat File Edit Fiber Properties Help Clipboard contents Empty Optimal Sarcomere Length um Type 1 Fiber type name ss Recruitment rank 1 V0 5 LO s 05 10 5 pps 6 Vf 20 fmin f0 5 05 T 05 fmax 10 5 gt lf a so Comments super slow al typical A generic fast twitch A e g soleus slow twitch fiber fiber type type v x v v Additional Comments All three fiber types were presented in Cheng et al J Neurosci Meth 2000 They were created by adapting the analogous feline fibertypes as described in that paper to fit data from triceps surae hamstrings forearm and hand musc
64. o solve var ious optimization problems EA has some major drawbacks however First we do not have certainty that the search will converge to an optimal solution This because it is a stochastic search method Second EA is very time consuming We easily speak in terms of thousands of fitness evaluations In our case one fitness evaluation means one simulation Luckily they are very well suited for trivial parallel implementations where we can use multiple computer processors simultaneously Other methods and operators can be used to extend the basic algorithm and improve its performance Input Fitness function f Output Best solution found population Initialize f random while terminal constraint is not met do parents select population f children recombine parents children mutate children population reinsert parents children end return Best solution population Algorithm 3 The basic Evolutionary Algorithm EA 4 5 Extending the basic EA The basic evolutionary algorithm is able to solve many complex problems Their performance can be improved by good choices for the GA operators and their parameters When the fitness functions are cheap we can afford creating large populations and allowing large evolution time Sometimes we have very expen sive fitness functions and we have to use more sophisticated EA A simulation is then very costly and with better EA we hope to get better results with less comput
65. objective Under the microscope we will see that a muscle fiber is made up of parallel myofibrils In figure 3 1 we can see the structure of a skeleton muscle A remarkable observation is the striated pattern that can be seen A myofibril is made up with small contractile units in series called sarcom 27 Periosteum covering the bone Tendon Fascia Skebtal Fasciculus Endomysium Muscle _ _ ares fiber cell 4 Endomysium Striations Sarcolemma Sarcoplasm P Filaments E a Myofibrile Figure 3 1 The structure of a skeletal muscle The muscular part is grouped into muscle bundles or fascicles Each fascicle contains a number of muscle cells or fibers A muscle fiber is made up of parallel myofibrils The myofibril in turn is made up of sarcomeres in series Image taken from 13 28 Nuclei i sti inert te il SERES a heb ett ih pee i beeee Muscle fiber HEH yi ell me o amenna Figure 3 2 Skeletal muscle fibers under a light microscope We can notice a striated pattern The white rings are formed by actine filaments while the myosin filaments form the darker rings The nucleus is an organelle that contains the fiber s genetic information Image taken from 13 eres Each sarcomere is built up with thick and thin filaments called myosin and actin filaments respectively Myosin filament has a dark color while actine filament has a lighter color Both filam
66. on with relTol At time n we have absTol relTol x max a t 1 5 Optimization problems We build a musculoskeletal model with different types of questions in mind Most of these questions can be translated into optimization problems These questions contain keywords like best most less worst etc For example the best kick the most economical kick Optimization is a different field in mathe matics One type of optimization is to find best or worst output for any kind of input We rate the output with an objective function F Y This function F tells us how good the results are We name the function that maps the input to F h h U F see figure 1 3 Our problem is to find an optimal input Uses such that h Oopt lt h U or Ass gt A Y if we want to minimize or maximize respectively It is very likely that h has many local optima We will implement an evolutionary algorithm that should be able to solve this type of optimization problems lWe can experience this in martial arts practice too We always find a better way to execute a technique U dX dt f X U y F Y X bime Figure 1 3 Schematic view of a dynamic system The output Y depends on the state X and the input U and is rated with objective function F To get the state of the system we should use numerical integrators Function h is very complex and with each simulation we get to know one more point of h Our optimization problem is optimizing h
67. ondence with real systems By doing this we can get clear and easy to understand models 2 1 4 2 Simulation with Simulink To simulate the Simulink model we should specify an integrator and its accu racy Simulink has several ready to use built in numerical integrators The rtol x Region in which rtol determines acceptable error ak ee a P Region in which atol determines acceptable error State Figure 1 2 The specifications of regions where the approximated state variable will be accepted Source 2 2 main types of solvers are fixed and variable time step In dynamic sys tems variable time step methods are recommended For real time simulations Simulink restricts us to use only fixed step solvers When no real time results are needed the recommended choice is the ODE45 variable time step solver with relTol 1078 The ODE45 solver is based on the Runge Kutta method of order 4 and is an explicit method It is a variable time step solver and its accuracy is set by specifying the relative tolerance relTol and absolute absTol The relTol specifies in which region x is accepted and absTol specifies the threshold when x is near zero State x is only accepted if all errors e satisfy ei lt max a x relTol absTol 1 10 otherwise the solver will use a smaller time step In figure 1 2 we see the different regions If we set absTol to automatic Simulink will start with an absTol 1076 and readjust it during simulati
68. oneuron act t is also the input for our musculoskeletal model In our optimization process we typically want to find act t e Length loi the length of mtc and depends on the skeletal configu ration Because the skeletal model outputs joint coordinates q the loi q should be calculated The state variables The state of the whole muscle model is determined by the length of the fascicles lcg and a set of variables that determines the active state of the fascicles The model uses 4 variables to specify the state of a fascicle The output variables The output is the force of the tendon Fsp in Newtons and has a positive sign If we want to couple the muscle models with the skeletal model we should take this into account 3 6 Using Virtual Muscle 3 1 5 Virtual Muscle 3 1 5 comes with two programs BuildFiberTypes With BuildFiberTypes we can build fibers by specifying their specific coefficients The relations are normalized to maximal force Fmaz and optimal state variables Because fibers can be scaled to fascicle level this program can be seen as BuildFascicles BuildMuscles with BuildMuscles we can build mtc models and export them as Simulink blocks Data for the mtc we want to model are needed 3 6 1 Building the fibers The muscle fiber data we used came with Virtual Muscle 3 0 It contains co efficients for super slow SS slow S and fast F human muscle fibers The coefficients are obtained after lots of real muscle experime
69. rent choices Furthermore they are parameterized In ap pendix C we see an overview of the operators that we have implemented These operators are implemented from 14 4 3 Terminology We assume that our audience are familiar with EA Below are some terminology we used Population A set of potential solutions for our problem With EA we want to evolve the population to a better one and so one until termination Mostly we are only interested in the best solution when the search is terminated Sometimes we are also interested in other good solutions local optima The size of the population is denoted with N Region Several publications have shown that it is better to split the entire population into a number of sub populations or regions More in Solution A solution is an array of variables that represents the real solution In biological terms we can compare a solution with a chromosome 1 Variable A solution is an array of variables When we used the term chromosome for a solution then the term gene is the equivalent for variable The three most common types of variables are binary valued multi valued and real valued Objective function a function that measures how good a solution is IEA and GA are inspired by processes in nature Some EA specialists nowadays want to break the link with biological terms 53 4 4 The basic Evolutionary Algorithm The basic EA algorithm is listed algorithm 3 and is quite powerful t
70. riangular shape see figure 4 11 We use our best solution from experiment 2 and calculate the activation function with 4 5 act t p t a t 4 5 i ll where i 2 t iw iw lt t lt i zY t 1 ut 4 9 ot iwy 4 sw StS i w 0 elsewhere Then we used this activation function to simulate our thigh model The solver is the well known ODE45 with relTol le The details of the simulation are e stepmin be 18 st Pmax 2e very large 68 Imax 0 1 degree e F 3 8369 interpolation e F gt 4 9767 real output values There is a big difference between F and F2 This could be the result of the relative large stepmax 4 We can prevent this by specifying stepmaxz manually In figures 4 12 and 4 13 we have some interesting plots of the used activation and the found 0 Here we see that specifying stepmaz is very important The step size used by the solver depends on the state of the system The state of the system in turn depends on the activation When relatively big time steps are used parts of the activation will be missed In this case the largest time steps st pmax 0 02 were used The width of a pulse is w 0 0125 and some pulses will be missed The simulation is done with a different activation than what we have in mind This is a wise lesson for the future We should specify stepmaz small enough otherwise another activation is used for the simulation 4 8 Conclusion We h
71. rst 4 7 2 Experiment 2 f 80 This experiment is the same as the previous one except we changed the fre quency to f 80 With f 80 the search space is much larger than with f 40 around 10 times and maybe we can find a better solution We also let our computer do the search for one night and the result after 4 restarts around 8 5 hours is given in figures 4 7 and 4 8 The number of restarts is lower than the case where f 40 This can be explained with the increased search space if we take f 80 The found objective is smaller than with f 40 0 0180 vs 0 0224 and the plot of the found vs time looks better If we look at the found activation 60 Activation of best solution 1 e U6 eS M M a eee eee 0 9 F J 0 8 F J 0 7 4 0 6 4 0 4 4 0 17 4 m a 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 time s a Hip angle vs time 12 I T T T q deg D fi L 1 L L 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 4 time s b Figure 4 5 Best result of experiment 1 a The found solution activation function We can not see a clear pattern in it b The found 6 vs time It is close to what we want and it look quite good if we also take numerical errors in mind 61 Fitness evolution of all regions 0 35 T T T T T region 1 region 2 region 3 0 37 region 4
72. s model is an array of pulses From this pulse array the activation will be calculated in the pulse2act converter Ready 100 jode4s 7 Figure 4 4 The extended thigh model that we named sim_act The input are the found pulses solutions and the pulse2act block will calculate the activation function at needed time points With ODE23 time s RelTol pulses 1 pulses 0 pulses rand 1078 2 8224 1 7862 1 6384 1074 1 7267 1 8196 1 6757 1075 1 7357 1 7251 1 7108 10 6 1 7580 1 7340 1 7459 Table 4 1 The time it takes for our computer to evaluate a fitness function The frequency is set to f 40 The pulse array is thus an array of length 40 pulses 1 means all variables are 1 pulses 0 means all variable are 0 pulses rand means the variables are randomly chosen to be 0 or 1 Once the array of random variables are generated it is used over and over again for the benchmarks In this table the ODE23 solver is used The time is the average of 10 simulations There is not much time difference between different relTol When we used our EA we should count on around 1 8 seconds for an objective function evaluation This is a serious bottleneck 59 With ODE45 time s relTol pulses 1 pulses 0 pulses rand 10 3 1 8420 1 3410 1 7662 1074 1 8205 1 4211 1 7656 1075 1 8680 1 3953 1 8378 1076 1 8874 1 3657 1 9026 Table 4 2 Same as 4 1 but with ode45 In most cases
73. s look more realistic In the next section we will use a multi body software package called SimMechanics to implement our skeletal models 17 Hip angle vs time 0 T T 0 5 1 g 1 5 2 og 2 H 25 tf 3 3 5 0 1 2 3 4 5 6 T 8 g 10 time s Hip angular velocity vs time 8 T T T 6 f i H 3 33 z fii F Hi t gt 0 k t Pte 33 2 h 2 t AHE 5 6 7 8 9 10 time s Figure 2 5 Top hip angle as function of time is the angle between the positive x axis and the thigh Bottom hip angular velocity 0 as function of time Initial state 0 0 0 0 0 0 Because we have a free rotational case and no loss of energy 0 will go from 0 rad to 7 rad back and forth 18 2 6 Modelling with SimMechanics In the previous section we have implemented the free falling thigh with Matlab Although it only took a few lines of code things become more complex if we want to add more features and or bodies Another way to study the dynamics of multi body system is the use of dedicated software packages Such a package is SimMechanics In 7 SimMechanics is chosen to be the best multi body systems software for integration with a muscle software package called MMS More infor mation about SimMechanics can found at http www mathworks com products simmechanics 2 6 1 What is SimMechanics and why The way of modelling with SimMechanics differs from the way we are
74. s or wild animals Sometimes with weapons and sometimes barehanded Through centuries people have discovered that some movements give better re sults than other Martial arts techniques can be seen as special movements meant to defend or attack Martial arts are the arts of adapting these tech niques to ones self and using it creatively Martials arts nowadays are mostly practiced for better health which is also a way to protect ourselves 1 2 Building a musculoskeletal model We are interested in the motion of a human body The human body is a very complex biological structure and consists of many systems The system we are interested in is the locomotor system The locomotor system consists of a skele tal system and a muscular system that drives the skeletal system Both are very dynamic systems The locomotor system is therefore known as musculoskeletal system To bring down complexity we will only model part of the total muscu loskeletal system the leg with its skeletal muscles for example In chapter 2 we will see how we can build a skeletal model in Simulink with SimMechanics In chapter 3 we will use Virtual Muscle to build skeletal muscles and add them to the skeletal model The number of skeletal muscles depends on the complexity of the skeletal model After this we can start simulation and eventually extend the model The input is an activation signal for each skeletal muscle These signals stand for the neural activation of the muscles
75. s we bet on more horses and thus hopefully have more chances for winning Furthermore we believe that with different strategies the subpopulation will help each other The most obvious example is the case where we use different muta tion strategies In the beginning the strategy with large mutation will do better while after some time and migrations the strategy with small mutations will do better The algorithm for EA with different strategies is given in algorithm 5 Input Fitness function f various evolution strategies EA strategy Output Best solution found regions Initialize f random while terminal constraint is not met do foreach regions do region evolve region EA strategy end regions migrate regions end return Best solution population Algorithm 5 The different strategies EA which is an extension of the differen regions EA 55 4 6 Project Static Thigh In the end of the previous chapter we wondered if we could keep the thigh still The reason for this is to find a so called operational point In this point the system should be in rest We can save the state of this point and start other experiments from this point In project Static Thigh we want to raise the thigh from an initial state and keep it still at 10 degrees for a short period of time The initial state was chosen to be the state where the thigh is at rest 6 0 and perpendicular to the ground 0 0 This can be seen as an optimizat
76. sci cles get activated is called recruitment Slow twitch fibers have lower activation thresholds and fascicles with slow twitch fibers are usually smaller than those with fast twitch fibers In slow movements like typing walking the fast twitch fibers are not activated They are only activated when explosive movements are needed This can be explained with the economic metabolism of the slow twitch fibers 3 2 1 The cross bridge theory The basic theory behind muscle contraction is the cross bridge theory Huxley 1957 For a basic explanation we zoom in at sarcomere the level There 3This is known as the Henneman size principle 30 sarcomere is the basic contraction unit and has an optimal length of around 2 7um It is made up with actine and myosin filaments named after the main protein molecules they are build with and titin filaments Actine and myosin filaments are also named as thing and thick filament respectively according to their diameter A cross bridge is a temporary bound between thick and thin filaments Once bound it will pull the thin filaments towards the center of the sarcomere and let loose The titin filaments keep the thick filaments in the center of the molecule The number of cross bridges that can be formed depends on several factors The length of the sarcomere the cross bridges can only be formed if there is overlap between the thick and thin filaments The cross bridge theory assumes that cross bridges are
77. specifications can be done in many ways Once the position of the origins of the frames are set we should also specify their orientations The output of our revolute joint is the difference of the orientation of the ground frame and the proximal frame In our case we can let the orientation of the CS es to be default In the 3D case the orientation of the CS of the center of mass depends on how the inertia tensor is given 5 Add a joint sensor The joint sensor can sense joint motions and joint reaction forces 2 7 2 Simulation Now that our model is built we can start the simulation First we need to specify the initial condition for our joint The zero state of the joint is determined in the specification for the thigh s CS s During simulation numerical integration should be done SimMechanics can be fully used with Simulink and we can use its implemented numerical solvers directly Output can be generated via Simulink s output blocks With SimMechanics we can literally view our out put via animation Animating the system s motion can significantly slow down the simulation however Then we should fill in the simulation time Simula tion is started by pressing the play button All these are done via the GUI of SimMechanics For serious experiments it is easier to run the model program matically With Matlab we can write scripts that simulates the model under different circumstances Free falling thigh We named our model ushiro_s
78. stem The hip only allows rotation and thus our system has one degree of free dom DoF We choose an inertial frame such that its origin coincides with the proximal end of the thigh See figure 2 1 The body s configuration is given by 0 The other two DoF are constrained by the following equations fa pcos 2 11 Yom psin 0 a There are two forces acting on the thigh One is the gravitational force F and the other is the internal or joint reaction force F F ensures the pure rotational 13 Figure 2 3 The forces that act on the thigh Fy is known and acts on the CM along the y axis a is the moment arm of F with respect to the hip F is the internal reaction force from the trunk motion around the hip and is exerted by the fixed trunk These forces are drawn in figure 2 3 While F mg is known F should be calculated Because the hip is fixed and we do not know F we write Euler s equation with respect to the proximal end of the thigh XO M Inrit By doing this we do not need to know F explicitly because its moment contri bution with respect to the hip is zero We have X T p x F mgpcos and our equation of motion is thus mgp cos Inipd If we are interested in F Hay Pay we can use Newton s equations and equation 2 11 to compute it Newton s equations are 5 Fr Mzicom 2 12 XF miem With X F F and gt F F F we have FF mz EN 2 13 Fr mycm Mg om and
79. t act 0 8 the fascicle then becomes activated and starts to fire action potentials at its lowest frequency With act 1 0 the stimulation is max imal and the fascicle will fire action potentials at its highest frequency With act 0 8 and above we expect the thigh to raise That it will fall afterwards and raise again and so on we honestly did not expected We might explain this with the force length and force velocity relations of the fascicle During the raising period the length will decrease and the contraction velocity will increase All this will result in less force and gravity will take over When falling however the fascicle length will increase and we also have eccentric contraction 45 Hip angle vs time various mass and inertia 120 T T T T T T o normal halve twice 100 F M 80 F 4 theta deg D o T time s Muscle moment vs time various mass and inertia 90 T UF i T T T T saai normal halve twice H 80 707 netM Nm 407 30 F 1 2 14 16 18 2 time s Figure 3 13 We changed the mass and thus the inertia of the thigh The lighter the thigh the higher it will raise This agrees with our expectations 46 Activation vs time various activation function sinus E06 0 4 F 0 2 4 0 L 1 L L L L L L 1 0 0 2 0 4 0 6 0 8 1 1 2 1 4 1 6 1 8 2 time s a Hip angle vs tim
80. th binary valued variables It has the following parameters Number of mutants the number of children where mutation operator will be applied to Mutation rate p each variable of a solution will mutate with chance p The recommended value is 1 1 where l is the number of variables of a solution 84 Bibliography 11 12 13 L F Shampine Numerical Solution of Ordinary Differential Equations Chapman amp Hall March 1 1994 Simulink 6 Using Simulink The MathWorks Inc 2006 http www mathworks com Grant R Fowles George L Cassiday Analytical Mechanics 6th edition Saunders College Publishing 1990 Forward Dynamics of Multibody Systems A Recursice Hamiltonian ap proach Vrije Universiteit Brussel Brussel September 2005 M A Me Conaill J V Basmajian Muscles and Movements a basis for human kinesiology Baltimore Maryland 21202 USA 1969 SimMechanics For Use with Simulink User s Guide Version 2 The MathWorks Inc Natick MA 2006 http www mathworks com P Montazemi R Davoodi Comparison of Dynamic Engines for Muscu loskeletal Modeling Software MSMS A M Institute for Biomedical Engi neering University of Southern California 2005 David A Winter Biomechanics and motor control of human movement Hoboken New Jersey 2005 B M Nigg W Herzog Biomechanics of the Musculo skeletal System 2nd Edition John Wiley amp Sons Ltd West Sussex PO19 iUD England 1999 Ernest
81. the thigh will feel The formulation is given in equation 2 13 The the rotational motion Due to gravity that acts in the y direction F is larger 16 944kg the gravity force equals 166 22N The y component reaches peaks that are much larger than the gravity force Joint reaction forces are interesting for those who want to study injuries than Fp With mass Minigh 22 id Name Length m p m Mass m kg Inertia Icm kg m m 1 Thigh 0 485 0 210 16 944 0 418 2 Shin 0 458 0 198 7 074 0 136 3 Foot 0 165 0 045 2 468 0 020 Table 2 1 Some realistic segment parameters id Name O0min deg Omax deg Stiffness Damping 1 Hip 5 110 100 40 2 Knee 150 0 1000 40 3 Ankle 5 45 1000 40 Table 2 2 Joint parameters The extreme angles are given with respect to the joint coordinates in a right handed coordinate frame The zero joint state is defined as the state where both bodies are aligned We can also change the solver and its settings The default solver is ODE45 with relative tolerance relTol 1073 Changing the solver with other relative tolerance goes with myopt simset Solver ode23 RelTol 1e 3 Zother solver reltol t x y sim ushiro_skel T myopt 2 8 Skeletal data To build a skeletal model we should use realistic data for the body segments and joints These data are obtained via special measuring
82. the right choice Virtual Muscle treats the muscle as a scaled sarcomere and we can easily build realistic muscle models if we have good macroscopic muscle data The if is a very big if because the skeletal and muscle models should fit together Finally it was time to perform optimizations We had implemented an evolu tionary algorithm EA with a wide variety of operators from scratch We used our EA to solve a problem that can be considered as a control problem In that problem we have modelled a thigh that is driven by a skeletal muscle and wondered if the thigh could be kept still at a certain degree Wether this is possible we did not know The EA produced results that are close to what we want A major drawback of using EA is the large number of simulations that is needed A simulation is very costly even for our simple thigh model A way 72 to handle this problem is to implement the model more efficiently on another platform The Simulink model can then be used to verify the program A more rigorous way is to use more computers The EA is well suited for trivial parallel programming Although we did not find a martial arts technique we have made a first step at how we can do that During this thesis our respect for the human body and hu man mind has grown more and more The human body as a piece of wonderful engineering and the human mind as a great wonder One should keep both of them in good shape by practicing martial arts for examp
83. tly or mostly neutralized We are mainly interested in the force production behavior of the muscle In this chapter we first look at the structure of the skeletal muscle and how it contracts Then we look at how muscle are modelled mathematically After that we build a thigh muscle model and add this to our thigh skeletal model 3 1 Structure of skeletal muscles Skeletal muscles are connected to the skeletal bones via tendons on both ends The proximal attachment point and the distal attachment point are named origin and insertion respectively Different connective tissue sheaths having their origins on both tendons run from end to end through the muscle subdi viding the muscle in units subunits till the muscle cell or fibers These sheaths are thus continuous and prevent the muscle from pulling out during contractions The largest unit is the skeletal muscle itself It is surrounded by a fascia and a connective tissue sheath called epimysium epi above and my muscle Then the muscle is subdivided into muscle bundles or fascicles by a connective tissue sheath called perimysium peri around Within a fascicle we have the muscle fibers or muscle cells kept away from each other by endomysium endo internal Each muscle fiber is surrounded by a delicate membrane called the sarcolemma Unlike other cells the muscle cell is shaped long If we want to zoom in we have to use a light microscope that can reach a magnification of 1000 times the
84. tthresh old Fetivation Figure A 3 The recruitment block It models the recruitment of the muscle units In this block the activation is turned into stimulation stim for each fascicle or muscle unit The numbering of input is done from top to bottom The switch is set to recruitment threshold u If input 1 is larger than the recruitment threshold it will pass else the value will be set to zero so it will not be recruited 75 Force N max Force 0 eh Convert FO to N Sum Forces FO Sum all Units FO Fpe1 0 01 u 1 23 0 046 log exp u 23 2302 1 17 0 046 1 Type SS Motor Units VL Dem uxaill Units Frequencies 10 5 Velocity LO8 0 02 exp 21 0375 u 0 70222 1 Length LO Figure A 4 The most complex block S the mtc block The PE1 block comes into action when the fascicle exceeds a certain length The PE2 block models the thick filament compression It comes into action when the fascicle decreases to a certain length The CE block is fairly complex and we will look at it later on The calculated force is first normalized to the optimal force Fo F Fo ER AP FL FXAPE2 Sum all Units vbtor Unit 1 PCSA1 Demuxall Units 1 0 35 1 exp abs u 0 19 Fvilengthen Figure A 5 The CE block revealed It has shortening length shortening speed Fpe2 and fenv fire frequency as input and the normalized fascicle force as output The CE force is a product of th
85. um Further searching can be a waste of computer resources By specifying s the search is stopped if there is improvement after s generations C 1 1 Region Parameters Region size n the population is divided into regions n is the number of solu tions per region Isolation time g the regions will evolve independently from each other for a number of generations This number is specified with g the isolation time C 2 Initializing the population The initialization is done with randomly generated solutions 82 C 2 1 Terminal constraints In our case the search space is enormous Furthermore the costs for evaluating our fitness function are very expensive We can stop the search based on the following criteria Goal When the goal is reached the search can be stopped This is an ideal situation We should remark that we can not always define a goal Time s When the searching time exceeds a user specified amount of time the search will be stopped Costs fitness evals A more objective stop criterion is the amount of fitness evaluations Streak There is a chance that the EA is stuck in a local optimum Streak is the number of generations that the search does not improve C 2 2 Selection For each generation we should select a number of members to produce offsprings Truncation nmates nmates is number of parents for the next generations The minimum is 2 Tournament nmates sp sp is selection pressure needed for linear
86. xtensor 20 0 15 0 5000 3 Rectos femoris hip flexor knee extensor 8 1 34 4500 4 Vastii knee extensor 9 3 16 13500 5 Gastrocnemius knee flexor ankle extensor 5 5 37 6 4000 6 Soleus ankle extensor 5 5 23 56 8000 Table 3 1 Data for different leg muscles taken from a cycling program from de VU Amsterdam Optimal force F N the maximal isometric force It is calculated as the product of the PCSA and specific tension The specific tension for hu mans is about 31 85 Optimal tendon length Lo r the tendon length when optimal muscle force acts on it Tendon has a very adaptive behavior When its resting length Lo r is given we can approximate it with the optimal length 1 04Lo 7 Maximal muscle tendon length Liar mec em the maximal length of the whole mtc Activation threshold U fractional activation level threshold If the input activation U exceeds this level all fascicles or units will be recruited and fire action potentials At U 1 all fascicles will fire at the maximum fire frequency Maximal fascicle length Lmaz Lo the maximal fascicle length is calculated with fone To Gathering the muscle data is one of the most difficult steps we encountered The Yamaguchi table 12 is the most well used table If we want to incorporate the muscle models with the skeleton model we should also know their length prop erties as function of the skeletons configuration In brief the skeleton model and the mtc s shou
87. y to compare both solutions is by looking at the modelled tendon force that belongs to these solutions In figure 4 9 we see the plots 4 7 3 How reliable is our best solution With the 2 experiments we have found 2 solutions Based on the objective of these solutions we declare the found solution of experiment 2 f 80 as winner Because we used numerical results there will always be numerical er rors rounding off numerical integration interpolation etc involved Can we trust our results The simulations in the EA searches used the ODE45 variable time step solver with relative tolerance relTol 1076 The reason for choos ing a variable time step solver is efficiency The solver will only use smaller time steps when needed For the calculation of the objective function however we used interpolated 0 This will lead to numerical errors and the question is how bad it is Therefore we examined the best found solution thoroughly We feed this solution to our model and simulate the motion a few times Each time the ODE45 solver is used but with different relative tolerances The calculated 0 instead of the interpolated 0 is then studied The objective function is now approximated with N 1 i 0 where to 0 75 and ty 1 Relevant quantities of these simulations are given in table 4 3 In this table we see that Fz 0 0176 relTol 1e and does not differ much from F 0 0180 which is good news It means that the interpolation errors ar
88. ystems to implement the musculoskeletal system Throughout the thesis it became very clear that using Simulink is about the easiest way to model a musculoskeletal system Especially with the help of SimMechanics and Virtual Muscle 3 1 5 Because of the complexity we modelled only parts of the musculoskeletal sys tem We started with the skeletal system first The skeletal system is divided into several segments and rigid bodies were used to represent these segments This model is also known as an articulated model and is a multi body system The complexity of our musculoskeletal system depends mainly on the size of the skeletal system The more bodies it contains the more muscle models should be added We used SimMechanics to model the skeletal system With SimMe chanics it was very easy to model and extend the skeletal system With the choice of modelling the skeletal system as a system of rigid bodies we moved into the field of multi body dynamics Here the dynamics are well under stood In the world of muscle modelling things are different Scientists still do no know how the muscle behaves precisely We have chosen Virtual Muscle to model the behavior of muscles Virtual Muscle is an empirical skeletal muscle model It can predict the muscle s behavior at sub maximal activation relatively well compared with other muscle models Because the maximal voluntary ac tivation of human muscles is still sub maximal we believe that Virtual Muscle is

Modelling Martial Arts Techniques

Contents

Download Pdf Manuals

Related Search

Related Contents