|Home | About | Journals | Submit | Contact Us | Français|
Global optimization algorithms (e.g., simulated annealing, genetic, and particle swarm) have been gaining popularity in biomechanics research, in part due to advances in parallel computing. To date, such algorithms have only been applied to small- or medium-scale optimization problems (< 100 design variables). This study evaluates the applicability of a parallel particle swarm global optimization algorithm to large-scale human movement problems. The evaluation was performed using two large-scale (660 design variables) optimization problems that utilized a dynamic, 27 degree-of-freedom, full-body gait model to predict new gait motions from a nominal gait motion. Both cost functions minimized a quantity that reduced the knee adduction torque. The first one minimized foot path errors corresponding to an increased toe out angle of 15 deg, while the second one minimized the knee adduction torque directly without changing the foot path. Constraints on allowable changes in trunk orientation, joint angles, joint torques, centers of pressure, and ground reactions were handled using a penalty method. For both problems, a single run with a gradient-based nonlinear least squares algorithm found a significantly better solution than did 10 runs with the global particle swarm algorithm. Due to the penalty terms, the physically-realistic gradient-based solutions were located within a narrow “channel” in design space that was difficult to enter without gradient information. Researchers should exercise caution when extrapolating the performance of parallel global optimizers to human movement problems with hundreds of design variables, especially when penalty terms are included in the cost function.
Optimization methods are frequently used to develop human movement simulations that reproduce experimentally measured motions or predict new motions for which experimental data are not available [1–19]. These simulations typically utilize complex musculoskeletal computer models possessing multiple degrees of freedom (DOFs). Since 10 or more design variables are often used to parameterize the control of each DOF [1,3,17,20], the resulting movement optimization problems can possess hundreds of design variables, increasing the likelihood of encountering local minima. These large-scale problems can also have a high computational cost due to iterative evaluation of the cost function and constraints. Though ineffective against local minima, gradient-based optimization algorithms have commonly been used to solve large-scale problems, primarily because of their rapid convergence properties. However, even with these algorithms, the complexities of present-day biomechanical models can necessitate thousands of function evaluations to achieve convergence [9,21].
To increase throughput, researchers have developed and evaluated several parallel algorithms for biomechanical optimization problems. Gradient-based optimization algorithms were the first to be parallelized for human movement problems [4,5,22], and more recently, global optimization algorithms such as genetic , simulated annealing , and particle swarm [23–26] have been parallelized as well. Global algorithms are advantageous since they decrease the risk of being trapped in a local minimum, but the price is significantly higher computational cost compared to a gradient-based algorithm. To date, only gradient-based algorithms have been used for large-scale (> 100 design variables) human movement problems (Table 1, top half) such as the three-dimensional jumping and gait problems with over 800 design variables solved by Anderson and Pandy [4,5]. Since global algorithms have been used successfully to solve small-to medium-scale human movement optimization problems (Table 1, bottom half), it is tempting to assume that parallel versions of these algorithms will work well on large-scale human movement problems. However, the reasonableness of this assumption has not been evaluated.
The objective of this study was to evaluate the ability of a parallel particle swarm global optimization algorithm to solve large-scale human movement problems. Algorithm performance was compared to that of a commonly used gradient-based optimization algorithm. Two large-scale human movement problems possessing 660 design variables were developed to perform the evaluation. Both problems involved predicting novel gait motions using inverse dynamic optimization of a three-dimensional, 27 degree-of-freedom (DOF), full-body gait model calibrated to match gait data collected from a patient with knee osteoarthritis (OA). The optimizations predicted how gait modifications would reduce the peak knee adduction torque subject to several “reality constraints” on the predicted motion. The results of the study provide insight into the limitations of global optimization algorithms for solving large-scale human movement problems.
A parallel particle swarm algorithm was chosen as the global optimization algorithm for attempting two large-scale human movement problems. This algorithm has been tested extensively and found to perform well on a number of small- and medium-scale benchmark problems . For comparison purposes, both large-scale human movement problems were also solved using the Levenberg-Marquardt nonlinear least squares algorithm (NLS)  available in Matlab (The Mathworks, Natick, MA). This gradient-based algorithm was selected since both optimization problems were formulated as tracking problems that minimized sums of squares of errors, the situation for which nonlinear least squares algorithms are designed. In general terms, gradient-based algorithms calculate a search direction in multi-dimensional design space, step in that direction until a minimum is found, and then iterate the entire process until no further improvement can be made. The following sections describe the parallel particle swarm algorithm and the two large-scale human movement problems used for evaluation.
Particle swarm global optimization (PSO) is a class of derivative-free, population-based computational methods introduced by Kennedy and Eberhart in 1995 . In the original PSO algorithm, particles (design points) are distributed throughout the design space and their positions and velocities are modified based on knowledge of the best solution found thus far by each particle in the “swarm.” Attraction toward the best-found solution occurs stochastically and uses dynamically adjusted particle velocities. Particle pseudo-velocities (Eq. 1) and positions (Eq. 2) are updated as shown below:
where represents the current position of particle i in design space and subscript k indicates a (unit) pseudo-time increment. The point is the best-found position of particle i up to time step k and represents the cognitive contribution to the search velocity . The point is the global best-found position among all particles in the swarm up to time step k and forms the social contribution to the velocity vector. Random numbers r1 and r2 are uniformly distributed in the interval [0, 1]. Parameters c1 and c2 are the cognitive and social scaling parameters, respectively, whose values are generally chosen such that their sum equals 4 . Choosing c1 = c2 = 2 causes the particles to overshoot the attraction points and half the time, thereby maintaining separation in the group while allowing a greater area to be searched than if the particles did not overshoot. The value of the variable w, set to 1 at initialization, is reduced dynamically based on the cost function improvement rate, thereby reducing the search area gradually (see  for further details).
Two parallel PSO algorithms have been reported in the literature. One algorithm is called parallel synchronous PSO (PSPSO) , where the master processor updates particle positions and velocities only after all function evaluations have been completed by all slave processors. The other algorithm is called parallel asynchronous PSO (PAPSO) [26,29], where the master processor updates particle positions and velocities continuously as function evaluations are completed by the slave processors. Since PAPSO tends to find a slightly lower cost function value than does PSPSO for the same number of function evaluations , a parallel PAPSO algorithm was used in the present study.
Ten PAPSO runs were performed for both large-scale movement optimization problems. Each run was performed on a homogeneous Linux cluster of 20 identical PCs, utilized 20 particles along with other standard algorithm parameter values (see  for details), and iterated for 1 million function evaluations to ensure convergence. One of the particles was always assigned to match the nominal experimental gait motion, while the other 19 particles were given random initial guesses within the initial bounds of the design space. During each optimization, particles were allowed to “stray” beyond the initial bounds. Similar to previous optimization algorithm evaluations studies [13,24], no tuning of algorithm parameter values was performed for any of the optimizations. Execution time for each PAPSO optimization was approximately 2.3 hours.
Previously reported experimental gait and isolated joint motion data collected from a single patient with knee OA (male, age 37 years, height 170 cm, mass 69 kg, alignment 5° varus) were used as the initial guess for the two benchmark human movement optimization problems . Institutional review board approval and informed consent were obtained. The gait and isolated joint motion data were used to calibrate the joint and inertial parameter values in a dynamic, full-body, three-dimensional (3D) gait model. The model possessed 27 degrees-of-freedom (DOFs) composed of gimbal (3 DOFs), universal (2 DOFs), and pin (1 DOF) joints, where each DOF was chosen as a generalized coordinate in the model (Fig. 1). The position and orientation of the pelvis in the laboratory coordinate system were defined by three translational and three rotational DOFs, resulting in six pelvis residual loads (three forces and torques) during inverse dynamic analyses. For the lower body, each hip was modeled as a gimbal joint, each knee as a pin joint, and each ankle as two non-intersecting and non-orthogonal pin joints. For the upper body, the back was modeled as a gimbal joint, each shoulder as a universal joint, and each elbow as a pin joint. Ground reaction forces and torques applied to the feet were treated as unknowns to be determined during periods when the foot was known to be in contact with the ground . Details on the experimental data collection and model calibration process can be found in .
The calibrated model and nominal gait data were used in two large-scale human movement optimization problems. Both problems involved predicting gait changes that reduce the left knee adduction torque subject to reality constraints that produced realistic-looking gait motions. The first problem predicted how increasing the toe out angle on each side by 15 deg would alter the left knee adduction torque (henceforth called the “toe out gait optimization”) , while the second problem predicted changes in gait motion that would decrease the left knee adduction torque without altering the toe out angle (henceforth called the “modified gait optimization”). The reality constraints minimized changes in trunk rotations, pelvis residual loads, and foot paths (with increased toe out angle where necessary) away from the nominal gait data.
The two gait prediction problems were formulated as inverse dynamic optimizations. Each optimization performed repeated inverse dynamics analyses, where each analysis used the patient-specific gait model to calculate joint torque and pelvis residual load outputs given the current guess for the joint angle and ground reaction inputs. The joint angle and ground reaction curves measured experimentally were parameterized as a function of time using a combination of polynomial and Fourier terms , thereby allowing analytic calculation of joint velocities and accelerations. To match the experimental gait data accurately, a cubic polynomial with eight Fourier harmonics (i.e., 20 coefficients) was selected, producing root-mean-square (RMS) errors between experimental and fitted curves on the order of 0.3 mm, 0.2°, 4 N, and 0.7 Nm . The design variables for the two optimization problems were defined as the coefficients of the parameterized curves. By changing the joint angle and ground reaction design variables, the optimizations were able to predict novel gait motions that satisfied the dynamics equations [6–8,30]. For both problems, the following motion and ground reaction curves were allowed to vary: all pelvis translations, all back, pelvis, hip, knee, and ankle rotations, and all ground reaction forces and torques (33 curves in all for a total of 660 design variables). Shoulder and elbow rotations were prescribed to match the experimental data. Ground reaction forces and torques were set to zero for time frames when the foot was known to be off the floor .
Since the PAPSO and NLS algorithms can only solve unconstrained optimization problems, the constraints for both movement prediction problems were handled using a penalty method. The form of the augmented cost function is indicated in Eq. (3) below:
In this equation, nframes is the number of time frames used in the problem (= 101), f is a specific time frame, j is a translational or rotational joint axis number, and s is a side. The weights wTrunk, wPelvis, and wFootPath (= 10) for the reality constraints were initialized to 1 and increased by factors of 10 until RMS errors in trunk orientation, pelvis residual loads, and foot path were less than 5 mm, deg, N, or Nm. TLAdd is the left knee adduction torque, which was removed for the toe out gait optimization and minimized for the modified gait optimization. Other constraints were changes in the following quantities away from their nominal experimental trajectories: hip flexion/extension, abduction/adduction, and inter/external rotation torque and angle (ΔTHip, ΔqHip), knee flexion/extension torque and angle (ΔTKnee, ΔqKnee), ankle flexion/extension and inversion/eversion torque and angle (ΔTAnkle, ΔqAnkle), anterior-posterior and medial-lateral center of pressure position (ΔCoP), pelvis translations and rotations (ΔqPelvis), trunk rotations (ΔqTrunk), pelvis residual forces and torques (ΔTPelvis), and foot translations and rotations (ΔqFoot). Though foot translations and rotations and trunk rotations with respect to the laboratory coordinate system were not generalized coordinates in the model, these quantities were calculated from the current values of the generalized coordinates.
In contrast to the NLS algorithm, the PAPSO algorithm converged to an improved solution only for the toe out gait optimization problem. For both problems, the best final cost function value found by 10 PAPSO runs was larger than the final cost function value found by one NLS run (Table 2). With the exception of foot path rotation error for the PAPSO toe out gait optimization, all optimizations successfully satisfied the three reality constraints (trunk orientations, pelvis residual loads, and foot path translation/rotation,), maintaining RMS errors to within 2 mm, 2 deg, 4 N, and 4 Nm.
For the toe out gait optimization, the PAPSO algorithm predicted little change in the first knee adduction torque peak and a 9.3% reduction in the second peak, while the NLS algorithm predicted little change in the first peak and a 32.2% reduction in the second one (Fig. 2). The PAPSO algorithm achieved only a 10 deg increase in toe out angle whereas the NLS algorithm achieved the full 15 deg increase. The corresponding lateral shift in the CoP was larger for the NLS algorithm than for the PAPSO algorithm (Fig. 3). To achieve an increased foot progression angle, both algorithms predicted increased hip abduction, greatly increased hip external rotation, slightly decreased ankle flexion, and increased ankle eversion during stance phase.
For the modified gait optimization, PAPSO was unable to produce any adduction torque reductions, while NLS predicted reductions of 32.4% and 29.4% in the first and second peak, respectively (Fig. 4). These reductions were accompanied by a small lateral shift in the CoP (Fig. 5). The predicted kinematic changes were decreased pelvic obliquity, slightly increased hip, knee, and ankle flexion, and increased hip internal rotation, resulting in a medial shift of the left knee center.
This study evaluated whether a global particle swarm optimization algorithm could solve two large-scale human movement problems possessing 660 design variables. Both problems utilized a three-dimensional, 27 DOF, full-body gait model whose joint and inertial parameters were calibrated to gait data collected from a single patient with knee OA. The cost function for both problems sought to reduce the peak knee adduction torque either by increasing the foot progression angle (i.e., the toe out gait optimization) or by minimizing the adduction torque directly without any changes in foot path (i.e., the modified gait optimization). Penalty terms were included to ensure that the predicted gait motions were physically realistic. The PAPSO algorithm achieved an improved solution only for the toe out gait optimization problem, whereas a gradient-based NLS algorithm achieved significant and realistic improvements for both problems. These findings suggest that global optimization algorithms proven to work well for small- to medium-scale problems may not work well for large-scale human movement problems. Researchers should not assume that an optimization solution is the global minimum simply because a global algorithm was used to generate it, and they should exercise caution when extrapolating the performance of global optimizers to problems involving hundreds of design variables.
There are at least three reasons why the PAPSO algorithm performed worse than did the NLS algorithm on these two benchmark problems. First, no tuning was performed of PAPSO algorithm parameters. It is possible that our PAPSO algorithm could have performed better had time-consuming algorithm parameter tuning been performed. Such tuning was not performed in the present study since few researchers take the time to do it (i.e., researchers want to find algorithms that work well “out of the box”), and since studies that evaluate optimization algorithm performance normally do not tune the algorithm parameters to the specific problems being investigated [13,24]. Second, the shape of the design space made it difficult to locate the global minimum without the use of gradient information. Though the design space was smooth for both sample problems, the minimum was located in a small deep hole in design space (Fig. 6), making it extremely difficult for a particle to migrate into this hole using only stochastic updates. Third, only a single set of penalty weights was used for both optimizations problems rather than ramping up the weights slowly over a sequence of optimization runs . A sequence of optimizations with increasing penalty weights was not performed because of the excessive computation time required by this approach. Furthermore, the selected weights posed no problem for the NLS algorithm, suggesting that the PSO algorithm may have convergence difficulties for problems involving penalty terms. Differences in optimization algorithm performance cannot be attributed to differences in problem formulation, since both algorithms used exactly the same formulations. In fact, the PSO algorithm possesses the advantage of insensitivity to design variable scaling, which gradient-based algorithms with finite-difference gradients do not possess .
Since the global optimization results reported in this study were based on a single optimization algorithm, it is not clear that similar poor performance would be observed for other global algorithms. However, this will likely be the case based on previous benchmark problems solved with a synchronous version of our PSO algorithm and several other global algorithms . At a minimum, researchers who desire to use a global optimization algorithm on a large-scale human movement problem should compare their results with those from a gradient-based algorithm as a check.
The best result for both optimizations was consistent with recent experimental observations. The toe out gait optimization predicted that increasing the foot progression angle by 15 deg would reduce the second but not the first peak of knee adduction torque curve. When Guo et al.  increased the experimental toe out angle by 15 deg, they found little change in the first adduction torque peak and roughly a 40% reduction in the second one. When the subject in our study increased his toe out angle by approximately 15 deg, the experimental reduction in the second peak was approximately 30% with little change in the first one . By comparison, our NLS toe out gait optimization predicted a 32.2% reduction in the second peak and only a slight change in the first one, with both peaks being within ± 1 standard deviation of the mean values measured from our subject over three trials. The modified gait optimization predicted that internally rotating the hips so as to medialize the knees would reduce both peaks of the knee adduction torque curve. When Davis et al.  recently trained a healthy subject to walk with the hips internally rotated and the knees medialized, they observed a 28% reduction in the first knee adduction torque peak. Our NLS modified gait optimization predicted a comparable reduction of 32.4%.
In summary, this study evaluated the ability of a global particle swarm optimization algorithm to solve large-scale human movement problems. Though the PSO algorithm used in our study has performed well on small- to medium-scale benchmark and biomechanical optimization problems , it did not perform well on the two large-scale human movement problems investigated in this study, with a gradient-based nonlinear least squares algorithm performing better on both problems. Significant algorithm parameter tuning or use of a global-local hybrid algorithm may be necessary for PSO and other global optimizers to solve large-scale human movement problems. Based on the results of this study, the PSO algorithm is not recommended for solving large-scale human movement optimization problems possessing constraints or competing terms in the cost function. Researchers should exercise caution when evaluating the results of large-scale human movement optimization problems solved using global algorithms.
This study was supported by NIH National Library of Medicine grant R03LM007332, a Whitaker Foundation Biomedical Engineering Research Grant, and NIH National Center for Medical Rehabilitation Research grant R21HD053490 to B.J.F.