An Optimization-Based Differential Inverse Kinematics Approach for Modeling Three-Dimensional Dynamic Postures During Seated Reaching Movements Xudong Zhang Don B. Chaffin Department of Industrial & Operations Engineering University of Michigan Ann Arbor, MI 48109 Technical Report 96-15 November 1996

AN OPTIMIZATION-BASED DIFFERENTIAL INVERSE KINEMATICS APPROACH FOR MODELING THREE-DIMENSIONAL DYNAMIC POSTURES DURING SEATED REACHING MOVEMENTS TECHNICAL REPORT 96-15 Xudong Zhang and Don B. Chaffin Center for Ergonomics Department of Industrial and Operations Engineering The University of Michigan Ann Arbor, Michigan 48109-2117 November 1996

AN OPTIMIZATION-BASED DIFFERENTIAL INVERSE KINEMATICS APPROACH FOR MODELING THREE-DIMENSIONAL DYNAMIC POSTURES DURING SEATED REACHING MOVEMENTS Abstract Although optimization techniques have been extensively utilized to model human postures and movements, existing models are either static, or dynamic but restricted in the scale or dimension of the biomechanical system due to prohibitive optimization algorithms. An optimization-based differential inverse kinematics approach is described in this paper. The posture determination problem is viewed in the velocity domain as a problem of determining how a change of three-dimensional hand position is partitioned by changes in the excessive degrees of freedom. The inherent kinematic redundancy is resolved by the weighted pseudoinverse of the Jacobian which quantifies the mapping from the joint angle changes to the hand position changes. A set of weighting parameters visualizes the strategy of allocating energy-type effort among the body segments involved in a movement. An optimization procedure based on simulated annealing is employed to identify the parameter values such that the predicted profiles best emulate the observed profiles. The proposed approach is tested empirically by modeling three types of seated reaching movements distinguished by the hand motion direction. The results demonstrate that accurate representations of measured motion profiles can be achieved, and reveal some difference between and congruence within sets of parameters corresponding to different types of movements. Theoretical implications regarding pre-organized motor program, levels of movement coordination, and dynamic posture control mechanism are discussed. 1

2 Introduction Optimization has been a viable tool to resolve redundant problems pertinent to biomechanical modeling (Chao and An, 1978; Crowninshield and Brand, 1981). Two primary redundancies encountered in modeling musculoskeletal systems, static or dynamic, are kinematic redundancy and muscle force redundancy. They are both incurred by the same ill-posed situation: the number of unknowns, either some kinematic variables or muscle forces, exceeds the number of equations available. Optimization techniques generally presume an optimal strategy adopted by human beings in controlling postural kinematics or muscle forces, and, by mathematically modeling this strategy as an objective function, make the indeterminate equation system solvable. Kinematic redundancy occurs in posture modeling when a human being is represented as a multi-segment rigid linkage system with excessive degrees of freedom (DOF). The problem of resolving indeterminacy becomes more sizable and complex as a movement that consists of a succession or sequence of instantaneous postures is being modeled. For example, if every instant corresponds to one optimization routine, as the modeling is extended from a single posture to a movement, the scale of the overall problem is literally multiplied by the number of time frames composing the sequence. This complexity has effectively been one of the major factors accounting for the scarcity of optimization-based predictive models to describe normal human postures and movements. In contrast, the application of optimization techniques to estimating muscle forces has been a rather fruitful research domain-a considerable volume of predictive models have been generated and, to a certain extent, validated (Seireg and Arvikar, 1975; Chao and An, 1978; Hardt, 1978; Crowninshield and Brand, 1981; Bean et al. 1988; Hughes, 1991). The few existing optimization-based models for posture prediction have inherent limitations that have hampered their applications. Many of these models possessed at least one of the following shortcomings: they were either static, or two-dimensional (within the sagittal plane), or incorporated a very limited number of segments or degrees of freedom

3 (Park, 1974; Ayoub et al., 1974; Ayoub and Hsiang, 1992; Byun, 1992; Jung et al., 1995). Three-dimensional upper body posture models include the Boeman model (Ryan, 1970), the ones by Kilpatrick (1970) and, more recently by Jung et al. (1994). All relied on iterative optimiation routines, one for each time frame, to predict three-dimensional sequential static postures of the torso and the arms during seated operations. None of these modeling efforts, however, recognized two critical facts: first, it is not guaranteed that the sequences formed by conjoning discretely optimized static postures would be as smooth or "graceful" as real human motions; and second, there is a distinction, conceptually as well as physically, between a static posture and an instantaneous dynamic posture (Zhang and Chaffin, 1996). In addition, the objective functions of optimization procedures employed in these models were based on heuristics and subject to empirical adjustments. Interpretations of these objective functions, in physically meaningful terms, were difficult. Therefore, despite the fact that good modeling accuracy might have been ultimately achieved by these models, little insight was gained into the optimal strategy presumably used in human postural control. In dynamic models, the large amount of computation involved in dynamic optimization can be a prohibitive factor (Girard, 1991). For instance, an optimizationbased dynamic model was developed by Ayoub et al. (1979) to predict the path of simple upper extremity motions. A similar approach was employed (Ayoub and Hsiang, 1992) to produce a two-dimensional dynamic biomechanical simulation of material handling motions that were assumed to be planar and symmetric. These two studies both used and advocated dynamic programming as the mathematical tool for solution search. Due to the "curse of dimensionality" problem associated with dynamic programming (Bellman, 1957), the number of variables or states that can be incorporated in the model is constrained. As a consequence, either the dimension (2D or 3D) of the model or the number of body segments included in the model, will have to be compromised.

4 There is another class of research, primarily related to gait analysis and simulation, which deals with both kinematic redundancy and muscle redundancy (Hardt, 1978; Pandy et al., 1992; Yamaguchi et al., 1995; Zajac and Gordon, 1989). In these investigations, dynamics equations were established to represent the causal relationships between muscle action and motion generation; optimization was applied, quasi-statically or dynamically, with a muscle-stress-related objective function. However, regardless of the scale of the musculoskeletal system being modeled, the number of muscles always far exceeds the degrees of freedom of the system; for instance, there were 30 muscles but only 5 degrees of freedom in the study by Yamaguchi et al. (1995). Thus, these studies solved a muscle redundancy problem per se, and the motions of only a limited number of segments were the "by-products" of musculotendon dynamics. As this type of study still suffers from the dimensionality problem (Yamaguchi, 1990), the number of body segments included is usually limited. In this paper, an optimization-based differential inverse kinematics approach for modeling three-dimensional dynamic postures during right-handed seated reaching is presented. The biomechanical system concerned is represented by a linkage composed of four rigid segments: the torso, the right clavicle, the right upper arm and forearm. Seven degrees of freedom are incorporated in this linkage. The approach is differential in that the problem of multi-segment posture determination is solved in the velocity domain-it is viewed as a problem of determining how a change of 3D end-effector (the hand) position is partitioned among the involved, normally excessive, degrees of freedom. The mapping from changes in seven joint angles, each measuring one degree of freedom, to the changes in three hand coordinates, is mathematically described as a 3x7 Jacobian matrix. The weighted pseudoinverse of the Jacobian resolves the kinematic redundancy with an implicit objective of minimizing the sum of the weighted norm of angular changes. A set of weighting parameters therefore models, in a straightforward and parsimonious way, the strategy of allocating energy-type effort among participating segments or joints during a

