{'fraction_non_alphanumeric': 0.02753837173368128, 'fraction_numerical': 0.007657178175325344, 'mean_word_length': 4.979059741326217, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Optimization-Based Model for Full-body Reaching Movements

Background The development of a simulation model of full body reaching tasks that can predict endeffector 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 begin to explore cost functions 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 tools 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, center 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-scale human skeletal systems. Conclusions This method can reasonably predict the whole body reaching movements including final postures, joint power and movement of COM using simple algebraic calculations of inverse dynamics and forward kinematics.


Abstract Background
The development of a simulation model of full body reaching tasks that can predict endeffector 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 begin to explore cost functions 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 tools 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, center 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-scale human skeletal systems. This method can reasonably predict the whole body reaching movements including final postures, joint power and movement of COM using simple algebraic calculations of inverse dynamics and forward kinematics.

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. Various 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 endeffector 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 minimized 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 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].
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.

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 system 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 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 , or body center of mass . 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 performance measures such as total joint power, displacement of center of mass, or a combination of these two variables.
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 vastly 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.  at the end of movement if the optimization converges. This is a necessary though not sufficient condition for convergence. Actually, once the cost function is determined, the convergence depends mainly on the behavior of the optimization method.

Polynomial Type Controller
The time history of the optimal trajectory of each joint is a function of time which can be approximated by an n th -order polynomial The angular velocity and acceleration can be derived analytically from equation (3). These joint trajectories are then used as input for the inverse dynamics and forward kinematics calculations ( Fig.1). The advantage of using an inverse dynamics method is that these calculations do not require integration, and since there is a one-to-one mapping from joint space to Cartesian space in forward kinematics, the problem of kinematic redundancy is eliminated. Thus, by using a forward kinematics calculation, the problem is reduced to determining the coefficients of the polynomial. Furthermore, the whole body voluntary reaching movement can be partially described as the point where the end-effector starts for each movement of each joint. The higher the order of polynomial, the more freedom the polynomial has to approximate the joint motion, however, as the polynomial order increases the number of variables that needs to be tuned also increases.

Skeleton Dynamics
The general equations of full body motion can be written as are the vectors of joint angle, angular velocity, angular acceleration respectively; ) (θ I is segment mass inertia moment matrix; τ is vector of net joint moment; are gravity terms, Coriolis-centripetal-viscoelasticity, and external terms such as ground reaction forces. Joint viscoelasticity can be written , 0 , where K and B are joint stiffness and viscoelastic coefficient matrices respectively. One way to implement the dynamics system model is to use analytical methods. For systems with low DoF (i.e. DoF<3), the analytical expressions for the relationships between angular accelerations and joint torques can be easily written. However, a system with as few as two segments with 7 DoF, requires an analytical expression with more than 200,000 elementary operations (e.g. +,-,*, cos, sin) [2]. Although these expressions could be Once the initial values for the parameters of controllers have been input, the optimization algorithm will modify the parameters until the preset criteria or minimum parameter changes are satisfied, i.e.
where j is the index of iteration. In principle, any well developed nonlinear optimization algorithm can be used to find the optimal controller parameters. However, different algorithms will produce different results. Here we use the nonlinear Least-squares function [ ] is the gradient of criterion C respective to controller parameter p . Jacobian matrix ( J ) and the gradient of criterion ( ) are approximately calculated through parameter perturbation.

Initializing Simulation Model
As noted earlier the cost function is a composite function with a task error cost term and a physiological cost term. Simulations were conducted using the following criteria cost Power+COM). We also calculated the hand trajectory using minimal jerk criteria (min Jerk) as described by Flash and Hogan (1985) [3]. The minimum-jerk trajectory of end-effector was calculated by the following equation, when moving from location To initialize the optimization, the initial joint angles at (t=0) and (t=t f ) were set to the mean values from neutral standing posture derived from experimental data, and the controller parameters were set to zero. In this study, a 6 th order polynomial controller with two unknown variables, i.e. final posture ( f θ ) and the coefficient of polynomial ( 6 p ) is used to determine each joint excursion. Therefore, there are 72 parameters to be tuned for the 12 joints. To ensure physiological fidelity of the simulation, joint range of motion values (Eq.3) were input based on accepted norms of joint range of motion [16], and measures of joint viscoelasticity determined from the extant literature [17][18][19][20] (See Table 1).
The average height and weight of 15 healthy subjects were used in the simulations and massinertial characteristics of each segment were derived from the regression equations provided in the literature [21,22].
Final motor task error (i.e. distance between target and end-effector), parameter tolerance ( 1 ε ) and criterion tolerance ( 2 ε ) of termination of optimization algorithm were set as 2mm, 10 -6 , and 10 -6 respectively for all simulations. The constant weight 0 λ was set as 10 -8 , 10 -9 , 10 -10 (for high, middle and low targets) and 7×10 3 , 10 3 , 10 2 (for high, middle and low targets) for power and COM term respectively. For simplicity, joint weight matrix R was set as the unit matrix. The movement duration was set as 0.56s, 0.575s and 0.68s for high, middle and low targets respectively based on the experimental results from 15 healthy subjects. The solver for numerical computation was set with a fixed-step sample interval of 0.001s.

