% ****** Start of file DC_paper.tex ******
%
%   See the REVTeX 4 README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.0
%
% See the REVTeX 4 README file
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex DC_paper.tex
%  2)  bibtex DC_paper
%  3)  latex DC_paper.tex
%  4)  latex DC_paper.tex
%
% Some ASME recommendations:
% 1) Bibtex can be used to create reference citations. However, a separate .bib file CANNOT %    be submitted. The .bbl file must be read or pasted into the .tex file before
%    submitting the final TeX file. Note: the article titles must be captured or keyed in
%    the reference section if using bibtex.
% 2) when using EPS Format for figures
%    Use for vector graphics
%    Set resolution to 600 dpi
%    If necessary, set color mode to RGB-encoded at 8 bits/channel
%
% 3) Figure Sizing
%    Graphics should be the same size as they should appear in journal.
%    Single Column Width: 8 cm (3.25 in., 20 picas)
%    Double Column Width: 16 cm (6.5 in., 40 picas)
%    Maximum Length: 23 cm (9.5 in., 57 picas)]
%
% 4) Line weights should be 0.5pt -1.5pt in thickness. WATCH: line weights below 0.5 points %    will disappear.



%
%\documentclass[twocolumn,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}
%\documentclass[prb,preprint,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}

% Suggested by ASME:
% Note: although ASME recommends using the argument "prb" for the final publication, the
%       citations seemed much better not using it, so for the draft paper we will not use
%       this argument but we may add it for the final publication.
%\documentclass[prb,twocolumn,secnumarabic,preprint,amsmath,amssymb]{revtex4}
%\documentclass[prb,preprint,secnumarabic,amsmath,amssymb]{revtex4}
%\documentclass[twocolumn,,secnumarabic,amsmath,amssymb]{revtex4}
\documentclass[preprint,showkeys,secnumarabic,amsmath,amssymb]{revtex4}

% Some other (several out of many) possibilities
%\documentclass[preprint,aps]{revtex4}
%\documentclass[preprint,aps,draft]{revtex4}
%\documentclass[prb]{revtex4}% Physical Review B


%\documentclass[letterpaper,12pt,onecolumn]{article}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
\usepackage{dcolumn} % to align on a decimal point in columns
\usepackage{amsmath}

%\nofiles

\begin{document}

%\preprint{submitted to {\it Journal of Biomechanical Engineering}}

\title{Simulation and Optimization of Human Gait Using Direct Collocation}

\author{Marko Ackermann}
%\altaffiliation[Also at ]{Department of Biomedical Engineering, Cleveland Clinic Foundation.}
\author{Antonie J. van den Bogert}%
\email{bogerta@ccf.org}
\affiliation{%
Department of Biomedical Engineering, Cleveland Clinic Foundation,
9500 Euclid Avenue (ND-20),
Cleveland, OH 44195, USA}

%date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
Many applications in biomechanics require finding optimal neural controls for a computational musculoskeletal model, either to track measured human performances, or to predict system behavior by optimizing a performance criterion. Current approaches consist of efficient control-based algorithms for the former and shooting methods for the latter problem. In this paper we present a general approach based on direct collocation, a numerical method commonly used for optimal control problems in aeronautics and aerospace engineering. In direct collocation, the controls as well as the states are time discretized and treated as unknowns in a large-scale parameter optimization.  The method was used with a two-dimensional musculoskeletal model to solve the optimal tracking problem and to generate predictive simulations of gait to minimize a weighted sum of squared muscle activations.  The relationship between accuracy and computation time was quantified.  It was found that only 50 time nodes for half a gait cycle allowed simulation of the relevant details of human gait.  At this temporal resolution, perfectly periodic optimal control solutions could be obtained in well under one hour of computation time with standard desktop hardware.  The tracking solution was well within one standard deviation of measured data.  The predictive solution had most features of normal gait but failed to predict knee flexion in mid-stance and had unrealistic impact forces.  These results suggest further research on optimality criteria for human gait. We conclude that this approach is a viable alternative to control-based techniques in solving the optimal neuromuscular tracking problem, and proved to be especially promising in predictive simulations of gait, where it has the potential of substantially reducing computation cost and overcoming many of the numerical problems of the shooting methods that are presently used.
\end{abstract}

%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
                             % Classification Scheme.
\keywords{gait, simulation, musculoskeletal modeling, optimal control}% Use showkeys class option if keyword
                              %display desired
\maketitle

\section{Introduction \label{sec:intro}}

Many theoretical problems in biomechanics require finding optimal neural controls for a computational musculoskeletal model.  Within this framework, there are two main lines of research: neuromuscular tracking, where the goal is to minimize the difference between simulated and measured system behavior, and predictive simulation, where measurements are not used and the optimization is based on performance (e.g.\ maximize jump height) or cost (e.g.\ minimize metabolic energy for a given task).

Neuromuscular tracking is typically applied to the quantitative assessment of muscle coordination and internal loads on articular structures ~\cite{erdemir:mclean:herzog:2007}. In many cases, it might be appropriate to perform the well-established, and efficient inverse dynamics combined with an static optimization to solve the muscle force-sharing problem when kinematics and ground reaction forces are measured (e.g.~\cite{crowninshield:brand:1981}). In some cases, however, this conventional approach might fail in providing realistic muscle force histories that are compatible with muscle capabilities because it does not consider the dynamics involved in muscle force generation (e.g.~\cite{ackermann:schiehlen:2007}). Neuromuscular tracking methods solve this problem by simulating a forward dynamic model of the musculoskeletal system, while minimizing a cost function containing the ``distance'' between the measured and simulated data.  The numerical methods range from global optimization of stimulation parameters (e.g.~\cite{neptune:zajac:kautz:2004}), and more recently, efficient control-based algorithms such as computed muscle control (CMC)~\cite{thelen:anderson:delp:2003, thelen:anderson:2006} and neuromusculoskeletal tracking (NMT)~\cite{seth:pandy:2007}.

Predictive simulations are those that do not require measured human performances.  For gait, optimizations have usually been based on minimization of metabolic energy. This attempts to predict human performance for a given task solely on the basis of musculoskeletal system properties. Predictive simulations can help uncover underlying principles of neuromuscular coordination and have potential applications in predicting patient responses to specific surgical interventions (e.g.\ a tendon transfer procedure), and to prosthetic and orthotic devices.  Other potential applications are in the reconstruction of gait for extinct species~\cite{sellers:cain:wang:2005,sellers:manning:2007}.  The current state of the art in this field is represented by the work of Anderson and Pandy~\cite{anderson:pandy:2001b} who predicted a walking gait using a 3D musculoskeletal model with minimization of metabolic cost of transport.

In biomechanics of human movement, the conventional approach to predictive simulation consists of converting the optimal control problem into a parameter optimization by means of a temporal discretization of the neural excitations. The system is simulated by forward simulation of the musculoskeletal dynamics, i.e.\ by solving the initial value problem. This approach, referred to as shooting method here, requires many numerical integrations of the musculoskeletal system dynamics, usually with small time steps, leading to high computation cost, when using deterministic or stochastic optimization algorithms. Stochastic algorithms, as Simulated Annealing~\cite{neptune:wright:bogert:2000}, require many function calls to conduct a global search, while deterministic algorithms, used e.g.\ by Anderson and Pandy~\cite{anderson:pandy:2001b}, require derivative information which is estimated by finite-difference approximations, leading to the necessity of integrating the system equations for each perturbation of the design variables in order to compute the gradient of the cost function.  Indeed, Anderson and Pandy~\cite{anderson:pandy:2001b} report $10000$ hours of computational time for a walking simulation, Neptune et al.~\cite{neptune:wright:bogert:2000}, $640$ hours for a running simulation, and McLean et al.~\cite{mclean:su:bogert:2003}, $37$ hours for a side stepping simulation. In order to diminish computational cost, modelers are often forced to reduce model complexity or simplify control discretization, reducing the parameters describing the neural excitation histories to a few control points or even simply to an onset time, a magnitude and an offset time. The shooting method is also plagued by various numerical problems.  Because of the unstable nature of walking, the model tends to fall down, making cost function and final state extremely sensitive to changes in initial conditions and to neural excitations \cite{betts:1998}. The nonlinearity and stiffness of foot-ground contact result in a large number of integration steps and, if the swing leg contacts the ground, increased sensitivity of cost function and final state.  A periodic motion using the shooting method is hard to obtain under these circumstances.  For instance, inspection of the results in~\cite{anderson:pandy:2001b} shows that the periodic boundary conditions were not satisfied.

The large computation cost and numerical problems of the shooting method has stimulated recent efforts to develop alternative methods. For instance, Ackermann~\cite{ackermann:2007} and Ren et al.~\cite{ren:jones:howard:2007} propose inverse-dynamics based approaches to simulate gait that represent selected states by interpolating polynomials or Fourier series, respectively. The former method requires, however, the inversion of the musculoskeletal system dynamics, while the latter does not incorporate muscles. Furthermore, the implementation of these methods is difficult and problem-specific.

Direct collocation (DC) is an alternative approach to transform an optimal control problem into a parameter optimization that has been extensively used for trajectory optimization in aeronautics and aerospace engineering~\cite{betts:1998}. The first use of DC to solve the optimal control problem of human movement can be attributed to Kaplan and Heegaard's work on pedaling~\cite{kaplan:2000,kaplan:heegaard:2001,kaplan:heegaard:2002}. In DC, the controls as well as the state variables are discretized on a suitable time grid and treated as a large set of unknowns. With the time discretization, the system differential equations become a large set of algebraic constraints. The resulting large-scale, nonlinear programming problem can be solved by state-of-the-art optimization algorithms designed to cope with the large set of unknowns and to exploit the sparse nature of the resulting Jacobian matrix of the constraints and Hessian matrix of the cost function. Direct collocation offers several potential advantages over the shooting method in solving the optimal neuromuscular coordination problem, in particular for gait: (1) it could reduce computational cost by avoiding the necessity of numerically integrating the system differential equations; the differential equations are implicitly integrated as a whole as the algebraic constraints are fulfilled~\cite{kaplan:2000}; (2) it is better suited to solve the two-point, boundary-value problem arising from the periodicity requirement; and (3) it may present a more favorable scaling of computational time with model complexity~\cite{kaplan:heegaard:2001}, which is important, considering the complex, high-dimensional models of the musculoskeletal system required and currently used to study human locomotion, refer e.g.\ to \cite{delp:loan:2000}.