5 movement. Using a simulated annealing (SA) search method, sets of parameters for different types of movements performed by different individuals are determined, such that the predicted motion profiles best approximate the actual profiles. In doing that, inherent strategies hypothetically used by people in dynamic posture control can be quantitatively assessed and visualized. Methods Modeling Scheme Linkage Representation One of the most commonly used biomechanical models of the human musculoskeletal system (or some portion thereof) is a collection of interconnected rigid segments (Andrews, 1995). In this study, the system used to describe seated postures during reaching movements consists of the torso and the right upper extremity. As depicted in Figure 1, the biomechanical system concerned can be adequately represented by a chain-like linkage composed of four rigid segments: torso, right clavicle, right upper arm, and right forearm. Right Acr t Clavicle Suprasternale Right Wrist / F or1earm &To & Hi Forearm & Hand wight Elb Bottom of Spine Figure 1. A linkage representation composed of 4 segments: torso, right clavicle, right upper arm, and right forearm.

6 This linkage is constructed by connecting five hypothetical spherical or revolute joints: L5/S1, suprasternale, right shoulder (acromion), right elbow, and right wrist. Implicit to this representation are several simplifying assumptions: 1) the reaching movements addressed in this study are right-handed, with the left hand and arms not posing any restraint to the movement; 2) the left hand and arm are therefore not modeled, as the study is limited to single-target-driven reaching; 3) finger actions and hand configurations are also excluded from the model due to the excessive degrees of freedom and complexity with which they are associated; and 4) the right hand is sustained in a neutral, "normative" posture and can be considered as an extended part of forearm. Seven degrees of freedom are incorporated in the linkage to specify the major motions of the system: 2 DOF at the bottom of spine (L5/S 1) for torso flexion and lateral bending, 1 DOF at the suprasternale for clavicle rotation (or torso twisting), 2 DOF at the acromion for shoulder abduction and extension, and 2 DOF at the elbow for humeral rotation and elbow flexion. The seven degrees of freedom are measured by seven Euler angles. The definitions of these seven Euler angles are established with reference to five coordinate systems as illustrated in Figure 2: four local or segmental coordinate systems affixed to individual links and one base coordinate system as the general reference frame. Table 1 presents the nomenclature of these seven Euler angles. These names are provided here for the ease of physical interpretation, and since they are Euler angles, they may not comply with clinically or anatomically used angle conventions. Note that since the axial rotation of each link cannot be specified, torso axial rotation is modeled partly by a rotation of the clavicle with respect to the torso long axis. Similarly, in order to identify a possible change of forearm orientation caused by humeral rotation an extra DOF is modeled at the elbow which otherwise could be well represented by a 1-DOF revolute joint. Of further note is that as the suprasternale is modeled as an imaginary 1-DOF revolute joint to provide clavicle rotation, the angle subtended by the clavicle and the torso is assumed to be fixed (identified as a in Figure 2).

