Evaluation of Direct Collocation Optimal Control Problem Formulations for Solving the Muscle Redundancy Problem

Estimation of muscle forces during motion involves solving an indeterminate problem (more unknown muscle forces than joint moment constraints), frequently via optimization methods. When the dynamics of muscle activation and contraction are modeled for consistency with muscle physiology, the resulting optimization problem is dynamic and challenging to solve. This study sought to identify a robust and computationally efficient formulation for solving these dynamic optimization problems using direct collocation optimal control methods. Four problem formulations were investigated for walking based on both a two and three dimensional model. Formulations differed in the use of either an explicit or implicit representation of contraction dynamics with either muscle length or tendon force as a state variable. The implicit representations introduced additional controls defined as the time derivatives of the states, allowing the nonlinear equations describing contraction dynamics to be imposed as algebraic path constraints, simplifying their evaluation. Problem formulation affected computational speed and robustness to the initial guess. The formulation that used explicit contraction dynamics with muscle length as a state failed to converge in most cases. In contrast, the two formulations that used implicit contraction dynamics converged to an optimal solution in all cases for all initial guesses, with tendon force as a state generally being the fastest. Future work should focus on comparing the present approach to other approaches for computing muscle forces. The present approach lacks some of the major limitations of established methods such as static optimization and computed muscle control while remaining computationally efficient. Electronic Supplementary Material The online version of this article (doi:10.1007/s10439-016-1591-9 contains supplementary material, which is available to authorized users.


Mathematical expressions for muscle-tendon characteristics
Properties of muscle and tendon were described by dimensionless characteristics. All characteristics are at least second order continuous and the tendon force-length characteristic f t is third order continuous for compatibility with the interior point numerical optimization solver that uses second order derivative information. Third order continuity of f t is required, since the first derivative of the tendon-force length curve is used in formulations 1 and 3 (see below). f t is described by an exponential function: withl T normalized tendon length, k T = 35 tendon stiffness at 4% strain, and coefficients c 1 , c 2 , and c 3 . The active force-length characteristic f act is described by the sum of three Gaussian functions: withl M normalized fiber length and coefficients b 1i , b 2i , b 3i , and b 4i for i = 1 . . . 3. The passive force-length characteristic f pas is described as in OpenSim's Thelen2003Muscle by an exponential function: The force-velocity characteristic f v is described by a logarithmic function: 2 Full-form expressions of constraints imposing muscle dynamics Contraction dynamics was imposed using four different formulations. To this aim, four different functions were derived from the Hill-model described by equations 3-7 (see manuscript). Full-form expressions for these functions are given below.

Formulation 1
In formulation 1, contraction dynamics is imposed using normalized tendon forceF f 1 can be found by taking the first time derivative of eq. 3: Evaluating eq. S6 requires normalized tendon lengthl T and normalized tendon velocitỹ v T = dl T dt .l T can be solved from eq. 3 usingF T , which is an input to function f 1 : Tendon velocity v T can be related to muscle-tendon velocity v M T and fiber velocity v M by differentiating eq. 5 with respect to time: Differentiating eq. 6 with respect to time yields: from which we can obtain an expression for dα dt : Eq. S10 can then be substituted in eq. S8, which simplifies to the following expression: Hence, we find tendon velocity from: requiring fiber velocity v M and the cosine of the pennation angle cosα. Normalized fiber velocity can be solved from eq. 4 : l M can be solved from eq. 5 and 6. These two equations can be rewritten as: Squaring both sides and adding the resulting equations together eliminates α and results in the following equation for l M : which can be evaluated using l T (eq. S7). The cosine of the pennation angle cosα can be solved from eq. 6: F M can be solved from eq. 7: Summarized, f 1 is computed by subsequently evaluating eq. S7, denormalizingl T (l T =l T l S T ), evaluating eq. S17, normalizing l M (l M = l M l 0 M ), evaluating eqs. S18 and S19, normalizing F M ), evaluating eqs. S14 and S13, denormalizingṽ ) and evaluating S6 using normalized tendon forceF T and muscle activation a as inputs. Note that muscle-tendon length l M T and muscle-tendon velocity v M T are determined by the motion and are computed prior to solving the optimization problem. l M T is computed using OpenSim's Muscle Analysis and v M T is obtained by differentiating the cubic spline interpolation of l M T in MATLAB.

Formulation 2
In formulation 2, contraction dynamics is imposed using normalized fiber lengthl M as a state: dl f 2 can be found by taking the time derivative of normalized fiber length: Normalized fiber velocity can be solved from eq. 4 (see eqs. S13 and S14). Normalized fiber lengthl M is an input to f 2 . Computing muscle force F M from eq. S19 requires F T and cosα.
To compute tendon force F T from eq. 2, tendon length l T is required. Tendon length l T can be solved from eqs. 5 and 6. Squaring both sides of eq. 6 and substituting 1 − cos 2 α for sin 2 α yields: Solving this equation for l M cosα and substituting the result in eq. 5 yields: cosα is then obtained from eq. S18.

Formulation 3
In formulation 3, contraction dynamics is imposed using normalized tendon forceF T as a state and introducing u F , the scaled time derivative of the normalized tendon force, as a new control simplifying the contraction dynamic equations: The Hill model was then imposed as a path constraint: f 4 is obtained by substituting eq. 4 in eq. 7: Evaluating eq. S26 requiresl M , cosα andṽ M in addition to inputs a andF T . Fiber length l M can be computed by subsequently evaluating eq. S7, denormalizingl T and evaluating eq. S17. The cosine of the pennation angle cosα can then be computed from eq. S18. Normalized fiber velocityṽ M can be computed from the relation between tendon, fiber, and muscle-tendon velocity (eq. S11) To obtain tendon velocity v T , eq. S6 is solved forṽ Evaluating this equation requires dF T dt andl T . dF T dt can be solved from eq. S24 and input u F andl T was computed in the process of obtainingl M . Summarized, f 3 is computed by subsequently evaluating eq. S7, denormalizingl T , evaluating eq. S17, normalizing l M , evaluating eqs. S18 and S28, denormalizingṽ T , evaluating eq. S27, normalizing v M , and evaluating eq. S26.

Formulation 4
In formulation 4, contraction dynamics is imposed using normalized muscle fiber lengthl M as a state and introducing u v , the scaled time derivative of normalized muscle length, as a new control simplifying the contraction dynamic equations: The Hill model was then imposed as a path constraint: f 4 is obtained by substituting eqs. 3 and 4 in eq. 7: can be computed from eq. S23 and l M = l 0 Ml M and cosα can subsequently be computed from eq. S18.

Algorithm settings
All non-default GPOPS-II settings are specified in Table 2.

Additional results
Computed muscle activations and tendon forces for all muscles of the complex model during walking and running are reported in Figures 1-12. Muscle activations and tendon forces computed with static optimization and muscle dynamic optimization are shown. Only the results obtained with formulation 4 are reported, since all formulations gave very similar results when they converged to a locally optimal solution. Muscle activations computed using static and dynamic optimization during walking only differ for muscles with long tendons (Figures 1-3), whereas tendon forces computed using static and dynamic optimization during walking are very similar (Figures 4-6). However, both muscle activations and tendon forces computed using static and dynamic optimization during running differ considerably for some muscles (Figures 7-12). Note that we chose to report absolute muscle forces to represent the relative importance of the different muscles.
Ideal joint torques during walking and running for the complex model are reported in Figures  13 and 14. Ideal torques during walking are small, i.e. around zero for static optimization and below 0.7Nm for dynamic optimization. Ideal torques during running are larger, i.e. up to about 15Nm, and very similar for static and dynamic optimization.