Despite its potential advantages, DC has not yet been applied to neuromuscular optimal control problems in gait, perhaps because of methodological complexity~\cite{kaplan:heegaard:2001,kaplan:heegaard:2002}.  In this paper we present a general method for optimization of movement using DC, utilizing an existing musculoskeletal model and commercially available optimization software. The method was applied to neuromuscular tracking as well as predictive simulation of gait, and tests were conducted to determine the relationship between computation time and accuracy.


\section{Methods \label{sec:methods}}

\subsection{Musculoskeletal Model \label{sec:musc_skel_mod}}
The model of musculoskeletal dynamics has been described previously~\cite{gerritsen:bogert:hulliger:1998, hardin:su:bogert:2004}. The model of the skeletal system consists of seven rigid body segments and their mass properties: a trunk, and bilateral thighs, shanks, and feet. Joints are assumed to be ideal hinges, and there are no kinematic constraints between the feet and ground, resulting in a total of nine kinematic degrees of freedom. Eight muscle groups are included in each lower extremity: {\it Iliopsoas}, {\it Glutei}, {\it Hamstrings}, {\it Rectus Femoris}, {\it Vasti}, {\it Gastrocnemius}, {\it Soleus}, and {\it Tibialis Anterior}. Each muscle is represented by a $3$-element Hill-type model, using the equations from~\cite{mclean:su:bogert:2003} and muscle properties from~\cite{gerritsen:bogert:hulliger:1998}. This model has $50$ state variables in the state vector $\bm{x}$: $f=9$ generalized coordinates in $\bm{q}$, $f=9$ generalized velocities in $\bm{\dot{q}}$, $m=16$ muscle contractile element (CE) lengths in $\bm{l}_{ce}$, and $m=16$ muscle activations in $\bm{a}$, where $\bm{x} = {[\bm{q}^T \quad \dot{\bm{q}}^T \quad \bm{l}_{ce}^T \quad \bm{a}^T]}^T $. The differential equations describing the dynamics of the musculoskeletal system read as
\begin{subequations}
\label{eq:der_dep}
\begin{eqnarray}
% \nonumber to remove numbering (before each equation)
  \dot{\bm{q}} &=& \bm{I}_{f \times f} \, \dot{\bm{q}} \label{eq:der_dep_dq} \\
  \ddot{\bm{q}} &=& \ddot{\bm{q}}(\bm{q}, \dot{\bm{q}}, \bm{l}_{ce}) \label{eq:der_dep_ddq} \\
  \dot{\bm{l}}_{ce} &=& \dot{\bm{l}}_{ce}(\bm{q}, \bm{l}_{ce}, \bm{a}) \label{eq:der_dep_dl} \\
  \dot{\bm{a}} &=& \dot{\bm{a}}(\bm{a}, \bm{u}) \label{eq:der_dep_da} \, ,
\end{eqnarray}
\end{subequations}
where $\bm{I}$ denotes the identity matrix, Eq.~\ref{eq:der_dep_ddq} contains the equations of motion which were generated by SD/Fast (Parametric Technology Corp., Needham, MA), Eq.~\ref{eq:der_dep_dl} describes the muscles contraction dynamics, and Eq.~\ref{eq:der_dep_da} describes the muscle activation dynamics in which $\bm{u}$ are the $m=16$ control inputs, i.e. the neural excitations for each muscle.


\subsection{Ground Contact Model \label{sec:ground_mod}}

The interaction between feet and ground is modeled by means of $10$ spring-damper elements uniformly distributed along each foot sole. The vertical force for each contact element $j$, was modeled as
\begin{equation}\label{eq:fyj}
    f_{y,j} = a \, \delta_j^3 \, (1 + b \, \dot{\delta}_j) \, ,
\end{equation}
where $\delta_j$ is the ground penetration of element $j$, $a$ is a vertical stiffness parameter set to $a=5.0e7 \, Nm^{-3}$, and $b$ is a vertical damping parameter set to $b=1.0 \, m^{-1}s$.  These parameter values are consistent with test results by~\cite{aerts:declercq:1993}.

The horizontal force at contact element $j$ was modeled by an approximation of Coulomb friction~\cite{bogert:schamhardt:crowe:1989} using a logistic function:
\begin{equation}\label{eq:fxj}
    f_{x,j} = - \frac{1-\exp(-v_{s,j}/v_c)}{1+\exp(-v_{s,j}/v_c)} \, \mu \, f_{y,j} \, ,
\end{equation}
where $f_{y,j}$ is the vertical force at contact element $j$, $v_{s,j}$ is the sliding velocity at element $j$, and $v_c$ is a scaling factor set to $v_c = 0.05 \, m s^{-1}$. The friction coefficient $\mu$ was given a value of 1.0. The approximation is necessary due to numerical reasons, where a lower scaling factor $v_c$ implies a better approximation of true Coulomb friction.


\subsection{Problem Formulation\label{sec:opt_prob}}

The optimization problem for gait consists of searching for time histories of controls $\bm{u}(t)$ and states $\bm{x}(t)$ that minimize a scalar, integral cost function
\begin{equation}\label{eq:cost_gen}
    J = \int_{0}^{T} L(\bm{x}(t),\bm{u}(t)) \, dt \, ,
\end{equation}
where $T$ is the gait cycle period.
%$\bm{x}(t)$ contains the muscle activations $\bm{a}(t)$, the muscle contractile element lengths $\bm{l}_{CE}(t)$,  the generalized coordinates $\bm{q}(t)$ and the generalized velocities $\dot{\bm{q}}(t)$.
The solution must also satisfy the differential equations describing the dynamics of the musculoskeletal system
\begin{equation}\label{eq:dyn_gen}
    \dot{\bm{x}}(t) = \bm{f}(\bm{x}(t),\bm{u}(t),t) \, ,
\end{equation}
represented by Eqs.~\ref{eq:der_dep} in this work, and to additional inequality constraints typically given by bounds on the neural excitations
\begin{equation}\label{eq:bounds_u}
    0 \leq \bm{u}(t) \leq 1 \, ,
\end{equation}
and on the states as
\begin{equation}\label{eq:bounds_x}
    \bm{x}_{min} \leq \bm{x}(t) \leq \bm{x}_{max} \, .
\end{equation}
A periodic gait is ensured by
\begin{subequations}
\label{eq:per_const}
\begin{eqnarray}
% \nonumber to remove numbering (before each equation)
  \bm{\bar{x}}(T) &=& \bm{\bar{x}}(0) \, , \quad \textrm{and} \label{eq:per_const1} \\
  x_h(T) &=& x_h(0) + v_w T \, , \label{eq:per_const2}
\end{eqnarray}
\end{subequations}
where $\bm{\bar{x}}$ contains all the states excepting the horizontal position of the trunk's center of mass $x_h$, and $v_w$ is the walking speed.  In our optimizations, we prescribe $v_w$ as a task constraint, but $T$ is allowed to be optimized.

In this paper two different cost functions are employed, one to solve the optimal tracking problem, and the other to generate predictive simulations of gait that minimize``effort''. Effort is quantified by squared muscle activations weighted by an estimator of muscle volume as
\begin{equation}\label{eq:effort_term}
    J_{e} = \frac{1}{V} \frac{1}{T} \sum_{i=1}^{m} V_i \int_{0}^{T} a^2_i(t) \, dt \, ,
\end{equation}
where $m$ is the number of muscle groups, $a$ is the muscle activation, $V_i$ are weighting factors proportional to muscle volume as $V_i = L^{ce}_{opt,i} F_{max,i}$, where $L^{ce}_{opt}$ is the optimal muscle contractile element length and $F_{max}$ is the maximal muscle isometric force, and $V = \sum{V_i}$. The volume weighting was introduced to take muscle size into consideration, so that activation of a larger muscle will have higher cost than the same activation in a smaller muscle.

In order to solve the optimal tracking problem we adopt the following cost function
\begin{equation}\label{eq:cost_track}
    J_{tr} = J_d + J_e \, ,
\end{equation}
where $J_d$ is a term quantifying the distance of the simulated motion and ground reaction forces to the experimental, target data, while $J_e$, Eq.~\ref{eq:effort_term}, quantifies muscle effort. We quantify deviation from target experimental data as
\begin{equation}\label{eq:dev_term}
    J_{d} = \frac{1}{T} \frac{1}{n_{tr}} \sum_{j=1}^{n_{tr}} \int_{0}^{T} \left(\frac{y_j(t)-\bar{y}_j(t)}{\sigma_j(t)} \right)^2 dt \, ,
\end{equation}
where $n_{tr}$ is the number of tracked variables, $y_j$ is the simulated variable, $\bar{y}_j$ is the mean value of the target experimental value of the variable, and $\sigma_j$ is the experimental standard deviation of the corresponding variable. The effort term, Eq.~\ref{eq:effort_term}, is added to the deviation term in order to resolve the muscle force-sharing problem. The relative importance of the two terms was chosen to favor a good tracking quality, as shown in section~\ref{sec:track_norm}.

\subsection{Collocation Approach \label{sec:dc}}

Direct collocation transforms the optimal control problem, Eqs.~\ref{eq:cost_gen}-\ref{eq:per_const}, into a parameter optimization. It consists of a temporal discretization of the system differential equations, Eq.~\ref{eq:dyn_gen}, transforming them into a large set of algebraic equations. Because walking is symmetric, only half a gait cycle needs to be considered to fully characterize the entire walking cycle. This substantially reduces the computation time, at the expense of excluding asymmetrical gaits such as galloping. The time interval for half a gait cycle was divided into $n$ equally spaced nodes:
\begin{equation}\label{eq:nodal_times}
    0 = t_0 < t_1 < \cdots < t_k < \cdots < t_n = T/2 \, .
\end{equation}
The states $\bm{x}(t)$ and the controls, i.e.\ the neural excitations $\bm{u}(t)$, are accordingly discretized and collected in a vector containing all unknowns:
\begin{equation}\label{eq:opt_vars}
    \bm{\chi} = \left[\, \bm{x}_1^T \, \bm{u}_1^T \, \bm{x}_2^T \, \bm{u}_2^T \, \cdots \, \bm{x}_k^T \, \bm{u}_k^T \, \cdots \, \bm{x}_{n-1}^T \, \bm{u}_{n-1}^T \, \,  T/2 \, \right]^T \, ,
\end{equation}
where $\bm{x}_k = \bm{x}(t_k)$ and $\bm{u}_k = \bm{u}(t_k)$. Note that $\bm{\chi}$ contains also the duration of half a gait cycle $T/2$. Also note that the values of the controls and states at the $n^{th}$ node are not included. This is because the motion is periodic, so that kinematics, CE lengths, muscle activations and neural excitations in one half of the body at the end of the half gait cycle coincide with the ones in the other half of the body at the beginning of the gait cycle, or $\bm{x}_{n} = \bm{x}^*_{1}$ and $\bm{u}_{n}=\bm{u}^*_{1}$, where $\bm{x}^*_k$ and $\bm{u}^*_k$ contain appropriate permutations of the elements of $\bm{x}_k$ and $\bm{u}_k$, respectively, and where the horizontal position of the hip $x_h$ is shifted by $v_w \, T/2$ as $x^*_{h,k} = x_{h,k} + v_w \, T/2$. This operation is referred to as ``mirroring'' throughout this paper.