Experimental Protocols
Fifteen healthy subjects (7 males [21,22]. The experimental trajectory of end-effector including distance from target location, velocity, acceleration and jerk also were calculated from the same Simlink model.

Data Analysis
The paths/trajectories of the end-effector determined from the simulation model were compared qualitatively to the means and standard errors of path/trajectories of end-effector from the experimental data. To quantitatively assess the simulation model, for each target height, t-tests were used to compare the predicted COM displacement and predicted total joint power to the experimental data.

Optimization
The optimization for each trial only takes a few hours of CPU time on a personal computer (see Table 2), however, there was no clear effect of cost functions or target location for the CPU time required to perform the simulations.
Examination of total power squared and final COM squared (i.e total cost) in Table 2 shows that the values of these measurements are significantly reduced when adding the physiological cost terms (e.g. min Power, min COM and min Power + COM) to the Error cost term even though we used an adaptive weight coefficient to balance the physiological cost terms. Comparing all four control strategies, final COM squared was smallest when the min COM strategy was used. However, total power squared was smallest for min Power + COM control strategy. The path of the end-effector determined by the different simulation methods and the experimental data (as described above) is shown for each plane and target height in Figure 3.

End-effector Trajectories
Examination of the experimental data for the path of the end-effector in the sagittal and frontal planes ( Figure 3A & B) reveals a fairly large path curvature, which is consistent with the path predicted by the simulation using each cost function. In contrast, using minimal jerk as a cost function predicts a straight line path of the end-effector in all planes. However, examination of the experimental data for the path of end-effector in the transverse plane indicates a much straighter path. Thus, for the transverse plane, the best fit appears to be provided by either minimal jerk or the minimal task error cost function. In contrast, minimal power and COM cost functions predict a fairly large curvature and do not provide a good fit of the experimental data. Nonetheless, as stated earlier, even if the path of the end-effector can be accurately predicted, the movement strategy or the posture adopted at target contact may not be consistent with experimental data. Thus, we next compare the predicted COM displacements to the experimental data for each target height.

Displacement of COM
The trajectories of COM as predicted by the simulation models using the different cost functions are compared to the trajectories of COM from the experimental data ( Figure 4).
While, the exact path of the COM does not appear to be well fit by any of the cost functions, the change in COM from initial posture to target contact appears to be consistent with the cost function that minimized joint power and COM displacement. In fact, t-tests revealed that using a cost function that minimized COM displacement and joint power was the only cost function that was not significantly different from the experimental data (Table 3).
Minimum task error predicted the largest displacement of the COM. Specifically, for the middle and low target reaching tasks, the predicted COM displacement exceeds 20 cm and 30 cm respectively, which would clearly cause the subject to fall forward or require a step to prevent falling. Therefore minimizing end-effector error alone does not provide a reasonable solution for the movement task with respect to COM displacement. Notice that when only task error is minimized, the posture adopted at target contact does not compensate for the forward displacement of the trunk and its effect on whole body COM ( Figure 5).
Qualitatively, the posture adopted at target contact by participants (i.e. experimental data) is best fit by the cost function that minimizes COM displacement and joint power ( Figure 5).
Total Joint Power Figure 6 shows the time series of joint power ∑  Table 4). For the high target and middle target it appears that minimizing COM displacement provides the best fit to the experimental data.
However, for the low targets, none of the cost functions provide a good fit to total energy expenditure.

Skeleton Dynamics
The underlying premise of these simulations models is that an inverse dynamics approach was the preferred method for predicting movement strategies. Inverse dynamics uses as inputs the joint motions as a function of time and differentiates them twice to yield the accelerations required to calculate joint torques and interaction forces needed to produce the motions. However, the motion functions must be checked for consistency to ensure they stay within geometrical constraints of linked rigid bodies at every time frame. For inverse dynamics with closed topologies, this process is complex and the computation load is quite large. The difficulty in handling closed topologies comes from indeterminacy which a generic property of the inverse dynamics itself [13]. In the present study, the full-body reaching tasks required the subjects to stand firmly on the force plates without any foot movement during the task. Therefore, there exists a closed topology within lower extremities. To avoid the consistency problems of geometrical constraints mentioned above, a one leg model was used in the inverse dynamics calculations used in the optimization and simulation with mass properties doubled for thigh and shank segments. For the experimental data, the same techniques were applied. Given the task constraints, a one leg approximation is acceptable for full-body reaching movements and consistent with our previous work [25].