7 Z2 Z1 Z4 Y4f 14(forearm & '" 2erY1 13 (upper arm) 11(torso) X3 ZO xo t/ general reference frame YO Figure 2. Coordinate systems established to derive the seven Euler angles which in conjunction with fixed link parameters describe the configuration of the linkage. Table 1. The nomenclature of seven Euler angles which determine the configuration of the linkage Angle Definition 1: Torso Flexion Rotation angle of the torso link with respect to X0 axis. 02: Torso Lateral Bending Rotation angle of the torso link with respect to YO axis. 03: Torso Twisting Rotation angle of the clavicle link with respect to Z1 axis. 04: Shoulder Extension Rotation angle of the upper arm link with respect to X2 axis. 05: Shoulder Abduction Rotation angle of the upper arm link with respect to Y2 axis. 06: Humeral Rotation Rotation angle of the forearm link with respect to Y3 axis. 07: Elbow Flexion Rotation angle of the forearm link with respect to Z3 axis.

8 Note: All the angles are 0 when a standard anatomical posture is assumed. Sign convention follows the right hand rule. Differential Kinematics and Jacobian The end-effector (the hand) position with respect to the general coordinate system can be expressed in terms of joint variables as the seven Euler angles, and link parameters including link lengths and link offsets (Denavit and Hartenberg, 1955). It is derived by concatenating a series of transformation matrices between two neighboring coordinate systems. The more detailed mathematical derivations are provided in Appendix A. Let r = [x y z]T be the hand position in three dimensions, then the non-linear, complex relation between r and the Euler angles can be expressed in an abstract form as: r=[fl(01,...,07) f2(01,.,07) f3(O...,07)]T =f(o) (1) In fact, equation (1) represents the kinematic relation dealt with in a static posture prediction model (e.g., Ryan, 1970; Jung, 1994) which attempts to determine the closedform non-linear solutions of joint angles from the hand location, strictly in displacement domain. Differentiating (1) with respect to time results in r =(0)0 = J(O)0 (2) where J is the Jacobian of f with respect to 0. In this particular case, the J is a 3x7 matrix with each element Jij =-i (i=1,2,3; j=1...7). (3) dOj What equation (2) represents is a linear relationship between the end-effector velocity and the joint angular velocity. This method of coping with velocity or rate rather than position to take advantage of the linearity is termed differential kinematics (Nakamura, 1991) and was first introduced by Whitney (1969). The Jacobian J which quantifies the mapping from the changes in joint angular coordinates to the changes in the hand coordinates plays a central role in this method. Optimization-Based Weighted Pseudoinverse

9 The problem of posture prediction using differential kinematics now can be stated as: given a change in hand position, how to determine the changes in joint angles. Mathematically, it requires solving equation (2) for O by attaining the inverse of the Jacobian. For redundant systems such as the one concerned in this study, the ordinary inverse of J is not defined. However, by hypothesizing that an objective function to minimize c= lwolI (4) is followed in solving redundancy, 0 can be determined as (=W [-JW-r.1 (5) Here, WO is a symmetric and positive definite 7x7 matrix, and # symbolizes the pseudoinverse of a matrix. The general form of a finite number of solutions that satisfy (2) is actually given by e=W' JW ' ] r+ Wo (I -[JWo JJWo')V (6) where I is a 7x7 identity matrix and V is an arbitrary vector of 7 dimensions. Among the infinite number of solutions, a unique solution with a minimum weighted-norm is provided by equation (5). Note that the minimization of (4) is equivalent to the minimization of 2 oTw (7) where W = WOTWo. Term (7) has been referred to as the instantaneous weighted system kinetic energy (Whitney, 1969). The diagonal components of the matrix Wo may be regarded as weighting parameters which model a strategy of allocating kinetic-energy-type effort among involved joints so as to minimize the weighted total effort. Therefore, a relatively smaller value of the parameter would indicate that the corresponding joint angle is more involved in the motion, and a "heavier" weight signified by a greater value would suppress or "penalize" any change in the associated joint angle. Empirical Testing Problem of Parameter Determination

10 In order to accurately predict seated reaching movements using the proposed modeling scheme, the parameters that characterize the strategy presumably used by people during a movement need to be identified. It is theorized that if the profiles predicted by a particular set of parameters, hopefully time-invariant, closely emulate the actual profiles of a movement, this set of parameters then describes, in a parsimonious and quantitative way, the inherent strategy. This suggests a technique which empirically determines the parameters that provide the best-fit with respect to the real human motion data. Such a technique involves the resolution of another optimization problem described as follows. Since motion data are captured by 3D measurement systems in discrete time frames at a certain sampling rate, kinematic quantities, hereinafter, will all be expressed in discrete forms with the subscript being the time frame number. Let et and Ot, both 7-dimensional vectors, represent respectively the instantaneous angles and angular velocities at an instant of time t. Equation (4) can now be re-written as: d A t t tew [j i -I]AO (8) At= At = [J,_wo \ ^ (8) Note that J,_l is in fact Jt_- (t-_1), a function of instantaneous angles at t-1 only. By factoring and eliminating the finite sampling time interval At at both sides, equation (8) becomes 0, - ) = WW]Art1. (9) {3t {3t-1 =WO[J-tWO t]#/, (9) Using equation (8) recursively, 0, can be computed as 0, = 01 + W [JkW ] Ark. (10) k=l In equation (10), as 01, Jk, and Ark are known or derivable, the seven weighting parameters in WO, presumed to be time-invariant, become the only variables influencing the behavior of 0,. Thus, if 0e represents the measured set of joint angles at time t, WO can be determined by minimizing T t- t Eie1 +Xw0 [JkWo Ark -O' (11) t=2 k=l

11 where T is the number of time frames used to completely describe a movement. Equation (11) now provides a means to evaluate the global average (across the entire movement) of the mean absolute error of 7 predicted angles with respect to the measured angles. Two alternative criteria also are used as the objective functions of optimization procedure to determine the parameters: Wo1 J-WO1] Art- (12) t=2 and T 01 + WO l[Jt-W-1w ] Art -O * (13) t=2 Function (12) evaluates the global average (across the entire movement) of the mean absolute error of the angular velocity values with respect to the actual values; function (13) evaluates the mean error of the predicted angles of the ending frame, which, due to progressive error accumulation, may be substantial even if (11) or (12) is minimum. There is a computational advantage to the latter two criteria, because the number of mathematical operations involved is considerably reduced as compared with (11). Nevertheless, a more comprehensive performance evaluation of a determined set of parameters (based on one particular criterion) is provided by a "cross-examination" using all the three criteria. The proposed modeling scheme provides the flexibility of implementing different scenarios of the parameters as variables in the optimization procedures. A simplifying yet plausible scenario is to hypothesize a segmental weighting strategy and use four parameters, one for each segment. The Wo can then be specified as w1 0 0 0 0 0 0 0 w10 0 0 0 0 0 0 w2 0 0 0 0 0 0 0 w30 0 0 (14) 0 0 0 0 w30 0 0 0 0 0 0 w4 0 0 0 0 0 0 0 w4

12 Note that angles of same segment have identical designate weight. A test of whether this Wo configuration or a more heterogeneous one results in a set of parameters with satisfactory predictions, based on above evaluation criteria, may lend support to different theories of how the control of and coordination between segments are actualized. Simulated Annealing The optimization problem of determining the best-fitting parameters, as posed by (1 1)-( 13), is difficult to solve. It is highly non-linear and a general knowledge of the function surface is not available. Gradient-based non-linear optimization methods are not applicable because the derivative of the objective function is analytically not achievable. Simulated annealing (SA) offers a heuristic alternative of obtaining good solutions to difficult optimization problems. It is an iterative stochastic search method, analogous to the physical annealing process whereby molten material is gradually cooled so that a minimum energy state is achieved. The simulated annealing approach has the virtue of being a powerful, relatively easy to implement, search technique for optimizing any function without much prior knowledge about its "terrain". Specifically, it is well suited for a class of problems such as the one concerned in this study, because the approach only requires a point's corresponding function value and no other information. A thorough treatment of simulated annealing can be found in Kirkpatrick et al. (1983) or Eglese (1990). Specifications of the simulated annealing algorithm used in this study to determine the weighting parameters are as follows. The three objective functions (11 )-(13) all serve as evaluation functions. The results are "cross-examined": a solution obtained by one objective function is evaluated using the other two functions; overall performance is then assessed to determine the final optimal set of parameters. A neighborhood function which specifies the move from one intermediate solution to the next is given by wi(n) = Wi-l(n- 1)+ Rand[-0.5,+0.5] (15)

13 where n is the iteration number and Rand[-O. 5, +0.5] generates a random number between -0.5 and 0.5. The starting "temperature" is 20 (dimensionless). Heating factor and cooling factor are 2 and 0.99 respectively. The search terminates as soon as either of the two criteria is satisfied: 1) the number of iterations exceeds 2000; 2) the function value is less than the stopping threshold-0.8 degree for functions (1 1) & (13), 0.2 degree/frame for function (12). Experimental Data for Model Testing Data for testing the proposed modeling scheme were acquired from an experiment (described in detail in Zhang & Chaffin, 1996) in which three types of seated reaching were performed by six young adult subjects. These three types of reaching tasks were distinguished by the directions of hand motions (i.e., Ar): anterior-posterior (AP), mediallateral (ML), and up-down (UD). Each type was performed at four different hand trajectory positions which were characterized by two of the following three dimensions: 1) the height-shoulder or hip height; 2) the radial distance from the right shoulder-60% or 120% of the arm length; 3) the offset angle from the sagittal plane-0 degree or 45 degrees, all with respect to an upright seated person. The former two dimensions were adjusted in accordance with the individual anthropometry. Thus, a movement can be described by its hand trajectory, for example, as anterior-posterior at shoulder height with 0 degree offset angle. The 12 movement trials (3 typesx4 positions) were completed by each subject using a normal, self-paced speed. Figure 3 graphically illustrates the experimental setup. The twelve hand trajectories approximately intersect each other at eight locations as labeled in the Figure. These intersection locations can be conveniently used to codify the movements. For instance, an anterior-posterior reach at the shoulder height with 0 degree offset angle was coded as AP13, whereas a medial-lateral reach at the hip height and close-in distance was coded as ML56. Reflective spherical markers were placed over palpable body landmarks identifying the right hand, right elbow, right acromion, suprasternale, right and left anterior superior