The system differential equations were discretized using two different schemes ~\cite{betts:huffman:1998, betts:1998}: Euler's discretization scheme (EDS) and the trapezoidal discretization scheme (TDS). The discretization of the system differential equations, Eq.~\ref{eq:dyn_gen}, results in a set of algebraic constraints called {\it defects} $\bm{\Delta}_k$, which, for the EDS case, are
\begin{subequations}
\label{eq:EDS}
\begin{eqnarray}
    \bm{\Delta}_k &=& \bm{x}_k - \bm{x}_{k-1} -  \Delta t \bm{f}_k = 0 \, , \quad k=2 \ldots n-1 \, , \label{eq:EDS1}\\
    \bm{\Delta}_1 &=& \bm{x}^*_1 - \bm{x}_{n-1} -  \Delta t \bm{f}^*_1 = 0 \, , \label{eq:EDS2}
\end{eqnarray}
\end{subequations}
where $\bm{f}_k = \bm{f}(\bm{x}_k,\bm{u}_k,t_k)$, $\Delta t = t_k - t_{k-1}$, and $\bm{f}^*_k$ is an appropriate permutation of $\bm{f}_k$ as done for $\bm{x}^*$. This discretization scheme corresponds to the computation of derivatives by means of a backward finite-difference approximation.  Furthermore, since the motion is periodic with $\bm{x}_{n} = \bm{x}^*_{1}$ and $\bm{u}_{n}=\bm{u}^*_{1}$ in our formulation of the problem, the periodicity constraints, Eqs.~\ref{eq:per_const}, are implicitly included in Eq.~\ref{eq:EDS2}. The TDS reads as
\begin{subequations}
\label{eq:TDS}
\begin{eqnarray}
    \bm{\Delta}_k &=& \bm{x}_k - \bm{x}_{k-1} -  \frac{\Delta t}{2} (\bm{f}_k + \bm{f}_{k-1}) = 0 \, , \quad k=2 \ldots n-1 \, , \label{eq:TDS1}\\
    \bm{\Delta}_1 &=& \bm{x}^*_1 - \bm{x}_{n-1} -  \frac{\Delta t}{2} (\bm{f}^*_1 + \bm{f}_{n-1}) = 0 \, . \label{eq:TDS2}
\end{eqnarray}
\end{subequations}
All algebraic constraints in Eqs.~\ref{eq:EDS} or Eqs.~\ref{eq:TDS} are collected in a vector of constraints
\begin{eqnarray} \label{eq:constrs}
    \bm{\Delta} = \left[ \bm{\Delta}_1^T \, \cdots \, \bm{\Delta}_k^T \, \cdots \, \bm{\Delta}_{n-1}^T \right]^T \, .
\end{eqnarray}

The lower and upper bounds on the neural excitations, Eq.~\ref{eq:bounds_u}, on the states, Eq.~\ref{eq:bounds_x}, and on the duration of half a gait cycle are replaced by bounds on the optimization variables as
\begin{equation}\label{eq:bounds_pars}
    \bm{\chi}_{min} \leq \bm{\chi} \leq \bm{\chi}_{max} \, .
\end{equation}
Notice that, by prescribing the walking velocity $v_w$, the ratio between horizontal displacement of the trunk's center of mass and the duration of half a gait cycle is constrained by Eq.~\ref{eq:EDS2} or Eq.~\ref{eq:TDS2}.

The two cost functions used in this work were discretized using the trapezoidal quadrature rule to estimate the integrals in Eq.~\ref{eq:effort_term} and Eq.~\ref{eq:dev_term}. In the case of uniformly distributed nodes, this leads to the following dicretized version of the ``deviation'' term, Eq.~\ref{eq:dev_term},
\begin{equation}\label{eq:dev_discr}
    J^d_{d} = \frac{1}{n_{tr} n} \sum_{j=1}^{n_{tr}} \sum_{k=1}^{n} \left(\frac{y_{j,k}-\bar{y}_{j,k}}{\sigma_{j,k}} \right)^2 \, ,
\end{equation}
and to the following one for the ``effort'' term, Eq.~\ref{eq:effort_term},
\begin{equation}\label{eq:effort_discr}
    J^d_{e} = \frac{1}{V} \frac{1}{(n-1)} \sum_{i=1}^{m} V_i \sum_{k=1}^{n-1} a^2_{i,k} \, .
\end{equation}
The dicretized version of the cost function for solving the optimal tracking problem, Eq.~\ref{eq:cost_track}, reads as
\begin{equation}\label{eq:track_discr}
    J^d_{tr} = J^d_{d} + J^d_{e} \, .
\end{equation}
The discretized version of the ``effort'' cost function is slightly modified by neglecting the values at the $n^{th}$ node, which coincide with the ``mirrored'' values at the first node. This simplifies the cost function and facilitates the computation and presentation of the gradient and Hessian of the cost function, while having no influence on the outcomes of the analysis presented here.

The ``effort'' cost function is a simple function of system states which allows for the computation of analytical gradient and Hessian, which reduces computation time, as described in Appendix~\ref{app:A}. The constraint Jacobian matrix $\frac{\delta \bm{\Delta}}{\delta \bm{\chi}}$ is also sparse (Appendix~\ref{app:B}) and was calculated through finite-difference approximations.

The SNOPT code for Matlab (Tomlab Optimization Inc., Pullman, WA) was used to solve the optimization problem, which searches for the optimal discretized states and controls in Eq.~\ref{eq:opt_vars}, which minimizes the cost function Eq.~\ref{eq:effort_discr} or Eq.~\ref{eq:track_discr}, subject to the dynamical constraints, Eq.~\ref{eq:constrs}, and to the lower and upper bounds on the states and neural excitations, Eq.~\ref{eq:bounds_pars}. SNOPT is a robust implementation of a state-of-the-art sequential quadratic programming (SQP) algorithm for solving large-scale, sparse nonlinear programming problems~\cite{gill:murray:saunders:2002}. For both optimization problems, the optimal tracking problem and the minimal ``effort'' problem, we provide the sparsity pattern of the constraints Jacobian, Eqs.~\ref{eq:Jac_subs}-\ref{eq:Jac_T} (Appendix~\ref{app:B}), while for the minimal ``effort'' problem we additionally provide analytical gradient, Eq.~\ref{eq:grad_effort}, and Hessian, Eq.~\ref{eq:Hess}, of the cost function (Appendix~\ref{app:A}).  This allows SNOPT to store and manipulate these large matrices efficiently. In order to assess the error due to the discretization we employ the approach proposed by Betts and Huffman~\cite{betts:huffman:1998}, with details in Appendix~\ref{app:C}.

In order to solve the optimal tracking problem we adopt Eq.~\ref{eq:track_discr} as the cost function and prescribe an average walking speed of $v_w = 1.1 \, m/s$. The tracked variables $y_j$ in Eq.~\ref{eq:dev_discr} are the hip, knee and ankle angles, and the vertical and horizontal ground reaction forces, with $n_{tr} = 5$. The target mean experimental values, $\bar{y}_j$, for these variables as well as the corresponding inter-subject standard deviations, $\sigma_j$,  are adopted from~\cite{winter:1991}. The minimal standard deviation for each variable was limited to $20\%$ of the maximal one in order to avoid division by zero in Eq.~\ref{eq:dev_discr}. As mentioned previously, the ``effort'' term, Eq.~\ref{eq:effort_discr}, is added to the cost function, Eq.~\ref{eq:track_discr}, to guarantee a unique solution for the muscle forces, i.e.\ to avoid the occurrence of an undetermined problem for muscle coordination. In order to study the accuracy of the solution, different number of nodes, $n=25, \, 50, \, 99 \, \text{and} \, 197$, as well as two different discretization schemes, the EDS, Eqs.~\ref{eq:EDS}, and the TDS, Eqs.~\ref{eq:TDS}, were used. Partially converged results obtained using the shooting method were used as initial guess for the tracking optimization using EDS and $50$ nodes. The latter solution was used as initial guess for the optimizations using EDS and $25$ nodes and EDS and $99$ nodes. The solution using EDS and $99$ nodes was used as initial guess for the optimization using EDS and $197$ nodes. For all the optimizations using TDS the solutions obtained using the EDS and the same number of nodes were used as initial guesses.

In order to generate predictive simulations of gait using direct collocation Eq.~\ref{eq:effort_discr} is adopted as the cost function measuring ``effort''. Notice that this corresponds to the ``effort'' term of the cost function used to solve the optimal tracking problem. An average walking speed of $1.1 \, m/s$ is prescribed. Because of the increased robustness of the EDS over the TDS, quantified by the number of successfully converged optimizations, we adopt the EDS for the predictive simulations. Here again the accuracy of the solution is studied by using different number of nodes, $n=25, \, 50, \, 99, \, 197 \, \text{and} \, 393$.  The solution for the optimal tracking problem using EDS and $25$ nodes was used as initial guess for the predictive optimization using $25$ nodes. The optimizations using higher number of nodes were successively performed using the solutions with the immediately lower number of nodes as initial guesses.

%%% Valuable literature passages for this section:

%% Kaplan (2000)
% "Collocation techniques, also called transcription, replace the differential equations of the optimal control problem with a set of algebraic constraints involving a large set of unknowns" (pg.30)
% "The integration of the differential equations is implicit in the form of the method lowering computational cost" (pg.30)
% "Direct methods are more robust and adaptable than indirect methods. Direct collocation methods with nonlinear programming have become practical techniques to solve the optimal control problem."(pg.30)
% "The method is well-suited to dealing with inequality constraints on the state and control variables since the arcs sequences are not required a priori." (pg.30)
% "The differential equations can be then turned into algebraic equality constraints, also called defects, $\delta$, in the unknown nodal state and control variable values," (pg.30)
% "Betts [21,22] compared three different discretization techniques for trajectory optimization problems. The collocation methods were trapezoidal, Hermite-Simpson and Runge-Kutta. The local truncation errors for these methods are respectively O(\delta t^3), O(\delta t^5)m and O(\delta t^5), with errors over the global While the accuracy of the trapezoidal was certainly less than that of the other methods, the number of function evaluations was lower and the method was more robust, meaning convergence could be achieved for significantly poorer initial guesses." (pg.34)
% We consider these methods the most fitting for biomechanical optimal control problems" (pg.34)