Comparison of Simulation and Experiments
The results of these simulations are in general agreement with the experimental data regarding the displacement of COM and the final postures adopted at different target locations. Perfect fits to the experimental data are not shown in this study because 1) fitting is generally associated with arbitrary parameter adjustments [2], i.e. adjusting some anthropometric parameters may help to improve the predictions [26]; 2) fitting quality may not be sufficient to estimate the validity of a model [27]. Figure 5b illustrates that minimizing Power and COM in the simulation model provided a reasonably close fit to the experimental data. However, it is also clear that there are differences between the simulation and the experimental data regarding location of the endeffector. This is particularly evident for reaches to the middle and low targets. These differences are most likely due to errors associated with motion capture (e.g. skin movement under sensors) and from model constraints (e.g. constraining shoulder translation). Shoulder translation which can dramatically affect location of the end-effector, particularly in forward and overhead reaching tasks [28,29]. The vertical and anterior-posterior translation could be up to ±3.8cm. However, even with these potential sources of error, both simulation and experimental results capture the major characteristics of the voluntary target reaching movement, i.e. curved path for end-effector, bell-shaped profile of velocity of end-effector and increasing displacement in COM and joint power with lowering the target height.

Controller
The minimum jerk solution for end-effector movement based on Euler-Poisson's theorem indicates that the optimal trajectory of end-effector must have the form of a 5 th -order polynomial [3]. Although a 6 th -order polynomial was used in our study, the value of the 6 thcoefficient is very small for minimum error trajectory (e.g.

Optimization Method
In this study, the Levenberg-Marquardt method with line search was used to find the optimal trajectory of the end-effector in terms of motor task error (end-effector error) or/and certain physiological strategies. It must be emphasized that this method belongs to the category of steepest descent. Therefore, the optimal trajectory is most sensitive to the selected criterion, i.e. the joint angles were calculated in such a way that the increments of joints correspond to the gradient (partial derivatives) of the criterion with respect to the joint angles. The gradient is the minimal change in joint angles that results in unit change of the criterion. This method is consistent with that proposed by Hinton in that each joint is moved autonomously, in proportion to how much moving that joint alone affects the end-effector-target distance. [30] In addition to using various criteria, any other optimization methods could be used and certainly different optimal trajectories for both end-effector and joint angles can be obtained. Trust region method [31], for instance, can be used to obtain the optimal trajectories for the posture comfort hypothesis if the joint's range of motion is known. A postural comfort hypothesis predicts that joint excursions in multi-joint tasks are in part determined by joint comfort [32][33][34]. Cruse (1986) pointed out that each joint has an associated discomfort function and that comfort costs influence the movement strategy chosen. The discomfort associated with an individual joint is highest near the joint's biomechanical range limits and lowest for some optimal configuration, which tends to be near the middle of the joint's range of motion [35]. Because certain joints (e.g. knee joint) are near their biomechanical range limits, particularly for reaches to the high target, it is reasonable to assume that the Trust region method would not provide reasonable results for these tasks.

Criteria
For full-body voluntary movements, at least two major performance criteria need to be accounted for, i.e. the end-effector must reach the target at the end of movement; and the whole body must maintain balance during the movement (the motion of COM is within the base of support). Other performance criteria, such as minimization of energy, may be necessary, especially for reaches to the low target. Different performance criteria may need to be adopted for different target reaching tasks or even for different periods of the reaching movements [26]. Clearly that more than one performance criterion is required to reasonably predict whole body reaching tasks. In fact, Ferry et al [26] had suggested that even for a simple arm raising task more than one performance criterion may need to be adopted.
Additionally, Parnianpour et al [6] reports that even for a movement task consisting of one segment with one degree of freedom may require more than one performance criterion. For end-effector path planning, the minimum motor task error may not be a necessary performance criteria because the location of the end-effector at each moment of time is known [29,36]. The final boundary condition in Cartesian space also can be put into optimization algorithm constraints. 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 planning, the path in Cartesian space is known; the question is how to solve a redundant inverse kinematics problem. Whereas in voluntary movement the path is unknown, there is no prior knowledge about the movement. In fact, all constraints can be combined into the movement performance criteria. For instance, the joint range can be measured with joint comfort and then combined into the criteria [37].

Conclusions
In summary, the proposed method in this paper is effective for the simulation of large-scale human skeleton systems, which can reasonable predict whole body reaching movements (i.e. final postures, movement of COM, joint power, and end-effector trajectories etc.). As applied, a combination of several control strategies such as minimizing end-effector error, joint power and COM and using the simple algebraic calculations of inverse dynamics and forward kinematics provided good fits to the experimental data. In the future different cost criteria should be examined and compared with even more complex movement tasks to further elucidate how the CNS plans and executes movements in a kinematics redundant system.

Figure 1 -Model of Optimized Controller
Block diagram for the optimal model of full body reaching movement. The outputs of controller are the joint angle functions of time. The inverse dynamics and forward kinematics models were used to calculate the physiological measurements and the location of endeffector. Optimization algorithm was used to minimize a certain criterion to produce the optimized parameters for the controller.