14 iliac spine (ASIS). The two ASIS markers were utilized to estimate the bottom of the spine assumed as their bisection. A MacReflexTM motion analysis system with four digital cameras was employed to measure the coordinates of these markers during the movements at a sampling frame rate of 25 Hz. A linkage representation as portrayed earlier in Figure 1 was constructed. Profiles of the seven Euler angles based on the definitions established with respect to the linkage system were derived, and, along with those of the hand coordinates, served as the empirical basis for model testing. Figure 3. An experiment in which three types of seated reaching movements were performed: anterior-posterior (AP), medial-lateral (UD), and up-down (UD). The labeled locations where the hand motion trajectory approximately intersect are used to codify the movements

15 Results Each optimization procedure for modeling one movement yields two results: the optimal objective function value and the corresponding solution. They represent, respectively, the best accuracy achieved by the proposed model and the weighting parameters identified for the particular movement. In general, the model performed well with the use of only four weighting parameters which was based on the hypothesis of a segmental effort-allocating strategy. Table 2 presents some descriptive statistics of the overall accuracy measures in terms of the three evaluation criteria for modeling with four parameters. These measures show considerable variation. In fact, the best achieved accuracy and the corresponding weighting parameters differ distinctively for the three types of reaching movements. Therefore, results are classified by the movement type and are presented in the following. The emphasis is placed on the use of four parameters (called a four-parameter scenario) when presenting these results because of its general success. Attention is also directed to examining whether seven different weighting parameters (called a seven-parameter scenario) provide a significant improvement in modeling accuracy, especially in those conditions where the four-parameter scenario fails to perform well. Table 2. Descriptive statistics of overall (across all the subjects and movements) accuracy measures in terms of three criteria for modeling with four weighting parameters Evaluation Criterion Mean Std. Dev. Minimum Median Maximum Cl (degree) 3.3 2.7 0.6 2.4 11.6 C2: (degree/frame) 0.57 0.38 0.15 0.45 1.9 C3: (degree) 4.8 4.3 0.78 3.5 20 Note: C1: the averaged (across the entire movement) mean error of 7 predicted angle values with respect to the actual angles; C2: the averaged (across the entire movement) mean error of 7 predicted angular velocity values with respect to the actual values; C3: the mean error of the 7 predicted angels of the terminal frame.

16 Modeling Accuracy Figure 4 consists of three box plots of the modeling accuracy that was achieved using the two different weighting parameter scenarios. Each is evaluated by the three criteria for the anterior-posterior (AP) reaching movements. A box plot graphically provides a collection of information about the data: the top, bottom, and line through the middle of the box correspond to 75th percentile (top quartile), 25th percentile (bottom quartile), and 50th percentile (median) respectively; the whiskers on the bottom extend from the 10th percentile and top 90th percentile; the square symbol represents the mean. All three plots in Figure 4 exhibit the same clear pattern: a better modeling was achieved with AP13 and AP57 as compared to AP24 and AP68, regardless of the parameter scenario. Recall that the former two were the motions with a 0-degree offset angle (hand trajectory in the sagittal plane) whereas the latter two were the motions with a 45-degree offset angle. Overall, for anterior-posterior reaching motions, the medians of the accuracy measures in terms of the three evaluation criteria are 5 degree, 0.86 degree/frame, and 7.3 degree respectively. Figure 5 illustrates the predicted versus actual profiles of all the seven angles for a "typical" movement trial with which the accuracy measures were close to the median values. Notice particularly a substantial discrepancy between the predicted and the observed for the elbow flexion which does not appear to be coordinated well with the motions of other joints. A comparison between two scenarios (left vs. right in Figure 4) reveals that, by the use of seven weighting parameters, modeling accuracy was improved as measured with Cl and C3 but not with C2. Figures 6 and 8 present similar box plots for medial-lateral (ML) and up-down (UD) reaching movements. Notice the scales of the ordinates in these box plots are half of those for anterior-posterior movements. A superior accuracy attained by the proposed modeling scheme for these two types of movements: even with the four-parameter scenario, the medians of the accuracy measures in terms of the three evaluation criteria are 1.5 degree, 0.35 degree/frame, 2.0 degree for medial-lateral motions, and 1.47 degree,

17 0.34 degree/frame, 2.06 degree for up-down motions. Again, an illustration of the predicted versus actual angular profiles of a movement, with which a "mediocre" modeling accuracy was achieved by the four-parameter scenario, is presented in Figure 7 for ML motions and Figure 9 for UD motions. In contrast to the AP movements, however, the improvement brought forth by the seven-parameter scenario appears to be trivial. An interesting pattern analogous to the one described above is exhibited by the box plots for the ML movements: ML12 and ML56 out-performed ML34 and ML78, regardless of the parameter scenario. Recall that the former two were the motions at the close-in locations whereas the latter two were the ones at the far-out locations. However, for the UD movements, no such pattern can be identified. Weighting Parameters While the accuracy measures the performance aspect of the model, the weighting parameters are parsed in the interest of making inferences, from the optimization procedure applied, regarding the strategic aspect of the movement or posture control, the consistency or variability across individuals, and so on. As the weighing parameters are scalable (i.e., only the relative ratio matters), the four parameters within each set that corresponds to a particular movement are normalized as percentages of their total to better visualize the partitioning among segments. Figure 10 illustrates these normalized percentages for all the subjects in three area plots for the three types of movements. The an across the subjects for each type of movements is presented as a separate column. These plots show that overall there exist some degree of congruency within each type of movements across the subjects, although individual differences yet appear to be substantial as reflected by the "zigzag" contours. However, the partition among the segment varies distinctively across different types of movements. For instance, the torso is very much involved in the anterior-posterior movements as its average percentage is relatively lower, but virtually motionless during the up-down movement as its

18 percentage is dominantly high; the reverse situation appears to be true for the clavicle movements. Criterion and Computation Considerations As described earlier, a set of parameters as a solution derived from the optimization procedure using one of the criteria was also cross-evaluated by the other two criteria. A final set of parameters was therefore determined from three solutions based on their respectivethree criterion values. In this determination process, it was found that solutions generated by C1 and C3 consistently out-performed the one generated by C2-unless the three solutions converged (only approximately in simulated annealing method), the one that was optimal in terms of C2 usually was not the optimal solution overall. For most cases, Cl and C3 yielded comparable solutions: each scored the best in its own category but not much worse than, if not as well as, the other in the other category. However, there is a considerable computational advantage to the criterion C3. Suppose the W" is given by a set of parameters for a movement of T time frames. Without saving the intermediate et values, the numbers of various matrix or vector operations needed to derive the Cl function T(T-l) T(T-1) value (referred to equation 11) are: ( - pseudoinverses; 3* - ) multiplications, 2 2 T( -1 ) + T - 2 additions, and T- 1 subtractions. For C3 function (referred to equation 13), these numbers are T- 1, 3(T- 1), T- 1, and 1 respectively. Note that, the larger the T is, the more advantageous the C3 is. For example, for a 25 frame motion, it took 200-300 seconds (Mathematica running on a 7100/66 Power Macintosh) to solve for a set of parameters based on C3, but 8-10 times more based on C1. Therefore, C3 can be used singly in the optimization procedures for solution search, and should be preferred when the computation resource or time is of concern.