%% Kaplan (2001)
% "we use a direct collocation method to discretize the differential equations. Direct collocation methods together with nonlinear programming have been effectively applied to optimal control problems in the field of trajectory optimization.", (pg.1078)

% neuromuscular control problem

%% Kaplan (2002)
% The equations of motion are never explicitly integrated; rather, satisfying the equality constraint equations is equivalent to integrating the differential equations." (pg. 2058)

%% Betts (1998)
% Issues in Direct Shooting:
% " Changes early in the trajectory propagate to the end of the trajectory. The net effect is that is that constraints can become very nonlinear with respect to variables, thereby making the optimization problem difficult to solve." (pg.200)
% " ... computational cost of evaluating the gradients information. A forward difference approximation ..." (pg.200)
% " ... small changes introduced early in the trajectory can propagate into very nonlinear changes at the end of the trajectory", (202)
% " Historically, transcription or collocation methods were developed for solving two-point boundary value problems such as those encountered when solving the necessary conditions with an indirect formulation." (pg.203)
% " the nonlinear programming problem that results from this formulation is large. ... Fortunately, the pertinent matrices for this NLP problem, namely, the Jacobian G and the Hessian H, are also sparse. ... . Consequently, exploiting sparsity to reduce both storage and computational time is a critical aspect of a successful implementation when using a direct transcription method" (pg.204)
% "The sparse NLP implements a sparse SQP method based on a Schur-complement algorithm suggested by Gill el at.. Jacobian and Hessian matrices are computed efficiently using sparse finite differencing as proposed by Coleman and More and Curtis et al." (pg.204)


\section{Results \label{sec:results}}

\subsection{Tracking of Normal Human Walking \label{sec:track_norm}}

In this section, the results obtained by solving the optimal tracking problem using direct collocation are presented. The left hand side of Fig.~\ref{fig:cost_error_track} shows the value of the tracking cost function, Eq.~\ref{eq:track_discr}, as well as the values of the deviation term, Eq.~\ref{eq:dev_discr}, and the ``effort'' term, Eq.~\ref{eq:effort_discr}, for the different number of nodes. The values of the cost function and of the corresponding terms tend to converge as the number of nodes increases. The curves suggest that a discretization using $50$ nodes for half a gait cycle is sufficiently fine for this specific problem.

%%%%%%%%%%%%%%%%%%%%%%%%
% To report values of cost function for different number of nodes choose either table or diagram
%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{table}
%  \caption{Cost function value, Eq.~\ref{eq:??}, as a function of the number of nodes for optimal tracking simulations of normal walking at $1.1$~m/s performed using both the Euler discretization scheme (EDS) and the trapezoidal discretization scheme (TDS). \label{tab:cost_track}}
%  \begin{ruledtabular}
%  \begin{tabular}{ccc}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
%           & \multicolumn{2}{c}{Cost Function} \\
%  \# nodes & EDS & TDS \\ \hline
%  25       & 0.1562 & 0.1490 \\
%  50       & 0.1459 & 0.1472 \\
%  99       & 0.1437 & 0.1434 \\
%  197      & 0.1427 & 0.1437 \\
%  \end{tabular}
%  \end{ruledtabular}
%\end{table}

%\begin{figure}
  % Requires \usepackage{graphicx}
  % Figure originally generated using the Matlab-function plot_cost.m and exporting the figure as eps setting the Width to 8cm, fixing the text font size to 10, and properly offsetting the figure objects
%  \includegraphics[width=8cm]{figs/cost_terms}\\
%  \caption{Cost function, Eq.~\ref{eq:track_discr}, and the corresponding deviation term, Eq.~\ref{eq:effort_discr}, and effort term, Eq.~\ref{eq:effort_discr}, obtained using different number of nodes and the EDS and TDS discretization schemes.}\label{fig:cost_terms_track}
%\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%
% To report relative errors for different number of nodes choose either table or diagram
%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{table}
%  \caption{Average (avg) and maximal (max) relative discretization error as defined in Eq.~\ref{eq:??} and %Eq.~\ref{eq:??}, respectively, as a function of the number of nodes for optimal tracking simulations of normal gait at $1.1$~m/s performed using both the Euler discretization scheme (EDS) and the trapezoidal discretization scheme (TDS). \label{tab:err_track}}
%  \begin{ruledtabular}
%  \begin{tabular}{ccdcd}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
%             & \multicolumn{4}{c}{Relative Discretization Error} \\
%             & \multicolumn{2}{c}{EDS}& \multicolumn{2}{c}{TDS} \\
%  \# nodes   &   \multicolumn{1}{c}{avg}     &   \multicolumn{1}{c}{max}   &    \multicolumn{1}{c}{avg}     &    \multicolumn{1}{c}{max}     \\ \hline
%  25         &   1.3e-3   &   1.60    &   8.1e-4   &   0.97     \\
%  50         &   1.0e-4   &   0.25    &   6.2e-5   &   0.15     \\
%  99         &   2.0e-5   &   0.10    &   7.9e-6   &   0.039    \\
%  197        &   2.3e-6   &   0.023   &   3.9e-7   &   0.0038   \\
%  \end{tabular}
%  \end{ruledtabular}
%\end{table}

%\begin{figure}
  % Requires \usepackage{graphicx}
  % Figure originally generated using the Matlab-function plot_errors.m and exporting the figure as eps setting the width to 8cm, fixing the text font size to 10, and properly offsetting the figure objects
%  \includegraphics[width=8cm]{figs/rel_error}\\
%  \caption{Relative discretization errors estimated using Eq.~\ref{eq:error_rel} for different number of nodes and the EDS and TDS discretization schemes. \label{fig:rel_error_track}}
%\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}
  \begin{tabular}{cc}
  \includegraphics[width=8cm]{figs/cost_terms} &
  \includegraphics[width=8cm]{figs/rel_error}\\
  \end{tabular}
  \caption{Cost function, Eq.~\ref{eq:track_discr}, and the corresponding deviation term, Eq.~\ref{eq:dev_discr}, and effort term, Eq.~\ref{eq:effort_discr}, obtained using the EDS and TDS on the left. Relative discretization errors estimated using Eq.~\ref{eq:error_rel} for different number of nodes and the EDS and TDS.  \label{fig:cost_error_track}}
\end{figure}

The right hand side of Fig.~\ref{fig:cost_error_track} shows the maximal and average relative errors estimated using Eq.~\ref{eq:error_rel}. The error decreases as the number of nodes increases, while the decrease is faster for the higher-order TDS compared to the EDS. Notice that the maximal error over all states and nodes is about two orders of magnitude higher than the average relative  discretization error. The maximal discretization error is mostly very localized and occurs consistently at the {\it defect} associated with the ankle angle acceleration. This occurs because of the low mass of the feet and because of the strongly nonlinear contact model, both directly influencing the computation of the {\it defects} related to the ankle angle accelerations.

Figure~\ref{fig:stick_track} shows a stick-figure of the overall model kinematics for a complete gait cycle, while Fig.~\ref{fig:kin_grf_track} shows the close agreement between the simulated and target kinematics and ground reaction forces. The solid, black lines represent the simulated kinematics, the solid, grey lines show the mean values of the target experimental data, and the dashed, grey lines represent one standard deviation from the mean values in~\cite{winter:1991}. The simulated results are within one standard deviation of the target experimental data throughout the whole gait cycle.

\begin{figure}
  % Requires \usepackage{graphicx}
  % Figure generated using the Matlab-function stick_plot.m and exporting the figure as eps setting the width to 8cm and then fixing the text font size to 10 (use 'result_bk_07_50_3.mat' as input)
  \includegraphics[width=8cm]{figs/stick_bk_07_50_3}\\
  \caption{Snapshots of a stick-figure of the model for a complete gait cycle obtained by solving the optimal tracking problem using $50$ nodes and the EDS. Time between snapshots is $0.0854 \, s$.}\label{fig:stick_track}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%
% For double column (final version) use the first pair of figures, for single column paper (preprint) use the single figure
%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{figure}
%  % Requires \usepackage{graphicx}
%  % Figure generated using the Matlab-function plot_kin_grf_act.m and exporting the figure as eps setting the %width to 8cm, the height to 12cm, and fixing the text font size to 10 (use 'result_bk_07_50_3.mat' as input)
%  \includegraphics[width=8cm]{figs/kin_track}\\
%  \caption{Hip, knee and ankle angles obtained solving the optimal tracking problem using $50$ nodes and the Euler discretization scheme, compared with target %experimental data of Winter~\cite{winter:1991}.}\label{fig:kin_track}
%\end{figure}

%\begin{figure}
%  % Requires \usepackage{graphicx}
%  % Figure generated using the Matlab-function plot_kin_grf_act.m and exporting the figure as eps setting the %width to 8cm, the height to 12cm, and fixing the text font size to 10 (use 'result_bk_07_50_3.mat' as input)
%  \includegraphics[width=8cm]{figs/grf_track}\\
%  \caption{Vertical and horizontal ground reaction forces obtained by solving the optimal tracking problem using $50$ nodes and the Euler discretization scheme, compared with target experimental data of Winter~\cite{winter:1991}.}\label{fig:grf_track}
%\end{figure}

\begin{figure}
  % Requires \usepackage{graphicx}
  % Figure generated using the Matlab-function plot_kin_grf_act.m and exporting the figure as eps setting the width to 8cm, the height to 12cm, and fixing the text font size to 10 (use 'result_bk_07_50_3.mat' as input)
  \begin{tabular}{cc}
  \includegraphics[width=8cm]{figs/kin_track} &
  \includegraphics[width=8cm]{figs/grf_track}\\
  \end{tabular}
  \caption{Hip, knee and ankle angles on the left, and vertical and horizontal ground reaction forces on the right obtained by solving the optimal tracking problem using $50$ nodes and the EDS, compared with the target experimental data in~\cite{winter:1991}.}\label{fig:kin_grf_track}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%

Figure~\ref{fig:acts_track} shows the muscle activation patterns predicted by solving the optimal tracking problem using Eq.~\ref{eq:track_discr} as cost function. Notice that the uniarticular, antagonistic muscles, the pairs {\it Tibialis Anterior} \& {\it Soleus} and {\it Glutei} \& {\it Iliopsoas}, present some co-contraction at onset and offset times, but the co-contraction is discouraged by the presence of the ``effort'' term, Eq.~\ref{eq:effort_discr}, in the cost function.

