Algorithmic differentiation improves the computational efficiency of OpenSim-based optimal control simulations of movement

Algorithmic differentiation (AD) is an alternative to finite differences (FD) for evaluating function derivatives. The primarily aim of this study was to demonstrate the computational benefits of using AD instead of FD in OpenSim-based optimal control simulations. The secondary aim was to evaluate computational choices including different AD tools, different linear solvers, and the use of first- or second-order derivatives. First, we enabled the use of AD in OpenSim through a custom source code transformation tool and through the operator overloading tool ADOL-C. Second, we developed an interface between OpenSim and CasADi to perform optimal control simulations. Third, we evaluated computational choices through simulations of perturbed balance, two-dimensional predictive simulations of walking, and three-dimensional tracking simulations of walking. We performed all simulations using direct collocation and implicit differential equations. Using AD through our custom tool was between 1.8 ± 0.1 and 17.8 ± 4.9 times faster than using FD, and between 3.6 ± 0.3 and 12.3 ± 1.3 times faster than using AD through ADOL-C. The linear solver efficiency was problem-dependent and no solver was consistently more efficient. Using second-order derivatives was more efficient for balance simulations but less efficient for walking simulations. The walking simulations were physiologically realistic. These results highlight how the use of AD drastically decreases computational time of optimal control simulations as compared to more common FD. Overall, combining AD with direct collocation and implicit differential equations decreases the computational burden of optimal control simulations, which will facilitate their use for biomechanical applications.


47
Combining musculoskeletal modeling and dynamic simulation is a powerful approach to study 48 the mechanisms underlying human movements. There exist two types of model-based dynamic 49 simulations. Inverse simulations identify controls (e.g., muscle excitations) underlying a given 50 movement whereas forward simulations generate the movement based on the controls. Forward 51 simulations are commonly combined with numerical optimization techniques to compute optimal 52 controls. In the resultant optimal control problems, an objective function is minimized subject to 53 constraints describing the musculoskeletal dynamics. Researchers

62
Over the last decade, the increase in computer performance and the use of efficient numerical 63 methods have equipped researchers with better tools for generating computationally efficient optimal 64 control simulations. In particular, direct collocation methods [3,5,7-10] and implicit formulations of many additional features (e.g., code generation) and also interfaces with NLP solvers designed to 123 handle large and sparse NLPs (e.g., IPOPT). CasADi provides a high-level, symbolic, way to construct an 124 expression graph, on which SCT is applied. The resultant expression graph can be code-generated to 125 achieve the computational efficiency of pure SCT. 126 The contribution of this study is threefold. First, we enabled the use of AD in OpenSim and 127 Simbody (henceforth referred to as OpenSim). We compared two approaches: we incorporated the based on second-order derivative information) or approximated with a limited-memory quasi-Newton 139 method (L-BFGS) that only requires first-order derivative information. We found that using AD through 140 Recorder was more efficient than using FD or AD through ADOL-C whereas the efficiency of the linear 141 solver and Hessian calculation scheme was problem-dependent.

143
Tools to enable the use of AD in OpenSim active variables, which are variables that may be considered as differentiable quantities at some time about their partial derivatives, ADOL-C introduced the augmented scalar type adouble whose real part

184
We made utilizable an OpenSim function within the CasADi environment by compiling the 185 function and its derivatives as a Dynamic-link Library (DLL) that then is imported as an external function 186 for use by CasADi (Fig 2). The function derivatives can be computed through ADOL-C (AD-ADOLC in Fig   187 2) or through Recorder (AD-Recorder in Fig 2).

188
Optimal control simulations to evaluate computational choices 189 We designed three example optimal control simulations to evaluate different computational 190 choices (see Table 1 for detailed formulations). The simulations are optimal control problems whose 191 general formulation consists of computing the controls , states , and time-independent ( ) ( )

192
parameters minimizing an objective functional: (2) and to algebraic path constraints: 196 which are equality constraints if . The optimization variables are typically bounded as =

197
follows: , 198 199  Bounds Dynamic constraints  In the first example, we perturbed the balance of 9 inverted pendulums, with between 2 and 223 10 degrees of freedom (DOFs), by applying a backward translation to their base of support. The optimal 224 control problem identified the joint torques necessary to restore the pendulums' upright posture 225 within one second while minimizing the actuator effort (i.e., squared joint torques) and satisfying the 226 pendulum dynamics.

227
In the second example, we performed predictive simulations of walking with a two-228 dimensional (2D) musculoskeletal model (10 DOFs, 18 muscles actuating the lower limbs, one ideal pelvis divided by gait cycle duration).

253
We transcribed each optimal control problem into a nonlinear programming problem (NLP) 254 using a third order Radau quadrature collocation scheme. We formulated each problem in MATLAB 255 using CasADi and IPOPT. We imposed an NLP relative error tolerance of 1 x 10 -6 and used an adaptive 256 barrier parameter update strategy. We selected a number of mesh intervals for each problem such 257 that the results were qualitatively similar when using a mesh twice as fine. We used 10 and 3 initial guesses for the pendulum and walking simulations, respectively. We ran all simulations on a single core 259 of a standard laptop computer with a 2.9 GHz Intel Core i7 processor.

260
Results analysis 261 We compared CPU time and number of iterations required to solve the problems using the 262 different computational choices. First, we compared AD, using the Recorder approach, with FD.

271
Using AD-Recorder was computationally more efficient than using FD or AD-ADOLC (Fig 3). The

288
The results were obtained using mumps and an approximated Hessian.

290
The solvers from the HSL collection were on average more efficient (faster with a similar 291 number of iterations) than mumps for the pendulum simulations, but the efficiency varied for the 2D 292 predictive and 3D tracking simulations (  The comparisons are expressed as ratios (mean ± one standard deviation; results obtained with 307 solver from the HSL collection over results obtained with mumps). The ratios are averaged over results 308 from different initial guesses. Ratios larger than one indicate faster convergence, fewer iterations, or 309 more time per iteration with mumps. The use of the solvers ma57 and ma97 led to memory issues for 310 the 3D tracking simulations and these cases were therefore excluded from the analysis. The 311 simulations were run using AD-Recorder and an approximated Hessian. 312 313 Using an exact Hessian was more efficient than using an approximated Hessian for the 314 pendulum simulations but not for the 2D predictive simulations (Fig 4). The exact Hessian required less

339
The 2D predictive simulations reproduced salient features of human gait but deviated from 340 experimental data in three noticeable ways (Fig 5; S2 Movie). First, the predicted knee flexion during 341 mid-stance was limited, resulting in small knee torques. Second, the simulations produced less ankle 342 plantarflexion at push-off. Third, the vertical ground reaction forces exhibited a large peak at impact.

343
The simulations converged in less than one CPU minute (average over solutions starting from three 344 initial guesses: 36 ± 17 s and 247 ± 143 iterations; results obtained with AD-Recorder, mumps, and an 345 approximated Hessian).

434
In this study, we enabled the use of AD when performing OpenSim-based optimal control 435 simulations. We showed that using AD drastically improved the computationally efficiency of such 436 simulations. This improved efficiency is highly desirable for researchers using complex models or 437 aiming to implement such models in clinical practice where time constraints are typically more 438 stringent than in research context. Overall, the combination of AD with other efficient numerical tools 439 such as direct collocation and implicit differential equations allows overcoming the computational 440 roadblocks that have long limited the use of optimal control simulations for biomechanical 441 applications. In the future, we aim to exploit this computational efficiency to design optimal 442 treatments for neuro-musculoskeletal disorders, such as cerebral palsy.