19 10 ~ 8 0} O 6 > 4 - AP13 AP57 AP24 AP68 0 degree offset 45 degree offset 4 parameter scenario AP13 AP57 AP24 AP68 0 degree offset 45 degree offset 7 parameter scenario 1.8 -1.6 -E 2 1.4 -| 1.2 -0) 0 1 -Q 0.8 -0.6 -0.4 0.2 -0 -20-. 18 16 14 12 '0 010 X 8 C 6 -4 -2 -n\ 0 1 4&- LT-i 0 — r- y I I I AP13 AP57 AP24 AP68 0 degree offset 45 degree offset 4 parameter scenario AP13 AP57 AP24 AP68 0 degree offset 45 degree offset 7 parameter scenario j, i~~~~ I E!3 0 41 L77J AP13 AP57 AP24 AP68 AP13 AP57 AP24 AP68 0 degree offset 45 degree offset 0 degree offset 45 degree offset 4 parameter scenario 7 parameter scenario Figure 4. Modeling accuracy measured by three evaluation criteria for anterior-posterior (AP) reaching movements at four different locations: a four-parameter scenario versus a seven-parameter scenario.

20 la c: 0 1 0O u o 0 0.0 o 00 5, -10 -20.. -- -30 -40 -50 -60 0 5 10 15 20 25 time frame -18 -- -- -20 -22 -24 -28 0 5 10 15 20 25 time frame -326 -_.__ 0 5 10 15 20 25 time frame 55 50 422 4C 30 3/ 20 0 5 10 15 20 25 time frame I II I I 41 It II a 41 0 14.14 4 11.1 20 15 lo 0 5_____ 0 5 10 15 20 25 time frame 120 1 00 C -I C 40 " 20 0 5 10 15 20 25 time frame 37.5 --- - - * I zD 35 --- DI - o 32.5 2 30 27.5 N '"v 5 10 15 20 25 time frame di 0 c,) Figure 5. Predicted (solid) versus actual (dotted) profiles of seven angles for an anterior-posterior motion trial with C1=5.25 deg, C2=0.85 deg/frame, and C3=7.26 deg.

21 5 L_ 4 -Ia 0) 3 ->2 1 -0 0.8 -E 0.6 -0.4 -" 0.2 - U.I I I I I I I I ML12 ML56 ML34 ML78 ML12 ML56 ML34 ML78 close-in far-out close-in far-out 4 parameter scenario 7 parameter scenario I I m m m 0 t:; "7-j L-r-j ML12 ML56 ML34 ML78 ML12 ML56 ML34 ML78 close-in far-out close-in far-out 4 parameter scenario 7 parameter scenario 8 -0 - w 6 -4 2 - I I m -0 -I L- __j 0 9i] Eip v I I I I I I ML12 ML56 ML34 ML78 ML12 ML56 ML34 ML78 close-in far-out close-in far-out 4 parameter scenario 7 parameter scenario Figure 6. Modeling accuracy measured by three evaluation criteria for medial-lateral (ML) reaching movements at four different locations: a four-parameter scenario versus a seven-parameter scenario.

22 I IC 0 C o To i -22 -24 -26 -28 -30........ 0 5 10 15 20 time frame -17 -17.5-..... -18 -18.5 -19 ' ' 0 5 10 15 20 time frame -10 -12 -14 -16 -18 0 5 10 15 20 time frame 82.5 80 _ 77.5 75 72.5 67.5 0 5 10 15 20 Do la 0 a 0 1. a 2 C to 0.C 0 5 10 15 20 time frame 50 45 40 8. -- 35 ~ 0 5 10 15 20 time frame 40. 31 20 10 -10 -20 0 5 10 15 20 time frame a '. 0 u I. o a 0 40 la 0.0 "S time frame Figure 7. Predicted (solid) versus actual (dotted) profiles of seven angles for a medial-lateral motion trial with C1=1.62 deg, C2=0.35 deg/frame, and C3=1.53 deg.

23 5 -4) L 4 -0) U - '1 13 ->2 U T 6I IUD UD51 UD62 UD73 UD84 UD51 UD62 UD73 UD84 close-in far-out 4 parameter scenario close-in far-out 7 parameter scenario 1 -0.8 -E 0.6 -0) 4) 0.4 -03 'O 80.2 -0, 0 - 1U ^^~~- ll L UD51 UD62 UD73 UD84 UD51 UD62 UD73 UD84 close-in far-out close-in far-out 4 parameter scenario 7 parameter scenario 8 -4) a 6 -'O 4) 0 2 - U a n S ffi I I I I I I I UD51 UD62 UD73 UD84 UD51 UD62 UD73 UD84 close-in far-out close-in far-out 4 parameter scenario 7 parameter scenario Figure 8. Modeling accuracy measured by three evaluation criteria for updown (UD) reaching movements at four different locations: a fourparameter scenario versus a seven-parameter scenario.

24 a 0 0 0 1* 0 *o.F C.0 10 to - o o X -26 -28 -30. -32 -34 0 5 10 15 20 time frame -22.. _ -24 -26 -28 -30 -32 -34. 0 5 10 15 20 time frame -9 -11 -' C -12 -1: 0 5 10 15 20 time frame 90 8s 80 7B 0 5 10 15 22 time frame De a w.0 a 5 '0 C.0 a 00._ C A 4,,, ---- - _ 2 0 '., ",,_ ' *'. ~. _ -2 -4 0 5 10 15 20 time frame 100 =. lot.._ 90 80 70 60 50. 0 5 10 15 20 time frame 34 32 2 30 ~. 28 26 24 22. 0 5 10 15 20 3 C- _- _ _ _ _ _ —_ I 0 L. time frame Figure 9. Predicted (solid) versus actual (dotted) profiles of seven angles for an up-down motion trial with C 1=1.43 deg, C2=0.32 deg/frame, and C3=1.79 deg.

25 10 E 0 *D 0 o3 E 0 c 13 S1S2 S3 S4 S5 S6 S1 S2 S3 S4 S5 S6 S1 S2 S3 S4 S5 S6 S1 S2 S3 S4 S5 S6 AP13 AP24 AP57 AP68 Oi forearm * upper arm D clavicle * torso Mean O forearm * upper arm O clavicle * torso Mean O forearm * upper arm. O clavicle B torso S1 S2 S3 S4 S5 S6 S1 S2 S3 S4 S5 S6 S1 S2 S3 S4 S5 S6 S1 S2 S3 S4 S5 S6 ML12 ML34 ML56 ML78 0 - L 0) co.o 0. 0).m c cm 0 N a r 20 10 S1 S2 S3 S4 S5 S6 S1 S2 S3 S4 S5 S6 S1 S2 S3 S4 S5 S6 S1 S2 S3 S4 S5 S6 Mean UD51 UD62 UD73 UD84 Figure 10. Normalized weighting parameters assigned to the four segments for different types of reaching movements: AP (upper panel), ML (middle panel), and UD (lower panel).