\begin{figure}
  % Requires \usepackage{graphicx}
  % Figure generated using the Matlab-function plot_kin_grf_act.m and exporting the figure as eps setting the width to 8cm, the height to 12cm, and fixing the text font size to 10 (use 'result_bk_07_50_3.mat' as input)
  \includegraphics[width=8cm]{figs/acts_track}\\
  \caption{Muscle activations obtained by solving the optimal tracking problem using $50$ nodes and the EDS.}\label{fig:acts_track}
\end{figure}


\subsection{Predictive Simulation of Human Walking \label{sec:pred_sim}}

In this section the results of predictive simulations using direct collocation obtained by minimization of ``effort'' are presented. The results for the cost function value, computation time, computation time per major iteration, and average relative dicretization error are summarized in Tab.~\ref{tab:cce_effort}. The value of the cost function converges as the number of nodes increases. Note that the predictive simulation with $25$ nodes for half a gait cycle, using the solution for the optimal tracking problem as initial guess, required only about $10.2$ minutes to converge to substantially different kinematics and muscle activation patterns, see Fig.~\ref{fig:kin_grf_effort}. Figure~\ref{fig:stick_effort} shows the overall kinematics of the model for a complete gait cycle, and Fig.~\ref{fig:kin_grf_effort} shows the predicted kinematics and ground reaction forces compared with experimental data for normal gait in~\cite{winter:1991}. Minimization of ``effort'', in opposition to deviation from experimental data leaded to a motion and ground reaction forces that are substantially different from normal patterns, including a straight-leg pattern during the stance phase, and high impact ground reaction forces following heel strike. The possible explanation and implications of these findings will be briefly discussed in the next section.

\begin{table}
  \caption{Cost function value, Eq.~\ref{eq:effort_discr}, total computational time, computational time per major iteration, and relative average discretization error, Eq.~\ref{eq:error_rel}, as a function of the number of nodes for predictive simulations of normal gait at $1.1 \, m/s$ performed using the EDS. \label{tab:cce_effort}}
  \begin{ruledtabular}
  \begin{tabular}{ccdcc}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  \# nodes    & cost function & \multicolumn{1}{c}{comp.\ time (hs)} & comp.\ time (hs/iter) &  rel.\ discr.\ error     \\ \hline
  25         &    0.01191    &        0.17     &      5.2e-4        &         2.7e-3    \\
  50         &    0.01083    &        0.079    &      2.1e-3        &         2.2e-4    \\
  99         &    0.01029    &        0.92     &      1.7e-2        &         4.1e-5    \\
  197        &    0.00980    &        1.1      &      1.1e-1        &         7.4e-6    \\
  393        &    0.00975    &        9.7      &      3.7e-1        &         1.5e-6    \\
  \end{tabular}
  \end{ruledtabular}
\end{table}

\begin{figure}
  % Requires \usepackage{graphicx}
  % Figure generated using the Matlab-function stick_plot.m and exporting the figure as eps setting the Width to 8cm and the sixing the text font size to 10 (use 'result_bk_06_50_from25.mat' as input)
  \includegraphics[width=8cm]{figs/stick_bk_06_50_from25}\\
  \caption{Snapshots of a stick-figure of the model for a complete gait cycle obtained by minimizing ``effort'' using $50$ nodes and the EDS. Time between snapshots is $0.0892 \, s$}\label{fig:stick_effort}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%
% For double column (final version) use the first pair of figures, for single column paper (preprint) use the single figure
%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{figure}
%  % Requires \usepackage{graphicx}
%  % Figure generated using the Matlab-function plot_kin_grf_act.m and exporting the figure as eps setting the %width to 8cm, the height to 12cm, and fixing the text font size to 10 (use 'result_bk_06_50_from25.mat' as input)
%  \includegraphics[width=8cm]{figs/kin_effort}\\
%  \caption{Hip, knee and ankle angles obtained by minimizing effort as defined in Eq.~\ref{eq:???} using $50$ nodes and the Euler discretization scheme, compared with experimental data of Winter~\cite{winter:1991}.}\label{fig:kin_effort}
%\end{figure}

%\begin{figure}
%  % Requires \usepackage{graphicx}
%  % Figure generated using the Matlab-function plot_kin_grf_act.m and exporting the figure as eps setting the %width to 8cm, the height to 12cm, and fixing the text font size to 10 (use 'result_bk_06_50_from25.mat' as input)
%  \includegraphics[width=8cm]{figs/grf_effort}\\
%  \caption{Vertical and horizontal ground reaction forces obtained by minimizing effort as defined in Eq.~\ref{eq:???} using $50$ nodes and the Euler discretization scheme, compared with experimental data of Winter~\cite{winter:1991}.}\label{fig:grf_track}
%\end{figure}

\begin{figure}
  % Requires \usepackage{graphicx}
  % Figure generated using the Matlab-function plot_kin_grf_act.m and exporting the figure as eps setting the width to 8cm, the height to 12cm, and fixing the text font size to 10 (use 'result_bk_06_50_from25.mat' as input)
  \begin{tabular}{cc}
  \includegraphics[width=8cm]{figs/kin_effort} &
  \includegraphics[width=8cm]{figs/grf_effort}\\
  \end{tabular}
  \caption{Hip, knee and ankle angles on the left, and vertical and horizontal ground reaction forces on the right obtained by minimizing ``effort'', Eq.~\ref{eq:effort_discr}, using $50$ nodes and the EDS, compared with experimental data in~\cite{winter:1991}. \label{fig:kin_grf_effort}}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%

Figure~\ref{fig:acts_effort} shows the activations predicted for the modeled muscle groups. The activation levels are substantially lower than the ones predicted for the optimal tracking problem in Fig.~\ref{fig:acts_track}. Very low activation levels are observed for the knee extensors, {\it Vasti} and {\it Rectus Femoris}. While high muscle activations are necessary for a close track of the target data in the optimal tracking problem, the minimization of ``effort'' leads to a modification of the motion aimed at lowering muscle activations.

\begin{figure}
  % Requires \usepackage{graphicx}
  % Figure generated using the Matlab-function plot_kin_grf_act.m and exporting the figure as eps setting the width to 8cm, the height to 12cm, and fixing the text font size to 10 (use 'result_bk_06_50_from25.mat' as input)
  \includegraphics[width=8cm]{figs/acts_effort}\\
  \caption{Muscle activations obtained by minimizing ``effort'', Eq.~\ref{eq:effort_discr}, using $50$ nodes and the EDS discretization scheme.}\label{fig:acts_effort}
\end{figure}


\section{Discussion \label{sec:diss}}

The results reported in section~\ref{sec:track_norm} demonstrate the applicability of direct collocation and nonlinear programming in solving the optimal tracking problem. The experimental data was tracked within one standard deviation throughout the whole gait cycle, refer to Fig.~\ref{fig:kin_grf_track}. It should be noted that the target experimental data was adopted from the literature and the model used does not necessarily fit the anthropometric features of the subjects. This could explain the overall high muscle activation levels required, Fig.~\ref{fig:acts_track}. These could be reduced by increased weighting of the effort term in the cost function, at the expense of increased tracking error. The computation time of $4.6$ hours required for convergence using $50$ nodes and the EDS is substantially lower than the required for solving the dynamic optimization using the shooting method with a comparable musculoskeletal model. For instance, McLean et al.~\cite{mclean:su:bogert:2003} report a computation time of $37$ hours for tracking simulations of a sidestep movement. Our approach is, however, computationally more costly than computed muscle control (CMC). Thelen and Anderson~\cite{thelen:anderson:2006} report a computation time of about $30$ minutes for walking simulations using a 3D model of the musculoskeletal system with more degrees of freedom and muscles than the model adopted here. In spite of this, direct collocation will be attractive for applications where the limitations of CMC are not acceptable or desirable. The limitations of CMC include the necessity of using a static cost function to solve for muscle redundancy, the necessity of choosing appropriate controller gains, which have no apparent physical interpretation, and the lack of straightforward tuning of the weighting of tracking and effort terms in the cost function. Moreover, CMC does not allow for the use of ground contact models, requiring the measurement of ground reaction forces, and this makes it impossible to use a tracking simulation as a starting point for predictive simulations. The neuromuscular tracking (NMT) method~\cite{seth:pandy:2007} overcomes the latter limitation.

The advantages of direct collocation become striking in predictive simulations of walking. The low computation times shown in Table~\ref{tab:cce_effort} are orders of magnitude lower than typical applications of the shooting method in predictive simulations of gait. Anderson and Pandy~\cite{anderson:pandy:2001b} report $10000$ hours of computation time for optimization of half a walking cycle, which can partly be explained by their more complex model, but we suggest most of the difference is due to the numerical method. Direct collocation effectively avoids many of the numerical pitfalls of the shooting method, arising from the inherent instability of gait, and due to strong nonlinearity of musculoskeletal models and their interaction with the environment which causes forward simulation to proceed with excessively small time steps. Direct collocation is also more suitable to the solution of the two-point boundary value problem arising from the periodicity constraints. The solution shown in Figs.~\ref{fig:stick_track}-\ref{fig:acts_track} is perfectly periodic.

Although already substantially more efficient than the shooting method, the implementation of the direct collocation presented here can be further improved to reduce computation time. The discretization mesh can be locally refined, instead of globally by seeking a distribution of nodes that minimizes maximal errors~\cite{betts:huffman:1998}. Another potential improvement lies in providing analytical expressions for the gradient and Hessian of the cost function and Jacobian of constraints. The structure of the ``effort'' cost function, Eq.~\ref{eq:effort_discr}, already allowed analytical expressions for gradient and Hessian.  This explains the faster computation times compared to the tracking optimizations in section~\ref{sec:track_norm}, for which the gradient and Hessian were numerically approximated. The tracking term of the cost function could be written analytically, with some effort, if necessary.

Analytical expressions for the constraints Jacobian would also improve the scaling of computation time with model complexity. Kaplan and Heegaard~\cite{kaplan:heegaard:2001} show that in a pedaling simulation the scaling of solution times with model complexity favors direct collocation over shooting, which suggests that our methods can be applied to more realistic and three dimensional musculoskeletal models. However, the derivation of analytical expressions for the constraints Jacobian is difficult. It requires symbolic manipulation of the large system differential equations, Eqs.~\ref{eq:der_dep}, and can be extremely tedious and problem-specific for large-scale models. Numerical estimation of derivatives by finite differences, although less accurate and more costly, is the more convenient and flexible approach, but complex step derivatives (e.g.\ \cite{martins:sturdza:alonso:2003}) or automatic differentiation (e.g.\ \cite{bischof:carle:corliss:1992}) might be viable alternatives.

