An Optimization-Based Model for Full-body Reaching Movements

- 1 - An Optimization -Based Model f or Full-bod y Reaching Movement s Daohang Sha 1 *, James S Thomas 2 * § 1 Department of Orthopaedic Surgery, Yale Medical School, New Haven, CT, USA 2 School of Physical Therapy , Ohio Universit y, Athens, OH, USA *These authors contributed equally to this work § Corresponding author Email addresses: DS: sha@ohiou.edu JST: thomasj5@ohiou.edu - 2 - Abstract Background The development of a simulation model of full body reaching tasks that can predict end- effector trajectories and joint excursions consistent with experimental data is a non-trivial task. Because of the kinematic redundancy inherent in these multi-joint tasks there are an infinite number of postures that could be adopted to complete them. By developing models to simulate full-body reaching movements in 3D space we can be gin to explore cost func tions that may be used by the central nervous system to plan and execute these movements. Methods A robust simulation model was developed using 1) graphic- based modeling t ools to generate an inverse dynamics controller (SimMechanics), 2) controller parameterization methods, and 3) cost function criteria. An adaptive weight coefficient based on the final motor task error (i.e. distance between end-effector and target at the end of movement) was proposed to balance motor task error and physiological cost terms (e.g. joint power). The output of the simulation models using different cost controller functions based on motor task error or motor task error and various physiological cost terms (e.g. joint power, ce nter of mass displacement) were compared to experimental data from 15 healthy participants performing full body reaching movements. Results In sum, the best fit to the experimental data was obtained by minimizing motor task error , joint power, and center of mass displacement. Simulation and experimental results demonstrated that the proposed method is effective for the simulation of large-sca le human skeletal systems. - 3 - Conclusions This method can reasonably predict the whole body r eaching movements including final postures, joint power and movement of COM using simple algebraic calculations of inverse dynamics and forward kinematics. - 4 - Background Full body reaching tasks require the central nervous system (CNS) to apportion motion to the legs, trunk and arms in order to successfully complete the task. The two primary constraints are that the end-effector makes contact with the target, and that the body’s center of mass remains within the base of support (i.e. at target contact). However, even with these constraints, there are an infinite number of joint configurations that can be used to complete the task due to the inherent kinematic redundancy in the human body. Given an infinite solution set, how does the CNS plan and execute coordinated movements in a kinematics redundant system. From a mathematical perspective, optimal control theory is an effective method used to solve redundant systems. In fact some have suggested that the CNS uses a method similar to optimal control to coordinate multi-joint movements [1, 2]. In optimal control, the apportionment of motion to the various joints is determined by an iterative process that attempts to minimize certain input criteria. Var ious criteria for the optimal control of human motor coordination have been proposed based on empirical findings, physiological phenomena or both. Flash and Hogan (1985) reported that smoothness of the end-effector trajectory in Cartesian space is optimized (i.e minimum end- effector jerk) [3], while Rosenbaum et al . (1985) found that smoothness of joint space is optimized (i.e. minimum joint angular jerk) [4]. Still others have proposed that joint torque [5], joint power [6]; or displacement of whole body center of mass (COM) [7] is minimized. Finally, Nakano et al , 1999 proposed that final end-effector error is minimiz ed as well as norm of torque change (i.e. the norm of 1 st order derivative of joint torque) [8]. In motor control simulations, the goal is to determine the apportionment of joint motions based on some optimized criteria (e.g. minimal joint torque) and this requires that the optimal controller input joint torques or muscle forces to produce joint motions [5, 9, 10]. This process is called forward or direct dynamics. This method requires an integration of the - 5 - differential equations of the multi-body systems, which can result in the entire algorithm being unstable, particularly when the gravitational and elastic forces are included in the dynamics model. While several methods have been proposed to improve the stability of these forward models [2, 9], they are still problematic in terms of the computation time required to perform these calculations, i.e. up to 10,000 hours of CPU time in desktop computers [5, 11]. This is because initial conditions must be checked for consistency with constraints after each of iteration of integration. An alternative approach uses an inverse dynamics model in which the optimal controller computes joint angular positions derived from a certain criterion. These joint angles are then used as input in an inverse dynamics model that computes the joint torques or internal forces [6, 12]. The inverse dynamics based method is much faster than the forward dynamics based method and eliminates the stability problem for the system dynamics due to its intrinsic algebraic structure [6]. However, a potential drawback of this method is that the inverse dynamics model only works on open topologies, i.e. open link chains[13]. While there are several ways to implement a dynamics system model, once the system’s degrees-of-freedom (DoF) increases, the required analytical expressions become unwieldy [2]. While these expressions could be determined using symbolic tools, a more efficient method is to use graphic-based tools for simulation of rigid body machines such as SimMechanics. SimMechanics requires the geometry of bodies and mass properties, possible motions, kinematics constraints, and the coordinate systems to initiate the model. It doesn’t require the user to develop the equations of motion independently [14]. Additionally , graphic tools allow real-time visualization during simulation, which is useful to validate the proposed model because the simulated movement should at least be qualitatively similar to what is observed experimentally [11]. - 6 - The purpose of this study is to develop an effective method/model to address how motor system performs full body reaching movements. In this paper, we present an optimal controller based on the parameterization optimization method using several different input criteria that satisfy the constraints of a full body reaching task. A method for using graphic based tools based on inverse dynamics to model full body reaching tasks is presented. Finally, the models are validated qualitatively and quantitatively to experimental data collected from healthy subjects performing full body reaching tasks. Methods Simulation Model In voluntary target reaching activities, the primary goal of the motor task is to ensure that the end-effector makes contact with the target. As illustrated in Figure 1, the first problem for the control sy stem is to complete the task with minimal end-effector error (task error). Cost functions generate parameters that are input to the polynomial controller, which outputs a set of joint angle trajectories. These joint angles are input to the inverse dynamics and forward kinematics models, which in turn provide output of joint torques, power, COM location and end-effector error. These are then used iteratively to compare with task constraint (i.e. end-effector reaching the target location) to refine the cost function. The cost functions are then minimized iteratively using an optimization algorithm (Figure 1). Each of the blocks of our model, as illustrated in Figure 1, is described in greater detail below. Cost function In general, the criteria or cost functions for voluntary target reaching movement can be written as (1) dt Ru u e e C f t T f T f ∫ + = 0 λ - 7 - where f e is the final motor task error / end-effector error vector with three elements representing errors in x , y , and z direction respectively (i.e. distance between end-effector and target at the end of movement); λ is a weighting coefficient that expresses the relative importance between the motor task error (first term) and the physiological cost (second term); R is a positive-definite matrix with proper dimensions, which indicates the importance of each joint involved in the full body reaching; u could be one of the following quantities or their combinations such as end-effector jerk 3 3 dt x d , joint torque τ , joint torque change dt d τ , joint power θ τ & × = P , or body center of mass ∑ ∑ = i i i C m x m x . All of these quantities are the function of parameters of the controller (see next section). It must be emphasized that all calculations of physiological cost functions include the effect of accelerations due to gravity . As presented above, the cost function has an end-effector error term and a physiological cost term. In this paper physiological cost refers to movement perfor mance measures such as total joint power, displacement of center of mass, or a combination of these two varia bles. The first term and second term are usually in different units (e.g. m 2 and (N.m) 2 ) and depending on the physiological cost term, their values can be vastl y different. Because this is a multiple objective optimization which requires the simultaneous optimization of more than one cost function, some trade-off between the criteria is needed to ensure a satisfactory movement prediction. Here we propose an adaptive weight coefficient that adjusts its value based on the final motor task error, i.e. f T f e e 0 λ λ = with constant 0 λ . So the criterion becomes (2) Now the constant 0 λ can be chosen so the value of the phy siolo gical cost term is equal to one when the motor task error satisfies the preset tolerance. When the two terms in brac kets       + = ∫ dt Ru u e e C f t T f T f 0 0 1 λ - 8 - are approximately equal then optimization convergence is assured. The constant can be calculated by 1 0 0 max −       = ∫ dt Ru u f t T λ . The maximum value of inte gration in the eq uation c an be ob tained by the primary simula tion with end- effector error co st as the only c riterion. Once the pr imary simulatio n has b een r un, and th en 0 λ can be dete rmined for eac h physio logical cost te rm (e .g. joint power, C OM). This en sures the end- effector reaches the targ et locatio n at the end o f move ment if the optimiza tion co nverg es. This is a necessary thoug h not suffic ient c ondition for conv erge nce. A ctually, once the cost f unction is deter mined, the conv ergenc e depen ds mainly on the behav ior of t he optimizatio n meth od. Polynomial Ty pe Controller The t ime histo ry o f the optimal trajec tory o f each join t is a fu nction o f time which can be appr oximated by a n n th -orde r poly nomial (3) The a ngula r velocity and ac celera tion can be derive d ana l y tically f rom equ ation (3) . These joint tra jector ies are then use d as inp ut for the inver se dyn amics and forwa rd kinem atics calc ulations (Fig. 1). The ad vantag e of usin g an inve rse dy namic s metho d is that the se calc ulations do not r equire integra tion, a nd since t here is a one- to-one mapping from jo int spac e to Ca rtesian sp ace in fo rwar d kinema tics, the proble m of k inema tic redunda ncy is elimina ted. T hus, by using a forwa rd kine matics calcula tion, the proble m is redu ced to deter mining the co efficients o f the p olyn omial. F urtherm ore, th e whole body voluntary rea ching m ovement can be partia ll y describ ed as the point w here the end -effector starts ( 0 = t ) mov ing (i.e . standing neutra l postur e wher e 0 θ is know n) to th e postu re ad opted at conta ct targe t within a certa in time ( f t ). Thus, the location of targe t in Ca rtesian spac e is know n but the final posture f θ is not. . ... ) ( 3 3 2 2 1 0 n n t p t p t p t p p t + + + + + = θ - 9 - Base d on equation (2) , with the initial cond itions , 0 ) 0 ( , 0 ) 0 ( , ) 0 ( 0 = = = θ θ θ θ & & & and f inal co nditions , 0 ) ( , 0 ) ( , ) ( = = = f f f f t t t θ θ θ θ & & & and us ing a 6 th orde r polynom ial as e xample, on e gets th e follow ing: and [ ] [ ] [ ] . 3 ) ( 6 1 , 3 ) ( 15 1 , ) ( 10 1 6 6 0 5 5 6 6 0 4 4 6 6 0 3 3 f f f f f f f f f t p t p t p t p t p t p − − = + − − = − − = θ θ θ θ θ θ (4) This meth od requ ires at lea st a 6 th order polynomia l to ensur e enoug h freedom for th e appro ximation of joint m ovem ents as is e vident f rom above . If the fina l posture is un known, the va riable f θ can b e tuned toge ther with othe r paramete rs such a s 6 p . Ther efore, there are only tw o variable s (i.e. fina l posture f θ and poly nomia l param eter 6 p ) that ne ed to be tun ed for e ach movem ent of each jo int. The h igher the order of polynomia l, the more f reedom the polyn omial has to approxima te the join t motion, h oweve r, as the po l yn omial ord er increa ses the numb er of v ariables that ne eds to be tuned also increa ses. Skeleton Dy namics The ge neral eq uations of fu ll body m otion can be written as (5) wher e θ θ θ & & & , , a re th e vectors of joint a ngle, angular velocity, angular acce leration respe ctively; ) ( θ I is segme nt mass inertia moment m atrix; τ is vecto r of net joint moment; ) ( θ G , ) , ( θ θ & V , ) ( t T are g ravity terms, Coriolis- centrip etal-visc oelastic ity, an d external terms such as gro und reac tion forces. Joint viscoe lasticity can b e written , 0 , 0 , 2 1 0 0 = = = p p p θ ) ( ) ( ) , ( ) ( t T g G V I + + + = θ θ θ θ θ τ & & & - 10 - as θ θ θ θ & & B K V + = ) , ( , where K and B are joint stiffness and viscoelastic coefficie nt matrices respective l y. One way to imple ment the dynamics system model is to use an alytical methods. Fo r systems with low D oF (i.e. DoF<3), the anal y tical expressions for the rela tionships between angular acc elerations and joint torques can b e easily written. Howe ver, a system with as few as two segmen ts with 7 DoF, requires a n anal y tical expression with more than 200, 000 elementary operations (e.g. +,-,*, cos, sin) [2]. Althoug h these expressions could b e determined using symbolic tools, it is still unimagina ble that how many e lementar y operations wou ld be needed for a full body model with 12 segments a nd 36 DoFs. Therefore, a more effic ient method is to use graphic-base d tools for simulation of rigid body machines such as SimMecha nics (The MathWorks, Inc .) or S I MM (Musculaographics, Inc.) . In this paper, a linked segment model for the inverse d yn amics of whole body motions in 3D space was developed using Ma tlab/Simulink and SimMechanics T oolbox. This linked seg ment model includes twelve segments (i.e. head, thorax, abdomen, pelvis, right/left ha nds, forearms, upper arms, thigh, a nd shank) and twelve joints (i.e. r ight/left wrist, elbow, shoulde r, cervical, thoracic , lumbar, hip, knee and a nkle). While this model was deve loped using only one leg because of inherent problems with solving inv erse dynamic s of closed loop systems using SimMechanics, the movements of the low er extremities are sma ll and nearly symme trical. Thus this single leg model should provide reasonab le results. To further simplify the model, each joint has only thr ee rotational degrees of freedom (D oF), (i.e. flexion/extension, internal/e xternal rotation and a bduction/adduction) within the sag ittal, frontal and transve rse planes. The total number of DoF for the model is 36. While others have included transla tion of the shoulder girdles in the ir models [15], we chose to use a more simplified model that could b e readily compare d to our experimental data. The inpu ts for the model (i.e. the joint a ngular trajectories a nd their derivatives) are - 11 - provided by the polynomial controller. T he outputs of the model (compu ted by SimMechanics), are the net muscle torques a nd forces for each joint, as well a s joint power and motion of COM, e nd-effector location ) ( θ g x = (forward kinematics) . These outputs then are used to reduce end-effector error and ph ysiolog ical cost through the use of a n optimization algorithm. Optimization A lgorithm The cost func tion is minimized subject to the equa lity constraint (n onlinear dynamics equation of motion, E q.5; initial and final bounda ry conditions, Eq. 4) and the inequality constraints (the limitation of joints) to obtain the optimal pa rameters of controllers, i.e . subject to where p is a parameter vector of controllers; therefore the cost f unction is also a function of the unknown c oefficients n i p i ,..., 1 , = and the fina l time f t ; 2 1 , ε ε ≤ ≤ ∆ C p are the stop conditions of the optima l algorithm. Once the initia l values for the paramete rs of controllers have been input, the optimization algorithm will modify the paramete rs until the preset criter ia or minimum parameter c hanges are satisfied, i.e. j j j p p p ∆ + = + 1 , where j is the index of iteration. I n principle, an y well deve loped nonlinear optimization algorithm can be used to find the optimal contro ller parameters. H owever, different algorithms w ill produce differen t results. Here we use the nonlinea r Least-squares function lsqnonlin in Matlab Optimal ToolB ox to perform the optima l parameter search . The Levenberg-Ma rquardt method with line searc h was used (More 1977), i.e. [ ] ) ( 1 j j T j j p C I J J p ∇ + − = ∆ − σ α , (7) ) ( min arg 2 1 , p C p C p opt ε ε ≤ ≤ ∆ = n i up i i low i ,..., 1 , = ≤ ≤ θ θ θ - 12 - where         ∂ ∂ = i t p C J with f t t ,..., 0 = and n i ,..., 1 = , is Jacob ian matrix; j α is the iter ation step which is d etermine d by line search; 0 > j σ is a positive constant; j j j p p C p C ∂ ∂ = ∇ ) ( ) ( is the grad ient of cr iterion C resp ective to contr oller par ameter p . Jacobia n matr ix ( J ) and the grad ient of cr iterion ( ) ( j p C ∇ ) are a pproxima tely calculated through parame ter perturb ation. Initializing Simu lation Model As note d ea rlier the c ost functio n is a c omposite function with a ta sk err or cost term and a phy siologica l cost te rm. Simulatio ns wer e con ducted usin g the follow ing criter ia cost func tions: 1) m inimize fi nal e nd-effec tor error only without any physiolog ical c ost ter m (i.e. min Er ror), 2 ) minimize final e nd-effe ct error and total joi nt powe r (i.e. min Powe r), 3) minimi ze final end-e ffector e rror an d body CO M displac emen t (i.e. min COM) , 4) minim ize final e nd-eff ector err or, total join t power and b ody COM displace ment (i. e. min Power +COM ). We a lso calcu lated th e hand tr ajectory using min imal jer k criteria (min Jer k) as de scribe d by Flas h and Hoga n (1985 ) [3]. The minimum- jerk tr ajecto ry of end -effector was c alculate d by the following equatio n,               +       −       − + = 5 4 3 0 0 6 15 10 ) ( ) ( T t T t T t x x x t x f , (8) when moving from loc ation 0 ) 0 ( x x = to f x T x = ) ( in T t = seco nds [3]. To initia lize the o ptimiza tion, th e initial jo int an gles at ( t= 0) an d ( t = t f ) were se t to the mean value s from ne utral stand ing po sture der ived fro m experim ental d ata, and the co ntroller para meters w ere set to zero. In this study , a 6 th orde r poly nomial co ntroller w ith two unknow n var iables, i.e. final p osture ( f θ ) and the coeff icient o f polyno mial ( 6 p ) is used to deter mine ea ch joint e xcursion. There fore, th ere are 72 param eters to be tuned for the 12 - 13 - joints. To ensur e phy siological fidelity of the sim ulation, joint ran ge of motion va lues (E q.3) wer e input b ased on a ccepted n orms of joint ra nge of mo tion [16], and me asure s of join t viscoe lasticity determ ined fro m the e xtant litera ture [17- 20] (See Table 1). The a verage height a nd weig ht of 15 he althy su bjects we re used in the simula tions an d mass- inertia l chara cteristics of ea ch seg ment were derived from the regre ssion equa tions pro vided in the lite ratur e [21, 22 ]. Fina l motor ta sk error (i.e. dista nce b etween target and e nd-eff ector), p arameter tolera nce ( 1 ε ) and criter ion tolera nce ( 2 ε ) of te rmina tion of o ptimization a lgor ithm we re set as 2mm , 10 -6 , and 10 -6 resp ective ly for a ll simulatio ns. The consta nt weig ht 0 λ was set as 1 0 - 8 , 10 -9 , 10 -10 (for high, middle a nd low targe ts) and 7×10 3 , 10 3 , 10 2 (for high, mid dle an d low targe ts) for power an d COM ter m respec tively. Fo r simplic ity, joint we ight matrix R w as set as the unit mat rix. The m ovem ent dur ation w as set as 0.56s, 0. 575s an d 0.68s for high , middle and low targe ts respec tively based o n the experime ntal re sults fr om 15 he althy subjec ts. The solver for numeric al com putation was set w ith a f ixed-step sample interva l of 0.001s. Experimental Protoco ls Fifte en healthy subje cts (7 male s and 8 female s with age 22.93 ± 1.7 9 year, we ight 6 8.59 ± 10.69 k g an d height 16 9.12 ± 7.74 cm) pe rforme d a serie s of reac hing tasks to three targets locate d in the mid-sag ittal plane at a fa st pace d speed (i.e . approxima tely 600 ms fro m initial posture to targ et con tact). Targ et locatio ns wer e standa rdized to the pa rticipant’s anthr opometr ics. Th e participan ts could , in theory , rea ch the ta rgets by flexing their tru nk 15, 30, or 60 deg rees, with their s houlder flexed to 90-degr ees and elbow exte nded w ithout any motion f rom the other joints [23 ]. Bef ore begin ning th e study , each p articipant w as infor med - 14 - of the experim ental pr otocol an d signed the co nsent form a pprove d by th e Ethics Co mmittee of the Ohio U niversi t y . Starting fro m an uprig ht standin g po sture, with each fo ot on a force p late, the partic ipant r eached with their right ha nd for the targe t. Subjec ts paused at the targ et for 1 seco nd and then retur ned to an uprig ht postur e. Three trials a t each targ et heig ht were perf ormed and the targ ets wer e presen ted from highest to lowes t. Motions o f the tr unk an d limb seg men ts were recorded at 120H z for 5 s econds using the M otionMo nitor Sy stem (Inno vative Sports Tr aining, I nc. Chicag o). T his sy stem can trac k the thr ee-dimens ional coor dinates of six-deg ree-of-fre edom magnetic sensors w ith a spa tial re solution of 1 .8mm in position a nd 0.5 deg in orientatio n (Asc ension TM , Floc k of Bir ds ® , Ascen sion Tech nology Corpor ation, B urling ton, VT, U SA). Th e mag netic sensor s were a ttached by Velcr o® stra ps to the lim b seg ments ( at the mid point betw een the joints) o f the r ight a nd left shank , thigh , arm, and fo rearm, as well as the thora cic verte bra (T1), lumbar vertebr a (L1), an d the sac rum. An Eu ler ang le seque nce was u sed to d erive the three dimens ional joint motions f rom the upper and lo wer extre mities bila terally , as we ll as the th oracic and lumb ar spine. These data wer e smooth ed with a 61-poin t fourth order Savitzky -Golay f ilter [24] a nd ser ved as input for inver se dyna mics calc ulations. T he sam e inve rse dy namics (Sim ulink mod el) mod el used in simulatio n was used to c alcula te the motion of COM, joint torqu es, an d joint pow er fo r each trial of each s ubject. The prop erties o f each segment, suc h as ma ss-ine rtial char acteristic s, size, an d loca tion of CO M etc ., were derived f rom anth ropome tric regr ession eq uations b ased on the mass a nd heig ht of eac h subject [21, 22]. T he expe rimenta l traje ctory of en d-effect or includin g dista nce fr om targe t location , velocity , accelera tion an d jerk also were ca lculated from the sam e Simlink m odel. - 15 - Data A nalysis The p aths/tra jectories of the en d-effect or determ ined f rom the simulation model were compa red q ualitatively to the m eans a nd stan dard erro rs of path /trajector ies of en d-effecto r from the expe rimenta l data. To quantita tively assess the s imulation m odel, for e ach targe t heig ht, t-tes ts were u sed to co mpare the pred icted COM dis place ment a nd predic ted tota l joint pow er to the exper imenta l data. Result s Optimization The o ptimization f or ea ch trial o nly ta kes a few hours of CPU time on a per sonal compu ter (see Table 2), howev er, there was no clear eff ect of c ost function s or targ et loc ation for th e CPU tim e req uired to pe rform the simu lations. Examina tion of to tal pow er s quared a nd final C OM squ ared (i.e tota l cost) in Tab le 2 shows that the values of these measure ments a re significa ntly red uced wh en adding the phy siologica l cost te rms (e.g. min Powe r, min C OM an d min Pow er + COM) to th e Err or cost te rm eve n thoug h we use d an a daptive w eight coe fficient to balan ce the phy siolog ical cost te rms. Com paring all four contr ol strateg ies, final CO M squa red w as smalle st when the min COM strate gy wa s used. H owever , total powe r squa red wa s smalle st for min Pow er + COM c ontrol s trategy . End-effector T rajectories Fig ure 2 illustr ates th e trajec tories o f the end -effector d etermin ed by 1) the simulatio n model using four d ifferent c ost function s, 2) c alcula ted minimum jerk, a nd 3) e xperime ntal da ta aver aged ov er 15 subje cts. Visua l inspec tion of this fig ure revea ls that f or eac h target h eight, the tra jector ies of the end-ef fector and their d erivatives a re qu ite similar to the e xperime ntal data regardle ss of th e cost crite ria use d in the simulation. T he be ll shape d velocity trace s are consiste nt with previou s finding s for two and three joint arm mo vemen ts. While th ese - 16 - traje ctory p lots sugge st that e ach sim ulation m ethod pro vides sim ilar re sults, exam ination o f the pa th of the end-e ffector g ives gre ater insig ht into the differ ences in these me thods. The p ath of the end-ef fector dete rmine d by the differe nt simula tion meth ods and the experim ental d ata (a s describe d above) is shown for ea ch plane and targ et heigh t in Figu re 3. Examina tion of th e exper imenta l data f or the path of the end-ef fector in the sag ittal an d fronta l plane s (Figur e 3A & B) reveal s a fairly large path cur vature, w hich is consiste nt with the pa th pre dicted by the simulation using eac h cost functio n. In contrast, using minim al jer k as a cost fun ction pre dicts a stra ight line path of the end- effector in all pla nes. How ever, examina tion of th e exper imenta l data f or the path of e nd-effec tor in the transve rse pla ne indica tes a m uch stra ighter pa th. Thu s, for the transver se plan e, the be st fit appea rs to be provid ed by e ither m inimal je rk or th e minima l task err or cost f unction. In contr ast, minim al powe r and COM cost f unctions pr edict a fairly la rge cur vature a nd do n ot provide a good fit of the experim ental da ta. None theless , as sta ted earlier , even if the pa th of the e nd-effec tor can be accura tely pre dicted, the movem ent stra tegy or the postur e ado pted at ta rget c ontact may not be consistent with expe rimenta l data. Thus, we next comp are th e predic ted CO M displac emen ts to the e xperime ntal da ta for eac h target h eight. Displacement of C OM The tr ajector ies of CO M as pre dicted by the simu lation mo dels usin g the different cost func tions are compare d to the trajec tories of C OM from the expe rimenta l data (F igure 4) . While, th e exac t path of the CO M does not app ear to be w ell fit by any o f the cost function s, the c hange in COM from initial p osture to targ et contac t appear s to be consistent w ith the cost f unction th at minimi zed join t powe r and C OM disp lacem ent. In fact, t-tests revea led that usin g a cost fu nction that m inimize d COM d isplac emen t and jo int power was the only cost f unction th at wa s not signif icantly differe nt from the e xperime ntal da ta (Table 3). Minimum task e rror p redicted the larg est disp lacement of the CO M. Spe cifica lly, for the - 17 - middle and low targe t reaching tasks , the pred icted COM displac emen t exceeds 2 0 cm a nd 30 cm re spective ly, whic h would c learly cause th e subject to fall f orward or require a step to prev ent falling . Ther efore minimizin g en d-effe ctor err or alone does not pr ovide a r easona ble solution f or the movem ent task with re spect to COM d isplace ment. N otice tha t when only task e rror is minimized, the postu re ad opted at ta rget con tact do es not c ompensa te for the forw ard displa ceme nt of the tr unk an d its effe ct on who le body COM (Figure 5) . Qualita tively , the p osture adopte d at ta rget con tact by participa nts (i.e. experimen tal data ) is best f it by th e cost fun ction that m inimizes COM displac ement and join t power (Figure 5). Total Joint Po wer Fig ure 6 sho ws the time serie s of joint power ∑ = 36 1 ) ( ) ( i i i t t θ τ & derive d from experime ntal data and from model sim ulations u sing th e variou s optimal control strateg ies. From this figu re, in ge neral, it a ppears th at the co st functi ons of m inimal (p ower + COM) a nd minima l power give the best q ualitativ e fit of the experim ental d ata fo r the h igh target only. Ho wever, for the mid dle and low ta rgets, only minimu m (pow er + CO M) app ear to pr ovide a g ood fit o f the da ta. To compare the cost f unction s to the e xperime ntal da ta qua ntitatively , we to ok the integ ral of join t powe r to ge t total energy (See Table 4) . For the high target and middle ta rget it appe ars tha t minimizin g CO M displa cement provides th e best fit to the experimenta l data . Howe ver, fo r the low targets, none of th e cost functions p rovide a good fit to total e nergy expend iture. Discussion Skeleton Dy namics The u nderly ing pre mise of the se simula tions mo dels is tha t an inv erse d y namics a pproach was th e prefe rred me thod for predic ting move ment stra tegies. I nverse dynamics uses as inputs the joint mo tions a s a fu nction of time a nd differ entiates th em twic e to yie ld the accele ration s require d to calc ulate jo int torque s and inte ractio n forces needed to p roduc e the - 18 - motions. Howe ver, th e motion fu nctions must be che cked for consisten cy to ensu re they stay within g eom etrical c onstraints o f linked rigid bodies at e very time fram e. For inv erse dyn amics with closed topolog ies, this proce ss is comple x and the compu tation load is q uite large . The difficulty in handlin g clo sed topo logies com es fro m indeterm inacy which a ge neric prope rty of the invers e dy namics itse lf [13]. I n the pre sent study , the f ull-body reach ing tasks requ ired the subjects to stand f irmly on the fo rce plates w ithout a ny foot moveme nt durin g the ta sk. The refore , there exists a closed topology within lower extremities. To av oid the consiste ncy p roblem s of ge ometrica l constrain ts mention ed ab ove, a one leg m odel was used in the in verse dynam ics calcu lations u sed in th e optimizat ion and simulation with mass prope rties do ubled f or thigh and shan k segme nts. For the expe rimenta l data, the same techn iques w ere ap plied. Give n the ta sk con straints, a one leg approximati on is a cce ptable for f ull-body reaching moveme nts and consiste nt with our previo us work [25]. Comparison of Simulation and Experiments The r esults of these simulations a re in g eneral agreemen t with the experim ental d ata rega rding th e displace ment of C OM an d the f inal posture s adop ted at d ifferent ta rget locatio ns. Per fect fits to the expe rimenta l data are not show n in this study becau se 1) fittin g is ge nerally a ssociate d with arb itrary parameter adjustmen ts[2], i.e. adjustin g so me anthr opometr ic para meters may help to improv e the p redic tions [26]; 2) fitting quality may not be suffic ient to e stimate th e valid ity of a model [27 ]. Fig ure 5b illu strates that minimi zing Pow er a nd COM in th e simula tion mo del pr ovided a reas onably clo se fit to the expe rimenta l data. Howeve r, it is a lso clear that there are diffe rences between the simula tion and the expe rimenta l data regarding location of the e nd- effec tor. Thi s is partic ularly evident f or reache s to the mid dle an d low ta rgets. The se - 19 - diffe rences ar e most likely due to erro rs associa ted with m otion ca pture ( e.g. skin m oveme nt under sensor s) and f rom model constra ints (e .g. constra ining shoulde r translation) . Shoulde r transla tion wh ich ca n dramatica lly aff ect loca tion of the end-ef fector, p articularly in forw ard and o verhea d reach ing tasks [28, 29]. Th e ver tical and a nterio r-posterio r transla tion cou ld be up to ± 3.8c m. How ever , even w ith these poten tial source s of err or, both sim ulation a nd experim ental r esults c apture the major characte ristics of the volu ntary ta rget reac hing movem ent, i.e . curve d path for end-eff ector, b ell-shape d prof ile of velo city of en d-effecto r and in creasing displac ement in C OM an d joint pow er w ith lowe ring the target heig ht. Controller The m inimum je rk soluti on fo r end-e ffector movemen t based on Eule r-Poisson’ s theore m indica tes that the optim al traje ctory of end-ef fector must h ave th e form of a 5 th -orde r poly nomial [3]. A lthoug h a 6 th -orde r poly nomia l was use d in our study, the value of the 6 th - coef ficient is v ery sm all for min imum er ror tra jectory (e.g. 5 2 6 10 1.9275 − × = p , 5 2 6 10 9119 . 3 − × = p and 5 2 6 10 4586 . 5 − × = p for high, middle and low targ et respectively). Thus, the minimum err or movement is similar to the move ment of minimum jerk of e nd- effector in ter ms of symmetric bell-shape veloc ity profile, sine-shape acceleration pro file, and parabolic jerk pr ofile. While for the minimum power a nd COM trajectory, the value of the 6 th -coeff icient is quite big relative to the minimum error trajectory (e.g. 0578 . 0 2 6 = p , 1341 . 0 2 6 = p an d 4489 . 0 2 6 = p for high, middle and low targ et respectively). The e ffects of the 6 th -coeff icients of controllers on the minimum pow er and COM tra jectories of end- effector are signif icant especially for middle and low target reaching mo vements. - 20 - Optimization Me thod In this study, the L evenberg-Marquardt method with line sea rch was used to find the optimal trajectory o f the end-effector in terms of motor task err or (end-effector error) or/an d certain physiologica l strategies. It must be emphasized that this method be longs to the category of steepest descent. T herefore, the optimal trajec tor y is most sensitive to the se lected c riterion, i.e. the joint ang les were calculated in such a way that the increments of joints corre spond to the gradient ( partial derivatives) of the criter ion with respect to the joint angles. The gr adient is the minimal chang e in joint angles that results in unit cha nge of the criterion. This method is consistent with that propose d by Hinton in that eac h joint is moved autonomously, in proportion to how muc h moving that joint alone affects the end-effector-targe t distance. [30] In addition to using various criteria, any other optimiz ation methods could be used and certainly different optimal trajectories for both end-effector and joint angles can b e obtained. Trust r egion method [31], for instance, ca n be used to obtain the optimal trajectorie s for the posture comfort hypothe sis if the joint’s range of motion is known. A postural comfor t h ypo thesis predicts that joint excursions in multi-joint tasks are in pa rt determined by joint comfort [32-34]. Cruse (1986) pointed out tha t each joint has an associated disc omfort function and that comfort costs influen ce the movement strategy chosen. The discomfort associated with an individual joint is hig hest near the joint’s biomechanica l range limits and lowest for so me optimal configuration, which te nds to be near the middle of the joint’s r ange of motion [35]. Because certain joints (e.g. knee joint) a re near their biomecha nical range limits, particularly for reaches to the high targe t, it is reasonable to assume that the T rust region method would not provide reasonable results for the se tasks. - 21 - Criteria For full-body voluntary movements, at least two major performance criteria ne ed to be accounte d for, i.e. the end-effector must reach the tar get at the end of movement; and the whole body must maintain balance during the movement (the motion of COM is within the base of suppor t). Other performance criteria , such as minimization of energy , ma y be necessary, esp eciall y f or reaches to the low target. Dif ferent performance criteria may need to be adopted f or different target reaching tasks or even for different periods of the reaching movements [26]. Clear ly that more than one performance criterion is req uired to reasonably predict whole bod y reaching tasks. In fact, Ferry et al [26] had sugg ested that even for a simple arm raising task more than one p erformance criterion may need to be a dopted. Additionally, Par nianpour et al [6] reports that even for a movement task consisting of one segmen t with one degree of freedom may re quire more than one performance criterion. F or end-effec tor path planning, the minimum motor task er ror may not be a necessary performa nce criteria because the location of the end-effec tor at each moment of time is known [29, 36]. The final boun dary condition in Carte sian space also can be put into optimization algorithm constr aints. However , the final posture (in joint space) may be required to be known a priori [28]. In principle, path planning shouldn’t belong in the voluntary movement category. In path pla nning, the path in Cartesian space is known; the question is how to solve a redundant inverse kinematics pr oblem. Whereas in voluntary movement the pa th is unknown, there is no pr ior knowledge about the movemen t. In fact, all constraints ca n be combined into the movemen t performance criter ia. For instance, the joint range can be measured with joint comfort and then combine d into the criteria [37]. - 22 - Conclusions In summary, the p roposed method in this paper is e ffective for the simulation of larg e-scale human skeleton sy stems, which can reasonable predict whole body reaching moveme nts (i.e. final postures, move ment of COM, joint power, and end- effector trajectories etc.). As applied, a c ombination of several control strateg ies such as minimizing end-ef fector error, joint power and CO M and using the simple alge braic calculations of inverse dy namics and forward kinema tics provided good fits to the experimenta l data. In the future dif ferent cost criteria shou ld be examined and compared with e ven more complex movement tasks to further e lucidate how the CNS plans and execute s movements in a kinematic s redundant system. Competing inter est s The authors decla re that they have no competing interests. Authors' contri butions DS developed the 3D simulation models and contributed substantially to the wr iting of this manuscript. JST was responsible for the concept of the reac hing task, data collection a nd anal y sis of experimental data . JST also contributed to the writing of the manuscript. JST was f unded by awards from The National Institutes of Hea lth and Ohio University Post-Doctora l Fellowship Program. - 23 - Acknowl edgement s This work was su pported by National Institutes of H ealth Grant grants RO1-HD04 5512 and from Ohio Unive rsit y Post-D octoral Fellowship Progra m Award. We would like to thank Nicole Van der Wiele and Stacey Moenter for their assistance in data collec tion. References 1. Todorov, E. an d M.I. Jordan, Optima l feedback control as a the ory of m otor coordina tion. Nat Neur osci, 2002. 5 (11): p. 1226-35. 2. Guigon, E., P. Ba raduc, and M. Desmurget, Co mputational Mo tor Contro l: Redund ancy a nd Invariance. J Neurophysiol, 2007. 97 (1): p. 331- 347. 3. Flash, T. an d N. Hogan, The coordin ation of arm moveme nts: an e xperimentally confirme d mathe matical mode l. J Neurosci, 1985. 7 : p. 1688- 1703. 4. Rosenbaum, D.A., et al., Planning re aches b y evaluating store d postures . Psycho l. Rev., 1985. 102 : p. 28-67. 5. Anderson, F .C. and M.G. Pandy, A Dyna mic Optimizatio n Solution for V ertic al Jumpin g in Three Dimensions. Comput Me thods Biome ch Biomed Engin, 1999. 2 (3): p. 201-231. 6. Parnianpour, A ., et al., Computation method fo r simula tion of trunk motion: towards a theoretic al based quantitative assessme nt of trunk pe rformanc e. Biomedical Engineering – A pplications, Basis & Communica tions, 1999. 11 (1): p. 27-39. 7. Massion, J., et al., Is the erect posture in m icrogravity based on the control of tru nk orientation or cente r of mass p osition? E xp Brain Res, 1997. 114 (2): p. 384-9. 8. Nakano, E., et al., Quantitative exam inations of inte rnal repre sentations fo r arm trajector y planning: minimum comm anded torque change m odel. J Neurophy siol, 1999. 81 (5): p. 2140- 55. 9. Menegaldo, L .L., A. de Toledo Fleury, and H.I. Weber, A 'cheap' op timal contro l approac h to estima te muscle fo rces in m usculoskeletal sy stems. J Biomec h, 2006. 39 (10): p. 1787-95. 10. P andy, M .G., F.C. Anderson, and D.G. H ull, A paramete r optimization approach for the optim al control o f large- scale m usculoskele tal systems. J Biomec h Eng, 1992. 114 (4): p. 450-60. 11. P andy, M .G., Compute r modeling an d simulation o f human m ovement. Annu Rev Biomed Eng, 20 01. 3 : p. 245-73. 12. Nagurka, M.L. and V. Yen, Four ier-base d optimal contro l of nonlinea r dyn amic system s. ASME Journal of Dynamic S y stems, Measurements, and Control, 1990. 112 : p. 17-26. 13. Ko, H. and N.I. Badler, Animating Human Locomotion with Inve rse Dyn amics. IEEE Computer Grap hics and Application, 1996. 16 (2) : p. 50-59. 14. SimMechan ics User's Guid e . Version 2 e d. 2004, Natick, MA: The MathWor ks, Inc. 15. Kim, J .H., et al., Pre diction and analy sis of human motion dy namics p erforming various ta sk. Interna ltional J ournal Human Factors Modelling and Simulation, 2006. 1 (1): p. 69-94. 16. Greene, W.B. and J.D. Hechman, The Clinica l Meaurem ent of Jo int mmotion . 1994: American Academy of Orthopaedic Surgeo ns. - 24 - 17. Amankwah, K., R.J. Triolo, and R. Kirsch, Effects of spin al cord injur y on lo wer-limb passive joint mome nts reveale d through a nonlinear vis coelastic model. J Rehabil Res Dev, 2004. 41 (1) : p. 15-32. 18. De Jager, M.K.J., Mathema tical Hea d-Neck M odel for Acc eleration Impacts . 1996, Eindhoven Unive rsity of Technology. 19. Gomi, H. and R. Osu, Task-de pendent v iscoelasticity of human m ultijoint a rm and its spatial ch aracteristic s for inte raction with environm ents. J Neurosci, 199 8. 18 (21): p. 8965-78. 20. Morone y, S.P., et a l., Load-displa cemen t properties of lo wer cer vical spine m otion segme nts. J Biomech, 1988. 21 (9): p. 769- 79. 21. P lage nhoef, S., F.G. Evans, and T. Abdelnour , Anatomical D ata for Ana lyzing Hum an Motion. E xercise and Sport, 1983. 54 (2 ): p. 169-178. 22. Zatsiorsky, V.M., Kinetics o f human M otion . 2002: Human Kinetics. 23. Thomas, J.S . and G.E. Gi bson, Coordina tion and t iming of spin e and hip joints durin g full body reaching tasks. Hum Mov Sci, 20 07. 26 (1): p. 124-40. 24. P ress, W.H., et a l., Numerica l Recipe s in FOR TRAN , in The Art of Sc ientific Computing . 1992: New York: Ca mbridge. 25. Thomas, J.S ., D.M. Corc os, and Z . Hasan, Kinematic a nd kinetic constraints on arm, trunk, an d leg seg ments in targe t-reach ing movemen ts. J Neurophysiol, 2005. 93 ( 1): p. 352-64. 26. Ferr y, M., e t al., Balance con trol during an arm raising movem ent in bipeda l stance: which bio mechan ical factor is c ontrolled? B iol Cy bern, 2004. 91 (2): p. 104-14. 27. Flash, T., The organ ization of hum an arm tra jectory c ontrol , in multip le musc le system s: Biomec hanics and mo vement or ganization , J.M. Winters a nd S.L.-Y. Woo, Editors. 1990, Spring er-Verlag: New York. p. 282-301. 28. P ark, W., et al., A computer a lgorithm for re presen ting spatial-te mporal st ructure of human m otion and a motion ge neralization method. J Biome ch, 2005. 38 (11): p. 2321-9. 29. Abdel-Malek, K., et al., Optimizatio n-based trajectory planning of the huma n upper body. Robotic a, 2006. 30. Hinton, G.E., Parallel co mputations fo r controllin g an arm . Journal of Motor Behavior, 1984. 16 : p. 1 71-194. 31. Coleman, T.F. and Y. Li, An inte rior, trust reg ion app roach for nonlinear minimizatio n subject to boun ds. SIAM Journal on op timization, 1996. 6 : p. 418- 445. 32. Cruse, H., Constraints fo r joint angle control o f the hu man arm. Biological Cybernetics, 1986 . 54 : p. 125-132. 33. Cruse, H., et al., On the cost functions for th e cont rol of the h uman arm move ment. Biol Cybern, 1990 . 62 (6): p. 519-28. 34. Rosenbaum, D.A., et al., Plans for object m anipulation , in Atte ntion and p erform ance XIV: Synergies in experime ntal psych ology, artific ial intellige nce, and cognitive neurosc ience , D.E. Mey er, Kornblum, S., Editor. 1993, MI T Press: Cambridge, MA. p. 803-820. 35. Engelbrecht, S.E., Minimum Principle s in Motor Co ntrol. J Math Psy chol, 2001. 45 (3): p. 497-542. 36. Zhang, X., A.D. Kuo, and D.B. Chaffin, O ptimization-b ased differ ential kine matic modelin g exhibits a velocity- control strate gy for dy namic p osture de termination in seated reaching movements. J Biomech, 199 8. 31 (11): p. 1035-42. 37. Yang, J., et al., Multi-obje ctive op timization for up per body posture p rediction , in the 10th AIA A/ISSM O Multidisciplina ry Analy sis and Optim ization Confe rence . 2004, In proceed ings: Albany, NY. - 25 - - 26 - Figures Figure 1 - Mode l of Optimized Co ntroller Block diagram f or the optimal model of full body reaching movement. The outputs of controller a re the joint angle functions of time. The inver se dynamics and forward kinematics models were used to calculate the phy siological measurements and the loc ation of end- effector. Optimization alg orithm was used to minimize a ce rtain criterion to produce the optimized parameter s for the controller. Figure 2 - End-Ef fector Trajec tories Trajector ies of the end-effector determine d b y 1) the simulation model using two different cost functions, 2) c alculated minimum jerk, and 3) exper imental data averag ed over 15 subjects are plotted for each target height (gray shadow areas represent the standar d error). The left panel is for the high target, middle pane l is for middle target, and the rig ht panel is for low targe t. The trajectories are all remar kabl y similar for ta rget distance, velocity, a nd accelera tion. The largest differences emerge for jerk (bottom row) where the experimental data are not we ll fit b y a ny of the models. The top row indica tes the distance between the end-effec tor and target location. Figure 3 - End-Ef fector Pats The path of the end-effector determined by the d ifferent simulation methods and the experimental data (as described above) for each plane a nd target height are plotted. The paths for plotted for A . sagittal plane all target heig hts B . frontal plane all target heights C . transverse plane for high target D . transverse plane for middle target E. transverse pla ne for low target. - 27 - Figure 4 - Center -of-Mass Movement Simulation and experime ntal motion of whole body COM in horizontal plane for high (top panel), middle ( middle panel) and low target ( bottom panel) respectively. Gray shad ow areas represent bi-dire ctional standard error. Large va riations of AP displacements ar e shown within subjects. Figure 5 - Posture Adopted at Targe t Contact Comparison of fina l postures between simulations with dif ferent cost functions (up per panel) and between simulation ( min power & CO M) and observed posture (lower pane l). Left panel is for high targe t, middle panel is for middle target, and right panel is for low targe t. Figure 6 – Tota l Joint Power Comparison of the tota l joint power (sum of a bsolute each joint power) from simulation with two kinds of crite ria for high (top panel), middle (middle pan el) and low target (bottom panel) respec tivel y . - 28 - T ables T able 1 - Inp ut Joint Dat a Range of motion an d viscoelastic coeffic ients of joints. Limitation (deg ) e Viscoe lastic coefficients Joint Plane upper lower Stiffness K (N.m.deg -1 ) Damper B (N.m.deg -1 .s) Flexion/extension 54.3 -12.2 1/6 a Int/external Rota tion 0.01 -0.01 Ankle Add/abduction 19.2 -19.2 1/15 a Flexion/extension 141.2 -0.01 1/20 a Int/external Rota tion 0.01 -0.01 Knee Add/abduction 0.01 -0.01 Flexion/extension 12.1 -121.3 1/3 a Int/external Rota tion 44.2 -44.2 Hip Add/abduction 25.6 -25.6 1 a Flexion/extension 62 -167 0.192 d 0.014 d Int/external Rota tion 69 -104 0.192 d 0.014 d L Shoulder Add/abduction 184 -0.01 0.192 d 0.014 d Flexion/extension 0.3 -140.5 0.1571 d 0.0122 d Int/external Rota tion 81.1 -75 0.1571 d 0.0122 d L Elbow Add/abduction 0.01 -0.01 Flexion/extension 35.3 -21.1 0.1047 d 0.0105 d Int/external Rota tion 0.01 -0.01 0.1047 d 0.0105 d L Wrist Add/abduction 74 -74.8 0.1047 d 0.0105 d Flexion/extension 141 -141 0.25 b Int/external Rota tion 93 -93 0.42 b Cervical Add/abduction 172 -172 0.33 b Flexion/extension 27 -27 0.25 c Int/external Rota tion 21 -21 0.42 c Thorax Add/abduction 4 -4 0.33 c Flexion/extension 43 -43 0.25 c Int/external Rota tion 19 -19 0.42 c Lumbar Add/abduction 8 -8 0.33 c Flexion/extension 62 -167 0.192 d 0.014 d Int/external Rota tion 69 -104 0.192 d 0.014 d R Shoulder Add/abduction 0.01 -184 0.192 d 0.014 d Flexion/extension 0.3 -140.5 0.1571 d 0.0122 d Int/external Rota tion 75 -81.1 0.1571 d 0.0122 d R Elbow Add/abduction 0.01 -0.01 Flexion/extension 35.3 -21.1 0.1047 d 0.0105 d Int/external Rota tion 0.01 -0.01 0.1047 d 0.0105 d R Wrist Add/abduction 74.8 -74 0.1047 d 0.0105 d a calculate d from Amankwah et al . 2004. b adopted fro m De Jager, 1996. c adopted fro m Moroney et al 1988. d adopted fro m Gomi and Osu 1998. - 29 - e adopted from Greene and Hechman 1994 Table 2 - Compu tational Costs Costs of optimal control strategies and CPU time s for each target loca tion are shown below. Optimal control stra teg y Target location Measureme nt min Error min Power min COM min Power & COM Total power sq uared ( ) ( ) dt f t T ∫ × × 0 θ τ θ τ & & (J 2 ) 3436 320.2 285.7 130 Final COM squar ed Cf T Cf x x (m 2 ) 0.03509 0.01309 0.000453 0.00784 High CPU time (hour) 1.84h a 3.90h b 1.96h b 5.43h a Total power sq uared 5608 1023 1672 353.4 Final COM squar ed 0.05101 0.02597 0.00047 0.003058 Middle CPU time 1.13h b 2.37 h b 2.66 h b 2.73h b Total power sq uared 11167 8638 13940 3705 Final COM squar ed 0.09823 0.0917 0.00155 0.01258 Low CPU time 2.33h b 4.76 h a 2.99 h b 3.11h a a Desktop compute r: Intel Xeon, 3.20G Hz and 3.19GHz, 2GB of RAM, Windows XP, Matlab 2006b b L aptop computer: Intel Pentium M, 1.86GHz, 1 GB of RAM, Windows X P, Matlab 2006b - 30 - Table 3 - COM compa risons Final COM displac ements (mm) derived from simulation mode ls are compared to experimental resu lts for each target location. Me an values from experimen tal data (± SEM) are also prese nted. H igh target Middle target Low target Anterior-posterior 172.3* 219.7* 311.2* Min task error Mediolatera l 73.6 * 52.5* 37.3 Anterior-posterior 113.6* 160.3* 299* Min Power Mediolatera l 13.3* 16.8* 47.9* Anterior-posterior 2* 17.9* 39.2* Min Com Mediolatera l 21.2* 12.3* 3.7(ns) Anterior-posterior 49.1 (ns) 60.2 (ns) 112 (ns) Min Power & COM Mediolatera l 2.1 (ns) 5 (ns) 5 (ns) Anterior-posterior 58.7±17. 9 67.8±20. 8 86.6±21.2 Experiment Mediolatera l -8.5± 6.7 -7.1± 6.8 3.5±7.1 * indicates p<.0 5 - 31 - Table 3 - Total energy comparisons Comparison of total e nergy (J) between model simula tions and their corresponding experimental resu lts of all joints are shown for each targ et location. Mean values for experimental data (± SEM) a re also presented. High target Middle target L ow target min Error 52.9 8* 61.64 (ns) 115.9 * min Power 21.06 * 35.46 * 104.5 * min COM 29.03 * 56.24 * 155.4 * min Power & COM 15.21 * 32.59 * 90.38 * Experiment 27.6±12. 44 50.86± 21.47 98.75±61 .33 Figure 1 End - effector error Criteria Forward Kinematics Polynomial Target location Parameters End - effector location + - Inverse Joint angles Optimization Joint power, C O M , e t c . Skeleton Initial Figure 1 Figure 2 Figure 2 Figure 3 Figure 3 Figure 4 Figure 4 Figure 5 Figure 5 Figure 6 Figure 6