26 Discussions Though optimization-based approaches to model human dynamic postures or movements have been widely appreciated, and attempted with a certain degree of success for a variety of biomechanical systems, the implementation of three-dimensional dynamic models with a desirable level of sophistication remains to be a challenge. Previous modeling approaches concerning movements can be generally classified as sequential static (Ryan, 1970; Kilpatrick, 1970), or dynamic (Ayoub et al., 1974; Ayoub and Hsiang, 1992; Yamaguchi, 1990). Both approaches postulated an optimal strategy explicitly formulated as an objective function. Methodologically, the former involved intricate optimization procedures, one for each time frame, whereas the latter usually resorted to dynamic programming or optimal control which suffers from the problem of dimensionality. It has been acknowledged in several previous efforts that high complexity and extreme computational demand are the primary obstacle to the quest for models that are three-dimensional, dynamic, and possess a sufficient number of degrees of freedom (Hsiang, 1992; Girard, 1991; Yamaguchi, 1990). The approach presented in this study offers a potential means of circumventing the above difficulties. Rather than working in the displacement domain or coping with nonlinear musculoskeletal dynamics (equivalent to working in acceleration domain), this approach resolves the inverse kinematics problem in the velocity domain to take advantage of the linearity (see equation 2). In our modeling approach, a set of tunable weighting parameters quantifies an effort-allocation strategy as an implicit objective function in optimal motion control. Once these parameters are given, joint angle profiles can be delivered via forward integrations (see equation 10, which consumes about 20 second CPU time to compute for a 25 frame motion using Mathematica on a 7100/66 Power Macintosh). This process is computationally much more efficient than, given a known objective function, to determine motion profiles using previous approaches aforementioned. For instance, Yamaguchi et al. (1995) reported a 3.5 hour CPU time (on a Silicon Graphics

27 IRIS) consumed to solve the kinematic redundancy of a 5-DOF biomechanical system in a 3 second motion, using a supposedly very efficient algorithm based on dynamic programming. It appears that a recent study by Jung et al. (1995) also recognized the availability of the velocity domain (differential) inverse kinematics method for posture modeling. In that study, while it is perplexing how angular profiles were generated (without specifying hand motion trajectory), the posture determination was clearly accomplished by iterative resolutions of an optimization problem minimizing the sum of deviations of joint angles from their neutral position. Therefore, at least in the respect of optimization procedures employed, their model still belongs to the sequential static category. Clearly, the success of our modeling scheme relies on the ability to determine the set of parameters such that resulting joint angle profiles well represent the observed ones. The present study introduced a comprehensive set of measures for evaluating modeling performance, devised an optimization procedure for determining the parameters (as variables), and identified good (though approximate) solutions using simulated annealing. The empirical evaluation using measured reaching motion data demonstrated a promising modeling performance by the proposed approach. Although there is no data available in the literature for evaluating what is a "good" or "bad" dynamic posture prediction model, the overall 3.3 degree mean error of predicted angles with respect to the measured angles (referred to Table 2) is intuitively satisfactory. Particularly, a superb modeling accuracy was achieved for medial-lateral (ML) and up-down (UD) movements with which a better inter-joint coordination appeared to be associated. A comprehensive accuracy comparison of the current modeling approach, applied to larger scale motion data, with a number of previous approaches (e.g., Ryan, 1970; Kilpatrick, 1970) is underway. The comparison will be conducted based on similar target-driven seated reaching movements. The general success of using a segmental strategy with four time-invariant ratecontrol parameters in modeling coordinated multi-segment movements leads to several

28 important theoretical interpretations. First, the fact that the parameters (and so as the inherent strategy they represent) are time-invariant lends support to the motor program theory of movement control. The generalized motor program notion asserts that the production of (at least some) movements has invariant characteristics within the structure of movement pattern (Schmidt, 1975; Young and Schmidt, 1991). Several lines of evidence that are supportive of the motor program notion came from animal deafferentation research, examination of reaction time behavior and EMG patterns during (mostly simple arm) movements (see Young and Schmidt, 1991 for an excellent review). The current study provides a new line of evidence, in a unique way, by extracting a set of invariant metrics that characterize common seated reaching movements. Although a direct link between the "strategy" in this work and the motor program seems groundless, the identified additional invariance presents a straightforward form of visualizing an abstract concept. Secondly, a segmental strategy (i.e., four-parameter scenario) claims that the energy-type effort is allocated among body segments rather than individual degrees of freedom. This strategy implies not only an inter-segment coordination (as reflected by fixed weights assigned to each segment) but also a within-segment "synergy" (as reflected by the same weight shared by two degrees of freedom within a segment). This work, for the first time, distinguishes these two levels of coordination-whether a segmental strategy can accurately model a movement indicates what level of coordination is achieved. Indeed, the modeling accuracy (referred to Figures 4, 6, 8) differentiates those more coordinated, less difficult movements from those less coordinated, more difficult ones (e.g., AP motions vs. ML or UD motions; 0-degree offset angle vs. 45-degree offset angle, close-in vs. far-out, etc.). The legitimacy of proposing above two levels of coordination is further justified by the result that a seven-parameter scenario improved substantially the modeling accuracy for AP movements (which a four-parameter scenario did not model well) but barely for the ML or UD movements. One comment that should be made here in this context is that the distinction between the two levels of coordination cannot be revealed in a

29 two-dimensional modeling effort, since each segment would possess only one degree of freedom. Thirdly, the nature of being a velocity domain approach and its considerable computational advantage over either the displacement or acceleration domain approach also has interesting implications. It may suggest a primary role of the first-order or rate control (as opposed to the O-order or second-order control) in determining dynamic postures. Using an analogy of CPU as the neural computing unit used by central nervous system (CNS) for motor control, we may be able to reason that optimization-based O-order or second-order control is simply too slow for solving kinematic redundancy and realizing a movement. The necessity of figuring out optimally positioned postures from moment to moment (as in sequential static models), or selecting one from an infinitive number of possible variations (as in models based on dynamic programming), is questionable. An alternative which facilitates much quicker information processing and decision making, as offered by the rate control mechanism, is perhaps more likely to be the one adopted by human beings. An additional appealing aspect of the proposed modeling scheme is that it can accommodate problems with a larger number of degrees of freedom. This is an especially important attribute, considering the trend towards increasingly complex, large scale biomechanical models (Yamaguchi, 1995). Seven degrees of freedom are possessed by the current model describing the upper torso and right arm, which already compares favorably with most existing relevant models. If more degrees of freedom were to be incorporated, the Jacobian matrix would become more sophisticated but yet derivable; in the search for weighting parameters using simulated annealing, the number of operations remains the same even though dimensions of involved matrices and vectors would be increased. Therefore, the complexity of applying current approach would conceivably increase by introducing more degrees of freedom, but not substantially, and in no way exponentially as those previous approaches relying on dynamic programming (Yamaguchi, 1990).