Another interesting possibility permitted by direct collocation is the use of a very low number of nodes. The implicit integration allows the calculations to "step over" unimportant details, unlike in explicit forward simulation methods which tend to perform extremely small time steps in musculoskeletal models. This is illustrated by a computation time of about $10$ minutes for the predictive simulation performed using $25$ nodes and the EDS (Table~\ref{tab:cce_effort}). This is an important feature which allowed us to start with a coarse mesh which is then successively refined, and is also important for applications where computation speed (interactivity) is more important than accuracy, for example as in physics-based character animation. Two decades ago, Witkin and Kass~\cite{witkin:kass:1988}, already developed a direct collocation approach for character animation, for which they used the term ``space-time constraints''.

A practical limitation of direct collocation is the necessity of an initial guess.  Based on our experience, it is important that the initial guess is close to satisfying the system dynamics equations. If the initial guess is not sufficiently dynamics consistent the nonlinear programming algorithm might encounter convergence problems. The primordial initial guess adopted here for finding the first dynamically consistent initial guess was obtained by using the shooting method, which failed to solve the optimization problem but still produced a dynamically consistent solution. Instead of using the shooting method, we did find initial guesses that leaded to convergence by ``trial and error'', but this strategy is not straightforward. An alternative strategy, not tried by the authors yet, consists in ``deactivating'' the periodicity constraints and using a single forward simulation followed by a direct collocation optimization where the periodicity constraint is added as a penalty term to the cost function. Once one has a dynamically consistent solution, this can be used as initial guess for all the subsequent optimizations. Large differences in body anthropometric features or muscle parameters might prevent the use of a previously obtained solution as a sufficiently dynamically consistent initial guess. In those cases the authors successfully used a strategy consisting in a gradual change in musculoskeletal properties combined with successive re-optimizations until the desired change is reached.  As in all gradient-based methods, such solutions are only known to be locally optimal and it is therefore advised to solve the same optimization problem via different initial guesses.

Occasional ``crashing'' of the SQP algorithm was observed, probably due to the high nonlinearity of the contact model. Simply restarting the optimization was mostly sufficient for convergence. A less nonlinear contact model reduced the frequency of such events without changing the outcomes significantly and this may be a good strategy to improve convergence before switching to a more realistic contact model. A fine tuning of the several parameters in the SQP solver could potentially prevent these occurrences but was not extensively performed. According to the authors' experience, in order to decrease computational cost and increase the likelihood of convergence in a predictive simulation the most effective strategy is performing the optimization first with a low number of nodes, say $25$, and using a low-order discretization scheme, the EDS, and then increase number of nodes or switch to higher order discretization schemes to achieve the desired accuracy.

The remarkable differences between the gait patterns predicted by minimizing ``effort'' and standard normal gait patterns deserve further attention. Figures~\ref{fig:stick_effort} and \ref{fig:kin_grf_effort} reveal some important deviations from normal gait patterns, specially a straight-leg pattern during the stance phase, in opposition to a knee flexion of about $20^\circ$ in normal gait. The straight leg at ``touch down'' might explain the high impact ground reaction forces following {\it heel strike}, which are also not consistent with normal gait patterns. This might be a side effect of reducing muscle activations at the knee extensors {\it Vasti} and {\it Rectus Femoris}, Fig.~\ref{fig:acts_effort}, both otherwise active during the weight acceptance phase in normal gait. If these predictions are valid, the minimal energy hypothesis in walking would have to be reconsidered. While energy minimization is undoubtedly an important performance criterion as corroborated by many experimental studies (e.g.\ \cite{bertram:ruina:2001}), other criteria such as reduction of peak joint forces, might play an important role even during normal walking. The possibly oversimplified model, and cost function, and the possible incompatibilities between the model adopted here and the tracked experimental data prevent, however, a more assertive conclusion. A deeper analysis of this topic does not belong to the scope of this paper but is the subject of ongoing research by the authors. Furthermore, although the results of the predictive simulations substantially reduced muscle effort, we cannot be sure that the solution found is a global minimum. Different local minima might exist, and a complete analysis would require exploring the solution space using several different initial guesses.

%\section{Conclusion \label{sec:concl}}
Direct collocation allows the solution of optimal control problems for human movement without excessive computation cost. We expect that this method will contribute to applications of predictive simulation in clinical applications such as the computer-aided design of surgical interventions, prostheses, and other mechanical interventions which aim to alter human gait.

\acknowledgements{This study is supported by the National Institutes of Health grant 1 R01 EB006735-01.}


\appendix

\section{Gradient and Hessian \label{app:A}}
% 5) explain general structure of gradient and sparse Hessian for effort cost function.
 The gradient of the the ``effort'' cost function, Eq.~\ref{eq:effort_discr}, is composed by sub-vectors corresponding to each node $k$ as
\begin{equation}\label{eq:grad_effort}
   \bm{G} = \partial J^d_e / \partial \bm{\chi} = [\bm{G}_1^T \, \cdots \, \bm{G}_k^T \, \cdots \, \bm{G}_{n-1}^T \quad \partial J^d_e / \partial T]^T \, ,
\end{equation}
where each sub-vector reads as
\begin{eqnarray}\label{eq:grad_effort_sub}
    \bm{G}_k = \frac{\partial J^d_e}{\partial \bm{\chi}_k} &=& \left[ \left(\frac{\partial J^d_e}{\partial \bm{q}_k}\right)^T \quad \left(\frac{\partial J^d_e}{\partial \dot{\bm{q}}_k}\right)^T \quad  \left(\frac{\partial J^d_e}{\partial \bm{l}_{ce,k}}\right)^T \quad \left(\frac{\partial J^d_e}{\partial \bm{a}_k}\right)^T \quad \left(\frac{\partial J^d_e}{\bm{u}_k}\right)^T \right]^T = \\
    &=& \left[ \bm{0}_{1 \times f} \quad \bm{0}_{1 \times f} \quad \bm{0}_{1 \times m} \quad \left(\frac{\partial J^d_e}{\partial \bm{a}_k}\right)^T \quad \bm{0}_{1 \times m} \right]^T \, ,
\end{eqnarray}
with
\begin{equation}\label{eq:grad_term}
    \frac{\partial J^d_e}{\partial \bm{a}_k} =  \frac{1}{V} \frac{2}{(n-1)} \left[ \, V_1 a_{1,k} \, \cdots \,  V_i a_{i,k} \, \cdots \,  V_{m} a_{m,k} \, \right]^T \, ,
\end{equation}
and
\begin{equation}\label{eq:der_T}
    \frac{\partial J^d_e}{\partial T} = 0.
\end{equation}
The Hessian for the chosen ``effort'' cost function, $J^d_e$, is a diagonal matrix with constant elements as
\begin{equation}\label{eq:Hess}
    \bm{H} = \frac{\partial^2 J^d_e}{\partial \bm{\chi} \partial \bm{\chi}} = \frac{\partial \bm{G}}{\partial \bm{\chi}} = diag\{[\bm{h}_1 \, \cdots \, \bm{h}_k \, \cdots \, \bm{h}_{n-1} \; 0]\} \, ,
\end{equation}
where
\begin{equation}\label{eq:Hess_sub}
    \bm{h}_k = \frac{1}{V} \frac{2}{(n-1)} \left[ \, \bm{0}_{1 \times f} \quad \bm{0}_{1 \times f} \quad \bm{0}_{1 \times m} \quad V_1 \, \cdots \, V_i \, \cdots \, V_{m} \quad \bm{0}_{1 \times m} \, \right] \, ,
\end{equation}

\section{Constraints Jacobian \label{app:B}}
% 6) explain general structure of constraints Jacobian;
Considering the subdivision of the constraints vector $\bm{\Delta}$ as in Eq.~\ref{eq:constrs} and of the optimization variables vector $\bm{\chi}$ as in Eq.~\ref{eq:opt_vars}, the sparse, nearly block-diagonal constraints Jacobian $\bm{C} = \partial \bm{\Delta} / \partial \bm{\chi}$ arising from the use of the EDS or TDS are formed by the following sub-matrices
\begin{equation}\label{eq:Jac_subs}
    \left( \frac{\partial \bm{\Delta}_p}{\partial \bm{\chi}_l} \right)_{(2f+2m) \times (2f+3m)} =
        \begin{cases}
        \frac{\partial \bm{\Delta}_k}{\partial \bm{\chi}_k} \, ,     & \quad \text{for  } p=l=k   \\
        \frac{\partial \bm{\Delta}_k}{\partial \bm{\chi}_{k-1}} \, , & \quad \text{for  } p=l+1=k \\
        \bm{0} \, ,                                                  & \quad \text{elsewhere} \, ,
        \end{cases}
\end{equation}
and
\begin{equation}\label{eq:Jac_T}
    \frac{\partial \bm{\Delta}}{\partial T} = \bm{1}_{(n-1)(2f+2m)  \times 1} \, ,
\end{equation}
where $\bm{1}$ denotes a matrix consisting of ones.

Considering the dependencies of the states derivatives on the states and controls in the specific model adopted, Eq.~\ref{eq:der_dep}, the sub-matrices of the constraints Jacobian in Eq.~\ref{eq:Jac_subs} have the following sparsity patterns, where ``$1$'' denotes dependency of the state derivative on the corresponding state or neural excitation, i.e.\ a nonzero element in the constraints Jacobian, and ``$0$'' denotes no dependency, i.e.\ a null element in the constraints Jacobian,
\begin{equation} \label{eq:Jac_sub1}
sparsity \left\{\frac{\partial \bm{\Delta}_k}{\partial \bm{\chi}_k} \right\} =
\left[ \begin{array}{ccccc}
\bm{I}_{f \times f} & \bm{I}_{f \times f} & \bm{0}_{f \times m} & \bm{0}_{f \times m} & \bm{0}_{f \times m} \\
\bm{1}_{f \times f} & \bm{1}_{f \times f} & \bm{1}_{f \times m} & \bm{0}_{f \times m} & \bm{0}_{f \times m} \\
\bm{A}_{m \times f} & \bm{0}_{m \times f} & \bm{I}_{m \times m} & \bm{I}_{m \times m} & \bm{0}_{m \times m} \\
\bm{0}_{m \times f} & \bm{0}_{m \times f} & \bm{0}_{m \times m} & \bm{I}_{m \times m} & \bm{I}_{m \times m}
\end{array} \right] \, ,
\end{equation}
for both the EDS and the TDS, where $\bm{A}$ is a sparse sub-matrix whose nonzero elements are defined by the degrees of freedom ``actuated'' by the corresponding muscles, while
\begin{equation} \label{eq:Jac_sub2_EDS}
sparsity \left\{\frac{\partial \bm{\Delta}_k}{\partial \bm{\chi}_{k-1}} \right\} = \left[ \, \bm{I}_{(2f+2m) \times (2f+2m)} \quad \bm{0}_{(2f+2m) \times m} \, \right]
\end{equation}
for the EDS, and
\begin{equation} \label{eq:Jac_sub2_TDS}
sparsity \left\{\frac{\partial \bm{\Delta}_k}{\partial \bm{\chi}_{k-1}} \right\} = sparsity \left\{\frac{\partial \bm{\Delta}_k}{\partial \bm{\chi}_k} \right\}
\end{equation}
for the TDS. Notice that, for $k=1$, $\frac{\partial \bm{\Delta}_k}{\partial \bm{\chi}_{k-1}}:=\frac{\partial \bm{\Delta}_1}{\partial \bm{\chi}_{n-1}}$ and the proper permutation of lines has to be performed according to the ``mirroring'' operation.

\section{Error Estimation \label{app:C}}

In order to assess the error due the discretization we employ the approach proposed by Betts and Huffman~\cite{betts:huffman:1998}. The first step in utilizing this approach is assuming a continuous approximation of the solution for the states $\bm{x}$ and for the neural excitations $\bm{u}$ as $\tilde{\bm{x}}(t)$ and $\tilde{\bm{u}}(t)$, respectively. The states at the nodes are interpolated using cubic Hermite splines so that the values and the first time derivatives of the approximating functions at the nodal times match the nodal values and derivatives as
\begin{eqnarray}\label{eq:approx_x}
    \tilde{\bm{x}}(t_k) &=& \bm{x}_k \quad \text{and} \\
    \frac{d }{dt} \tilde{\bm{x}}(t_k) &=& \bm{f}_k = \bm{f}(\bm{x}_k,\bm{u}_k,t_k) \, ,
\end{eqnarray}
for $k = 1 \ldots n$. The neural excitations at the nodes are linearly interpolated so that the approximated neural excitations at the nodal times match the computed nodal values as
\begin{equation}\label{eq:approx_u}
    \tilde{\bm{u}}(t_k) = \bm{u}_k, \quad \text{for} \quad k=1 \ldots n \, .
\end{equation}

The idea is to compute a better estimate of the states at the subsequent node by halving the time step and using the discretization scheme, Eqs.~\ref{eq:EDS} or Eqs.~\ref{eq:TDS}, twice instead of once. The difference between the states $\bm{\hat{x}}_k$ at the nodal time $t_k$ computed in this way and the nodal states $\bm{x}_k$ is used as an estimator of the error due to the discretization. This procedure leads to the following expressions for the improved estimate of the states at $t_k$,
\begin{equation}\label{eq:est_x_EDS}
    \bm{\hat{x}}_k = \bm{x}_{k-1} + \frac{\Delta t}{2} \left( \bm{f}(\tilde{\bm{x}}_{k-1/2}, \tilde{\bm{u}}_{k-1/2}) + \bm{f}_k  \right)
\end{equation}
for the EDS, and
\begin{equation}\label{eq:est_x_TDS}
    \bm{\hat{x}}_k = \bm{x}_{k-1} + \frac{\Delta t}{4} \left(\bm{f}_k + 2 \bm{f}(\tilde{\bm{x}}_{k-1/2}, \tilde{\bm{u}}_{k-1/2}) + \bm{f}_{k-1} \right)
\end{equation}
for the TDS, for $k=1 \ldots n-1$, where $\tilde{\bm{x}}_{k-1/2} = \tilde{\bm{x}}(t_{k-1} + \Delta t/2)$ and $\tilde{\bm{u}}_{k-1/2} = \tilde{\bm{u}}(t_{k-1} + \Delta t/2)$. Notice that the ``mirroring'' operation has to be performed for node $k=1$. The average local dicretization errors throughout the interval between nodes $k-1$ and $k$ are estimated as
\begin{equation}\label{eq:error_abs}
    \bm{\varepsilon}_{abs,k} \approx \frac{1}{2}\left| \bm{\hat{x}}_k - \bm{x}_k \right| \, ,
\end{equation}
for $k=1 \ldots n$. The relative, local discretization error is also introduced as
\begin{equation}\label{eq:error_rel}
    \varepsilon_{rel,i,k} = \frac{\varepsilon_{abs,i,k}}{|\bar{x}_{i}| + 1} \, ,
\end{equation}
for $k=1 \ldots n$ and $i = 1 \ldots (2f+2m)$, where $|\bar{x}_{i}|=\frac{1}{n-1}\sum_{k=1}^{n-1} |x_{i,k}|$.


%\newpage %Just because of unusual number of tables stacked at end
%\bibliography{apssamp}% Produces the bibliography via BibTeX.


%\bibliography{literature/mydatabase_LRI}% Produces the bibliography via BibTeX.


\begin{thebibliography}{30}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
  \def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
  \def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
  \def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{\citenamefont{Erdemir et~al.}(2007)\citenamefont{Erdemir, Mc{L}ean,
  Herzog, and van~den Bogert}}]{erdemir:mclean:herzog:2007}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Erdemir}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Mc{L}ean}},
  \bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Herzog}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.~J.} \bibnamefont{van~den Bogert}}, ``\bibnamefont{Model-Based Estimation of Muscle Forces Exerted During Movements},''
  \bibinfo{journal}{Clinical Biomechanics} \textbf{\bibinfo{volume}{22}},
  \bibinfo{pages}{131} (\bibinfo{year}{2007}).

\bibitem[{\citenamefont{Crowninshield and
  Brand}(1981)}]{crowninshield:brand:1981}
\bibinfo{author}{\bibfnamefont{R.~D.} \bibnamefont{Crowninshield}}
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{R.~A.} \bibnamefont{Brand}}, ``\bibnamefont{Physiologically Based Criterion of Muscle Force Prediction in Locomotion},''
  \bibinfo{journal}{Journal of Biomechanics} \textbf{\bibinfo{volume}{14}},
  \bibinfo{pages}{793} (\bibinfo{year}{1981}).

\bibitem[{\citenamefont{Ackermann and
  Schiehlen}(2007)}]{ackermann:schiehlen:2007}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Ackermann}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Schiehlen}}, ``\bibnamefont{Physiological Methods to Solve the Force-Sharing Problem in Biomechanics},'' in
  \emph{\bibinfo{booktitle}{Proceedings of Multibody Dynamics 2007 (ECCOMAS
  Thematic Conference)}}, edited by
  \bibinfo{editor}{\bibfnamefont{C.}~\bibnamefont{Bottasso}},
  \bibinfo{editor}{\bibfnamefont{P.}~\bibnamefont{Masarati}}, \bibnamefont{and}
  \bibinfo{editor}{\bibfnamefont{L.}~\bibnamefont{Trainelli}}
  (\bibinfo{address}{Milan, Italy, 25-28 June 2007}, \bibinfo{year}{2007}).

\bibitem[{\citenamefont{Neptune et~al.}(2004)\citenamefont{Neptune, Zajac, and
  Kautz}}]{neptune:zajac:kautz:2004}
\bibinfo{author}{\bibfnamefont{R.~R.} \bibnamefont{Neptune}},
  \bibinfo{author}{\bibfnamefont{F.~E.} \bibnamefont{Zajac}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.~A.} \bibnamefont{Kautz}}, ``\bibnamefont{Muscle Force Redistributes Segmental Power for Body Progression During Walking},''
  \bibinfo{journal}{Gait and Posture} \textbf{\bibinfo{volume}{19}},
  \bibinfo{pages}{194} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Thelen et~al.}(2003)\citenamefont{Thelen, Anderson, and
  Delp}}]{thelen:anderson:delp:2003}
\bibinfo{author}{\bibfnamefont{D.~G.} \bibnamefont{Thelen}},
  \bibinfo{author}{\bibfnamefont{F.~C.} \bibnamefont{Anderson}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{S.~L.} \bibnamefont{Delp}}, ``\bibnamefont{Generating Dynamic Simulations of Movement Using Computed Muscle Control},''
  \bibinfo{journal}{Journal of Biomechanics} \textbf{\bibinfo{volume}{36}},
  \bibinfo{pages}{321} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Thelen and Anderson}(2006)}]{thelen:anderson:2006}
\bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Thelen}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.~C.} \bibnamefont{Anderson}}, ``\bibnamefont{Using Computed Muscle Control to Generate Forward Dynamic Simulations of Human Walking from Experimental Data},''
  \bibinfo{journal}{Journal of Biomechanics} \textbf{\bibinfo{volume}{39}},
  \bibinfo{pages}{1107} (\bibinfo{year}{2006}).

\bibitem[{\citenamefont{Seth and Pandy}(2007)}]{seth:pandy:2007}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Seth}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.~G.} \bibnamefont{Pandy}}, ``\bibnamefont{A Neuromusculoskeletal Tracking Method for Estimating Individual Muscle Forces in Human Movement},''
  \bibinfo{journal}{Journal of Biomechanics} \textbf{\bibinfo{volume}{40}},
  \bibinfo{pages}{356} (\bibinfo{year}{2007}).

\bibitem[{\citenamefont{Sellers et~al.}(2005)\citenamefont{Sellers, Cain, Wang,
  and Crompton}}]{sellers:cain:wang:2005}
\bibinfo{author}{\bibfnamefont{W.~I.} \bibnamefont{Sellers}},
  \bibinfo{author}{\bibfnamefont{G.~M.} \bibnamefont{Cain}},
  \bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Wang}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{R.~H.} \bibnamefont{Crompton}}, ``\bibnamefont{Stride Lengths, Speed and Energy Costs in Walking of Autralopithecus Afarensis: Using Evolutionary Robotics to Predict Locomotion of Early Human Ancestors},''
  \bibinfo{journal}{Journal of the Royal Society Interface}
  \textbf{\bibinfo{volume}{2}}, \bibinfo{pages}{431} (\bibinfo{year}{2005}).

\bibitem[{\citenamefont{Sellers and Manning}(2007)}]{sellers:manning:2007}
\bibinfo{author}{\bibfnamefont{W.~I.} \bibnamefont{Sellers}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.~L.} \bibnamefont{Manning}},  ``\bibnamefont{Estimating Dinosaur Maximum Running Speeds Using Evolutionary Robotics},''
  \bibinfo{journal}{Proceedings of the Royal Society B} pp.
  \bibinfo{pages}{1--6} (\bibinfo{year}{2007}).

\bibitem[{\citenamefont{Anderson and Pandy}(2001)}]{anderson:pandy:2001b}
\bibinfo{author}{\bibfnamefont{F.~C.} \bibnamefont{Anderson}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.~G.} \bibnamefont{Pandy}}, ``\bibnamefont{Dynamic Optimization of Human Walking},''
  \bibinfo{journal}{Journal of Biomechanical Engineering}
  \textbf{\bibinfo{volume}{123}}, \bibinfo{pages}{381} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Neptune et~al.}(2000)\citenamefont{Neptune, Wright, and
  van~den Bogert}}]{neptune:wright:bogert:2000}