30 It would be a remiss not discussing the limitations of the proposed approach. One limitation of the model arises from its dependence on a specification of hand (endeffector) motion trajectory. This specification may in effect reduce the applicability of this modeling approach. An improvement can be made through the implement of a separate model which predicts the hand trajectory, provided the initial and terminal hand positions. For point-to-point extracorporeal simple reaching motions, the hand trajectory has been well researched. One of the most robust results reported by several experimental studies (Morasso, 1981, 1983; Flash and Hogan, 1985) is that the hand trajectory is essentially straight with a bell-shaped velocity profile. A model derived from a minimum jerk theory for predicting such hand motion trajectory is available and can be readily utilized (Flash and Hogan, 1985). Nevertheless, for more complex hand motions, the issue of characterizing their trajectories still remains a contentious research topic. Another limitation is that the best level of modeling accuracy graded for the less coordinated, more arduous movements, particularly the anterior-posterior (AP) movements. Though the accuracy was considerably improved using a modeling scenario with seven parameters, yet it was not as satisfying what was achieved for other types of motions even with a four-parameter scenario (see Figures 4, 6, 8). This limitation leads to two possible avenues for further research. One is to implement a modeling scheme with time-variant weight parameters, starting with piecewise time-invariant parameters. This generalization of the current approach would certainly entail a great deal of added intricacies in modeling. But it may prove to be worthwhile, if it helps us gain a more comprehesive insight by detecting possible strategic change(s) during the course of a variety of movements. The other investigative line worth exploring is, as alluded earlier in the discussion, to take advantage of the modeling performance difference for movement classification purpose. A tool that stems from the current modeling approach may be developed as a gauge to differentiate more coordinated vs. less coordinated, or normal vs. pathological reaching motion patterns.

31 References Andrews, J.G. (1995). Euler's and lagrange's equations for linked rigid-body models of three-dimensional human motion. In Allard, P. Stokes I.A.F., and Blanchi, J. (Eds.), Three-dimensional analysis of human movement (pp. 145-175). Champaign, IL: Human Kinetics. Ayoub, M.A., Ayoub, M.M., and Walvekar, A.G. (1974). A biomechanical model for the upper extremity using optimization techniques. Human Factors. 16. 585-594. Ayoub, M.M., and Hsiang, M.S. (1992). Biomechanical simulation of a lifting task. In Kumar, S. (Ed.), Advances in industrial ergonomics and safety IV (pp. 831-838). London; New York: Taylor & Francis. Bean, J.C., Chaffin, D.B., and Schultz, A.B. (1988). Biomechanical model calculation of muscle contraction force: a double linear programming.method. Journal of Biomechanics. 21 59-66. Byun, S.N. (1991). A computer simulation using a multivariate biomechanical posture prediction model for manual material handling tasks. Unpublished doctoral dissertation, University of Michigan, Ann Arbor, MI. Chao, E.Y., and An, K.N. (1978). Graphical interpretation of the solution to the redundant problem in biomechanics. Journal of Biomechanics. 21, 249-242. Crowninshield, R.D. and Brand, R.A. (1981). A physiologically based criterion of muscle force prediction in locomotion. Journal of Biomechanics. 14. 793-801. Denavit, J. and Hartenberg, R.S. (1955). A kinematic notation for lower pair mechanisms based on matrices. ASME Journal of Applied Mechanics. 22. 215 -221. Eglese, R.W. (1990). Simulated annealing: a tool for operation research. European Journal of Operation Research. 46, 271-281. Flash, T., and Hogan, N. (1985). The coordination of arm movement: an experimentally confirmed mathematical model. Journal of Neuroscience. 7. 1688-1703.

32 Girard, M. (1991). Constrained optimization of artificial animal movement in computer animation. In Badler, N.I., Barsky, B.A., and Zeltzer, D. (Eds.), Making them move: mechanics, control, and animation of articulated figures (pp. 209-232). San Mateo, CA: Morgan Kaufmann Publisher, Inc. Hardt, D.E. (1978). Determining muscle forces in the leg during normal human walkingan application and evaluation of optimization methods. Journal of Biomechanical Engineering, 100, 72-78. Hsiang, M.S. (1992). Simulation of manual material handling. Unpublished doctoral dissertation, Texas Tech University, Lubbock, TX. Hughes, R.E. (1991). Empirical validation of optimization-based lumbar muscle force prediction models. Unpublished doctoral dissertation, University of Michigan, Ann Arbor, MI. Jung, E.S., Choe, J., and Kim S.H. (1994). Psychophysical cost function of joint movement for arm reach posture prediction. In Proceedings of the Human Factors and Ergonomics Society 38th Annual Meeting (pp. 636-640). Santa Monica, CA: Human Factors and Ergonomics Society. Jung, E.S., Kee, D., and Chung, M.K. (1995). Upper body reach posture prediction for ergonomic evaluation. International Journal of Industrial Ergonomics, 16, 95-107. Kirkpatrick, S., Gelatt, Jr., C.D., and Vecchi, M.P. (1983). Optimization by simulated annealing. Science. 220, 671-680. Kilpatrick, K.E. (1970). Models for the design of manual work station. Unpublished doctoral dissertation, University of Michigan, Ann Arbor, MI. Morasso, P. (1981). Spatial control of arm movements. Experimental Brain Research, 42, 223-227. Morasso, P. (1983). Three dimensional arm trajectory. Biological Cybernetics. 48, 187 -194.

33 Nakamura, Y. (1991). Advanced robotics: redundancy and optimization. Reading, MA: Addison-Wesley. Pandy, M.G., Anderson, F.C., and Hull, D.G. (1992). A parameter optimization approach for the optimization control of large scale musculoskeletal systems. Journal of Biomechanical Engineering. 114, 450-460. Park, K.S. (1973). Computerized simulation model of posture during manual material handling. Unpublished doctoral dissertation, University of Michigan, Ann Arbor, MI. Ryan, P.W. (1970). Cockpit geometry evaluation (Joint Army-Navy Aircraft Instrumentation Research Report 700201). Seatle, WA: The Boeing Company. Schmidt, R.A. (1975). A schema theory of discrete motor skill learning. Psychological Review, 82, 225-260. Seireg, A., and Arvikar, R.J. (1975). The prediction of muscle load sharing and joint forces in the lower extremity during walking. Journal of Biomechanics. 8, 89-102. Whitney, D.E. (1969). Resolved motion rate control of manipulators and human prostheses. IEEE Transactions on Man-Machine Systems. 10, 47-53. Young, D.E., and Schmidt, R.A. (1991). Motor programs as units of movement control. In Badler, N.I., Barsky, B.A., and Zeltzer, D. (Eds.), Making them move: mechanics, control, and animation of articulated figures (pp. 129-155). San Mateo, CA: Morgan Kaufmann Publisher, Inc. Yamaguchi, G.T. (1990). Performing whole-body simulation of gait with 3-D, dynamic musculoskeletal models. In Winters, J.M., and Woo, S.L-Y. (Eds.), Multiple muscle systems: biomechanics and movement organization (pp. 663-679). New York: Springer-Verlag. Yamaguchi, G.T., Moran, D.W., and Si, J. (1995). A computationally efficient method for solving the redundant problem in biomechanics. Journal of Biomechanics. 8, 999-1005.