\bibinfo{author}{\bibfnamefont{R.~R.} \bibnamefont{Neptune}},
  \bibinfo{author}{\bibfnamefont{I.~C.} \bibnamefont{Wright}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.~J.} \bibnamefont{van~den
  Bogert}}, ``\bibnamefont{A Method for Numerical Simulation of Single Limb Ground Contact Events: Application to Heel-Toe Running},'' \bibinfo{journal}{Computer Methods in Biomechanics and Biomedical
  Engineering} \textbf{\bibinfo{volume}{3}}, \bibinfo{pages}{321}
  (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Mc{L}ean et~al.}(2003)\citenamefont{Mc{L}ean, Su, and
  {v}an~{d}en Bogert}}]{mclean:su:bogert:2003}
\bibinfo{author}{\bibfnamefont{S.~G.} \bibnamefont{Mc{L}ean}},
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Su}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.~J.} \bibnamefont{{v}an~{d}en Bogert}}, ``\bibnamefont{Development and Validation of a 3-D Model to Predict Knee Joint Loading During Dynamic Movement},''
  \bibinfo{journal}{Journal of Biomechanical Engineering}
  \textbf{\bibinfo{volume}{125}}, \bibinfo{pages}{864} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Betts}(1998)}]{betts:1998}
\bibinfo{author}{\bibfnamefont{J.~T.} \bibnamefont{Betts}}, ``\bibnamefont{Survey of Numerical Methods for Trajectory Optimization},''
  \bibinfo{journal}{Journal of Guidance, Control, and Dynamics}
  \textbf{\bibinfo{volume}{21}}, \bibinfo{pages}{2} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Ackermann}(2007)}]{ackermann:2007}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Ackermann}}, ``\bibnamefont{A Novel Approach to Generate Physiological Human Walking Pattern},'' in
  \emph{\bibinfo{booktitle}{Proceedings of the ASME 2007, IDETC/CIE 2007 (Las
  Vegas, USA, September 4-7), paper IDETC35014}} (\bibinfo{year}{2007}).

\bibitem[{\citenamefont{Ren et~al.}(2007)\citenamefont{Ren, Jones, and
  Howard}}]{ren:jones:howard:2007}
\bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Ren}},
  \bibinfo{author}{\bibfnamefont{R.~K.} \bibnamefont{Jones}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Howard}}, ``\bibnamefont{Predictive Modelling of Human Walking over a Complete Gait Cycle},''
  \bibinfo{journal}{Journal of Biomechanics} \textbf{\bibinfo{volume}{40}},
  \bibinfo{pages}{1567} (\bibinfo{year}{2007}).

\bibitem[{\citenamefont{Kaplan}(2000)}]{kaplan:2000}
\bibinfo{author}{\bibfnamefont{M.~L.} \bibnamefont{Kaplan}}, ``\bibnamefont{Efficient Optimal Control of Large-Scale Biomechanical Systems},'' \bibinfo{type}{PhD
  thesis}, \bibinfo{school}{Stanford University} (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Kaplan and Heegaard}(2001)}]{kaplan:heegaard:2001}
\bibinfo{author}{\bibfnamefont{M.~L.} \bibnamefont{Kaplan}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.~H.} \bibnamefont{Heegaard}}, ``\bibnamefont{Predictive Algorithms for Neuromuscular Control of Human Locomotion},''
  \bibinfo{journal}{Journal of Biomechanics} \textbf{\bibinfo{volume}{34}},
  \bibinfo{pages}{1077} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Kaplan and Heegaard}(2002)}]{kaplan:heegaard:2002}
\bibinfo{author}{\bibnamefont{Kaplan}} \bibnamefont{and}
  \bibinfo{author}{\bibnamefont{Heegaard}}, ``\bibnamefont{Second-Order Optimal Control Algorithm for Complex Systems},'' \bibinfo{journal}{International
  Journal for Numerical Methods in Engineering} \textbf{\bibinfo{volume}{53}},
  \bibinfo{pages}{2043} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Delp and Loan}(2000)}]{delp:loan:2000}
\bibinfo{author}{\bibfnamefont{S.~L.} \bibnamefont{Delp}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.~P.} \bibnamefont{Loan}}, ``\bibnamefont{A Computational Framework for Simulating and Analyzing Human and Animal Movement},''
  \bibinfo{journal}{IEEE Computing in Science and Engineering}
  \textbf{\bibinfo{volume}{2}}, \bibinfo{pages}{46} (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Gerritsen et~al.}(1998)\citenamefont{Gerritsen, van~den
  Bogert, Hulliger, and Zernicke}}]{gerritsen:bogert:hulliger:1998}
\bibinfo{author}{\bibfnamefont{K.~G.~M.} \bibnamefont{Gerritsen}},
  \bibinfo{author}{\bibfnamefont{A.~J.} \bibnamefont{van~den Bogert}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Hulliger}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{R.~F.} \bibnamefont{Zernicke}}, ``\bibnamefont{Intrinsic Muscle Properties Facilitate Locomotor Control -- A Computer Simulation Study},''
  \bibinfo{journal}{Motor Control} \textbf{\bibinfo{volume}{2}},
  \bibinfo{pages}{206} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Hardin et~al.}(2004)\citenamefont{Hardin, Su, and
  van~den Bogert}}]{hardin:su:bogert:2004}
\bibinfo{author}{\bibfnamefont{E.~C.} \bibnamefont{Hardin}},
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Su}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.~J.} \bibnamefont{van~den Bogert}}, ``\bibnamefont{Foot and Ankle Forces During an Automobile Collision: the Influence of Muscles},''
  \bibinfo{journal}{Journal of Biomechanics} \textbf{\bibinfo{volume}{37}},
  \bibinfo{pages}{637} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Aerts and de~Clercq}(1993)}]{aerts:declercq:1993}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Aerts}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{de~Clercq}}, ``\bibnamefont{Deformation Characteristics of the Heel Region of the Shod Foot During a Simulated Heel Strike: The
  Effect of Varying Midsole Hardness},''
  \bibinfo{journal}{Journal of Sports Sciences} \textbf{\bibinfo{volume}{11}},
  \bibinfo{pages}{449} (\bibinfo{year}{1993}).

\bibitem[{\citenamefont{van~den Bogert et~al.}(1989)\citenamefont{van~den
  Bogert, Schamhardt, and Crowe}}]{bogert:schamhardt:crowe:1989}
\bibinfo{author}{\bibfnamefont{A.~J.} \bibnamefont{van~den Bogert}},
  \bibinfo{author}{\bibfnamefont{H.~C.} \bibnamefont{Schamhardt}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Crowe}}, ``\bibnamefont{Simulation of Quadrupedal Locomotion Using a Rigid Body Model},''
  \bibinfo{journal}{Journal of Biomechanics} \textbf{\bibinfo{volume}{22}},
  \bibinfo{pages}{33} (\bibinfo{year}{1989}).

\bibitem[{\citenamefont{Betts and Huffman}(1998)}]{betts:huffman:1998}
\bibinfo{author}{\bibfnamefont{J.~T.} \bibnamefont{Betts}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{W.~P.} \bibnamefont{Huffman}}, ``\bibnamefont{Mesh Refinement in Direct Transcription Methods for Optimal Control},''
  \bibinfo{journal}{Optimal Control Applications and Methods}
  \textbf{\bibinfo{volume}{19}}, \bibinfo{pages}{1} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Gill et~al.}(2002)\citenamefont{Gill, Murray, and
  Saunders}}]{gill:murray:saunders:2002}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Gill}},
  \bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Murray}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.~A.} \bibnamefont{Saunders}}, ``\bibnamefont{SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization},''
  \bibinfo{journal}{SIAM Journal on Optimization}
  \textbf{\bibinfo{volume}{12}}, \bibinfo{pages}{979} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Winter}(1991)}]{winter:1991}
\bibinfo{author}{\bibfnamefont{D.~A.} \bibnamefont{Winter}},
  \emph{\bibinfo{title}{The biomechanics and motor control of human gait:
  normal, elderly and pathological}} (\bibinfo{publisher}{University of
  Waterloo Press}, \bibinfo{address}{Waterloo}, \bibinfo{year}{1991}),
  \bibinfo{edition}{2nd} ed.

\bibitem[{\citenamefont{Martins et~al.}(2003)\citenamefont{Martins, Sturdza,
  and Alonso}}]{martins:sturdza:alonso:2003}
\bibinfo{author}{\bibfnamefont{J.~R.~R.} \bibnamefont{Martins}},
  \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Sturdza}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.~J.} \bibnamefont{Alonso}}, ``\bibnamefont{The Complex-Step Derivative Approximation},''
  \bibinfo{journal}{ACM Transactions on Mathematical Software}
  \textbf{\bibinfo{volume}{29}}, \bibinfo{pages}{245} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Bischof et~al.}(1992)\citenamefont{Bischof, Carle,
  Corliss, Griewank, and Hovland}}]{bischof:carle:corliss:1992}
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Bischof}},
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Carle}},
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Corliss}},
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Griewank}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Hovland}}, ``\bibnamefont{ADIFOR -- Generating Derivative Codes from Fortran Programs},''
  \bibinfo{journal}{Scientific Programming} \textbf{\bibinfo{volume}{1}},
  \bibinfo{pages}{1} (\bibinfo{year}{1992}),
  \urlprefix\url{www-unix.mcs.anl.gov/autodiff/ADIFOR}.

\bibitem[{\citenamefont{Witkin and Kass}(1988)}]{witkin:kass:1988}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Witkin}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Kass}}, ``\bibnamefont{Spacetime Constraints},'' in
  \emph{\bibinfo{booktitle}{Computer Graphics (Proc. Siggraph 1988)}}
  (\bibinfo{year}{1988}), vol.~\bibinfo{volume}{22}, pp.
  \bibinfo{pages}{159--168}.

\bibitem[{\citenamefont{Bertram and Ruina}(2001)}]{bertram:ruina:2001}
\bibinfo{author}{\bibfnamefont{J.~E.~A.} \bibnamefont{Bertram}}
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Ruina}}, ``\bibnamefont{Multiple Walking Speed-Frequency Relations are Predicted by Constrained Optimization},''
  \bibinfo{journal}{Journal of Theoretical Biology}
  \textbf{\bibinfo{volume}{209}}, \bibinfo{pages}{445} (\bibinfo{year}{2001}).

\end{thebibliography}


\end{document}
%
% ****** End of file DC_paper.tex ******