34 Zajac, F.E., and Gordon, M.E. (1989). Determining muscle's force and action in multiarticular movement. Exercise and Sport Science Review. 17. 187-230. Zhang, X., and Chaffin, D.B. (1996). Task effects on three-dimensional dynamic postures during seated reaching movements: an analysis method and illustration. In Proceedings of the Human Factors and Ergonomics Society 40th Annual Meeting (pp. 594-598). Santa Monica, CA: Human Factors and Ergonomics Society.

35 Appendix A: Jacobian Derivation Notational Conventions A vector is expressed with a leading superscript donating the coordinate system or frame to which it is referenced, unless it is referenced to the general frame (frame [0]). Therefore, 1P represents a vector referenced to frame [1]; either ~P or P denotes a vector referenced to the general frame. A 4x4 homogeneous transformation matrix jT achieves the mapping of a vector from its description in frame [i] to its description in another frame [j] by 'P=-TJ P. This transformation requires a "1" to be added as the last component of the vectors. Transformations The transformation matrices between two neighboring coordinate systems are: cos 02 sin 01 sin 02 -cos 01 sin 02 0 COS 03 sin 03 2T= 0 0 -sin 03 cos 03 0 0 0 cos 01 sin 01 0 0 i: 0 - 1 0 0 cos04 sin 04 0 sin 02 -sin 01 cos 02 cos 01 cos 02 0 11 sin 02 l1 COS 01 COS 02 1 2 sin a cos 03 12 sin a sin 0 - 12 cos a 1 - sin 05 2T = sin 04 cos 05 -cos 04 cos 05 0 cos 05 sin 04 sin 05 -cos 04 sin 05 0 -13 sin 05 13 sin 04 cos 05 1 cos 06 cos 07 3 cosS 06 sin 07 43 -sin 06 0 - sin 07 cos 07 0 0 sin 06 cos 07 sin 06 sin 07 cos 06 0 l4 cos 06 COS 07 14 cos 06 sin 07 -4 sin 06 1

36 The transformation from frame [0] to frame [4] is then obtained by concatenating the above matrices as 4-T=T2T23T34T. Jacobian The hand position with respect to the general reference frame can be expressed as r=0r=4~ T.4r=4T [O 0 0 1 The first three elements of r as the three Cartesian coordinates with respect to the general frame are: x = 14 cosS6 sin 07(-cos 02 sin 03 cos 04 + sin 02 sin 04) + 14 cos 06 cos 07 (- cos 02 sin 03 sin 04 cos 05 - sin 02 cos 04 cos 05 - cos 02 cos 03 sin 05) +12 sin a cos 02 cos 03 - 3 cos 02 sin 03 sin 04 cos 05 + 12cos a sin02 + 1 sin 02 - l3sin 02 cos 04 cos 05 -13 cos 02cos 03sin 05 - 4 (cos 02 cos 03 cos 05 - cos 02 sin 03 sin 04 sin 05 - sin 02 cos 04 sin 05) sin 06 y = 14 cos 06 sin 07 (- sin 01 cos 02 sin 04 + cos 04 (cos 01 cos 03 - sin 01 sin 02 sin 03)) + 14 cos 06 cos 07(sin 01 cos 02 cos 04 cos 05 + sin 04 cos 05(cos 01 cos 03 - sin 01 sin 02 sin 03) - sin 05(cos 01 sin 03 + sin 01 sin 02 cos 03))12 cos a sin 01 cos 02 - 11 sin 01 cos 02 + 13 sin 01 cos 02 cos 04 cos 05 + 12 sin a cos 01 sin 03 + 12 sin a sin 01 sin 02 cos 03 + 13 sin 04 cos 05 (cos 01 cos 03 - sin 01 sin 02 sin 03) - 13 (cos 01 sin 03 + cos 01 sin 02 cos 03 sin 05) - 14(cos 05(cos 01 sin 03 + sin 01 sin 02 cos 03) + sin 01 cos 02 cos 04 sin 05 + sin 04(cos 01 cos 03 - sin 01 sin 02 sin 03)sin 05)sin 06 z = 14 sin 06 cos 07 (cos 01 cos 02 sin 04 + cos 04(sin 01 cos 03 + cos 01 sin 02 sin 03)) + 14 cos 06 07(- cos 01 cos 02 cos 04 cos 05 +sin 04 cos 05 (sin 01 cos 03 + cos 01 sin 02 sin 03) - (sin 01 sin 03 - cos 01 sin 02 cos 03 )sin 05) + 12 cos acos 01 cos 02 + 11 cos 01 cos 02 - 13 cos 01 cos 02 cos 04 cos 05 + 12 sin a sin 01 sin 03 - 12 sin a cos 01 sin 02 cos 03 + 13 sin 04 cos 05 (sin 01 cos 03 + cos 01 sin 02 sin 03) -/3 (sin 01 sin 03 - cos 01 sin 02 cos 03 sin 05) - 14(cos 05 (sin 01 sin 03 - cos 01 sin 02 cos 03) - cos 01 cos 02 cos 04 cos 05 + sin 04 (sin 01 cos 03 + cos 01 sin 02 sin 03)sin 05 )sin 06

37 The Jacobian matrix is derived as dx dx dx dx dx dx dx dO9 dO2 dO3 dO4 do5 d06 d07 =d=r y y dy dy dy dy do01 d2 dO3 d04 d05 d06 d07 dz dz dz dz dz dz dz d 0 do2 d03 do4 dO05 dQ6 d07

38 Appendix B: Simulated Annealing Algorithm Arguments To: starting temperature; Xo: initial guess of solution hf: heating factor; cf: cooling factor; fnb: neighborhood function; feal: evaluation or objective function; Texit: temperature to exit the search; Vexit: feval value to exit the search; Nexit: maximum number of iterations allowed. Variables nacpt: number of accepted solutions; nwacpt: number of accepted worse solutions; Xbest: best solution so far achieved; Vbest: best feval value so far achieved; T: current temperature; X: current solution; V: current value; i: number of iterations completed; converged: Boolean variable; heating: Boolean variable. Procedures (for maximization) Step 1. Initialization: T = To; X = X0; V = feval(Xo); Xbest = Xo; Vbest = V; i = 0; nacpt = 0; nwacpt = 0; converged = false; heating = true if hf>l, false otherwise. Step 2. Set i = i+l; If i = Nexit, go to Step 6. Step 3 Generate a candidate solution: Xcand = fnb (X); Evaluate Vcand = feval (Xcand); Compute the probability of accepting the candidate solution: P = Min[l, e(Vca-V)/T]. Step 4. Generate a random number between 0 and 1; If the number < P, nacpt = nacpt+1, X = Xcand, V = Vcand;

39 and if Vcand < X, nwacpt = nwacpt+1; If Vcand > Vbest, Xbest = Xcand, Vbest = Vcand. Step 5. If nwacpt > 2, heating = false. If heating = true, T = T*hf; otherwise, T = T*cf If T < Texit or Vbest > Vexit converged = true; If converged = false, go to Step 2. Step 6 Output Vbest, Xbest, i, T, converged, etc.