OpenSim Moco: Musculoskeletal optimal control

Musculoskeletal simulations are used in many different applications, ranging from the design of wearable robots that interact with humans to the analysis of patients with impaired movement. Here, we introduce OpenSim Moco, a software toolkit for optimizing the motion and control of musculoskeletal models built in the OpenSim modeling and simulation package. OpenSim Moco uses the direct collocation method, which is often faster and can handle more diverse problems than other methods for musculoskeletal simulation. Moco frees researchers from implementing direct collocation themselves—which typically requires extensive technical expertise—and allows them to focus on their scientific questions. The software can handle a wide range of problems that interest biomechanists, including motion tracking, motion prediction, parameter optimization, model fitting, electromyography-driven simulation, and device design. Moco is the first musculoskeletal direct collocation tool to handle kinematic constraints, which enable modeling of kinematic loops (e.g., cycling models) and complex anatomy (e.g., patellar motion). To show the abilities of Moco, we first solved for muscle activity that produced an observed walking motion while minimizing squared muscle excitations and knee joint loading. Next, we predicted how muscle weakness may cause deviations from a normal walking motion. Lastly, we predicted a squat-to-stand motion and optimized the stiffness of an assistive device placed at the knee. We designed Moco to be easy to use, customizable, and extensible, thereby accelerating the use of simulations to understand the movement of humans and other animals.

auxiliary dynamics, implicit 0 = φ(q, p) kinematic constraints y, x, λ, p) d t k = 1, . . . , K g L ≤ g(t, y, x, λ, p) ≤ g U path constraints y 0,L ≤ y 0 ≤ y 0,U y f ,L ≤ y f ≤ y f ,U initial and final states initial and final controls with respect to t 0 ∈ [t 0,L , t 0,U ] initial time The problem contains constraints for the system's multibody dynamics (involving the mass matrix M ; applied forces f app from gravity, muscles, etc.; and centripetal and Coriolis terms f inertial ) and any auxiliary dynamics. The auxiliary dynamics may be expressed as explicit ( fż ,ex ) or implicit ( fż ,im ) differential equations; the auxiliary state variables z are partitioned according to whether they are integrated using explicit differential equations (z ex ) or implicit differential equations (z im ). The system may contain position-level (holonomic) kinematic constraints to, for example, weld a foot to a bicycle pedal. Each constraint is enforced by forces exerted by tissue, bones, bodies, or other parts of the modeled system. These generalized constraint forces are applied in the constrained directions (e.g., the six degrees of freedom between the foot and pedal), and we introduce time-varying Lagrange multiplier variables λ to solve for these forces. The derivatives of the kinematic constraints φ with respect to the generalized coordinates yields the kinematic constraint Jacobian G; the transpose of this matrix converts the Lagrange multipliers into generalized forces along the system's degrees of freedom. See below for details on how Moco handles kinematic constraints.
Additionally, the problem contains K boundary constraints V k (with bounds V L,k and V U,k ) and algebraic path constraints g over the motion (with time-invariant bounds g L and g U ). The cost terms and boundary constraints may depend on initial and final time; states; controls; kinematic constraint multipliers (required for joint reactions); time-invariant parameters; and an integral, S c, j or S b,k , over the motion.

Kinematic constraints
Support for problems incorporating kinematic constraints is a key feature of Moco. Understanding how kinematic constraints are handled in multibody dynamics is valuable for understanding how Moco handles such problems. Consider a two-dimensional point mass system with coordinates q x and q y constrained to a parabola 0 = q y − q 2 x . To prevent the point mass from violating the constraint, we must apply a force perpendicular to the parabola. Each constraint has a corresponding scalar force variable, called a Lagrange multiplier λ. We must solve for the required magnitudes of these constraint forces, but the direction in which we apply each of these forces is determined by the derivative of the constraint equations. We gather the derivatives of the constraint equations in the kinematic constraint Jacobian matrix G. Each row in this matrix contains the derivative of a single constraint equation with respect to each degree of freedom, and the matrix has a column for each degree of freedom. For the parabola example, the Jacobian is (−2q x , 1). The transpose of this matrix, G T , contains columns that are vectors in the state space which are perpendicular to each constraint. For our single constraint, the vector (−2q x , 1) T is perpendicular to the parabola. The matrix-vector product between the Jacobian transpose and the Lagrange multipliers, G T λ, yields the vector of generalized forces (whose length is the number of degrees of freedom) necessary for enforcing the kinematic constraints. For our point mass example, the generalized forces yielded by the vector-scalar product (−2q x , 1) T λ keep the point mass on the parabola. To apply these forces to the multibody system, we include the G T λ term in the multibody dynamics equations of motion.
When simulating a multibody system with time-stepping forward integration, we first ensure the initial generalized coordinates and speeds satisfy the kinematic constraints φ(q) = 0 and their first derivativeφ(q, u) = 0 via a root-solve. During the integration, we solve for generalized accelerations and Lagrange multipliers that obey the multibody dynamics equations of motion and the second derivative of the kinematic constraints,φ(q, u,u) = 0. Numerically integrating the resulting generalized accelerations yields generalized coordinates and speeds that approximately lie on the constraint manifold defined by φ(q) = 0; to fix any errors in the constraints caused by numerical integration error, we project the generalized coordinates and speeds back onto the constraint manifold [1].
In direct collocation, we solve for the entire trajectory of the system-including the generalized coordinates, generalized speeds, and Lagrange multipliers-all at once. When expressing multibody dynamics as implicit differential equations, the generalized accelerations are also unknowns. How we solve for the trajectory of the system in the presence of kinematic constraints depends on the transcription scheme; see the remainder of this Appendix for details.

Moco's optimal control problem with prescribed kinematics
A common task in musculoskeletal biomechanics is to estimate the muscle and actuator behavior that drove an observed motion. We can solve this problem by minimizing the error between the observed motion and the simulated motion, as with Computed Muscle Control (using the "slow target") [2] or MocoTrack. Alternatively, we can prescribe the motion exactly, as with Static Optimization [3], electromyography-driven simulation [4], and the Muscle Redundancy Solver [5]. Consider a two-dimensional point mass with coordinates q x and q y for which we prescribe a circular motion via the functionsq x (t) = cos(t) andq y (t) = sin(t). We can either add these functions to the kinematic constraints φ(q), or we can substitute these functions into the equations of motion, thereby eliminating the variables q x and q y . With Moco, users can choose either the former approach through OpenSim's Coordinate, or the latter (and usually preferable) approach using PositionMotion, a new component that employs Simbody's Motion class. Prescribing kinematics by eliminating variables leads to a problem that is robust and fast-the nonlinear multibody dynamics are removed from 3/11 the optimization problem-but prevents predicting kinematic deviations from the observed motion.
When we prescribe kinematics in Moco by eliminating variables, we replace the problem in Eq 1 with the following: auxiliary states We replace the kinematic variables q and u with known quantitiesq andû. The system still depends on auxiliary state variables z and control variables x, and includes auxiliary dynamics. If none of the parameter variables affect the multibody system, then the multibody dynamics are reduced to a force balance: muscles and other force elements must generate the net generalized forces determined by the kinematics and external loads data. Whether the motion is prescribed by adding constraints or eliminating variables, OpenSim supplements the modeled force elements with Lagrange multipliers to ensure the prescribed motion is achieved. When using PositionMotion with Moco, we require that the prescribed motion's Lagrange multipliers are zero, thereby ensuring the motion is fully generated by the modeled force elements. The easiest way to prescribe kinematics in Moco is to use the MocoInverse tool, which uses PositionMotion internally.

Transcription schemes Trapezoidal transcription
The trapezoidal scheme transcribes the optimal control problem into a nonlinear program by approximating integrals using the trapezoidal rule. As a second-order scheme, trapezoidal transcription exhibits accuracy that improves four-fold when halving the mesh interval (i.e., time step).
We discretize the continuous variables t, y, x, and λ on a mesh of time points t i defined 4/11 by dimensionless time τ i , yielding n mesh intervals with durations h i : For conciseness, we define the following function: where trap i (F (η, p))) is a trapezoidal rule approximation of the area under the function F for mesh interval i, and η represents any subset of continuous variables. We define the kinematic and dynamic explicit differential equations as: The mass matrix M , the centripetal and Coriolis forces f inertial , and kinematic constraint Jacobian G are computed by Simbody (OpenSim's multibody dynamics engine) using order-N recursive algorithms 1 [1]. The applied forces f app can include any force elements supported by Simbody, and are often defined by OpenSim components that provide access to existing Simbody force elements (e.g., SmoothSphereHalfSpaceForce) or that define custom Simbody force elements (e.g., DeGrooteFregly2016Muscle). Simbody computes the constraint Jacobian G based on any Simbody kinematic constraints that the OpenSim model adds to the system. In Moco, the Lagrange multipliers λ are explicit optimization variables; this approach differs from that in time-stepping forward integrations, in which Simbody solves for generalized accelerations and Lagrange multipliers simultaneously. The result of the trapezoidal transcription, with multibody dynamics expressed as explicit differential equations, is the following nonlinear program: fu(t, y, x, λ, p) In this form, the problem can be solved directly by a nonlinear program solver. We introduce the algebraic (control) variable ζ as the derivative of auxiliary state variables whose dynamics are expressed with an implicit differential equation. When expressing the multibody dynamics implicitly, we remove the constraint involving fu, introduce generalized accelerations as an algebraic variable υ and enforce multibody dynamics in "inverse dynamics" form: The constant υ B is a large positive number (1000 by default). The dynamic, kinematic, and path constraints are enforced at a set of discrete time points, so the quadratic spline approximation to the continuous variables may violate the original continuous-time constraints between the discrete time points. For this reason, a mesh with more points leads to a more accurate solution.
Our implementation of trapezoidal transcription handles kinematic constraints, but not in the most robust fashion. We enforce φ but not its time derivatives; enforcing the constraints at only the position level yields an index-3 system of differential-algebraic equations, which are challenging to solve [6][7][8] (the differential-algebraic equations are "index-3" because the algebraic constraints φ must be differentiated three times to convert the system of differential-algebraic equations into ordinary differential equations). To improve numerical conditioning, we minimize the Lagrange multipliers associated with the constraints (with weight w λ ; see Eq 7).
When kinematics are prescribed (see "Moco's optimal control problem with prescribed kinematics"), multibody dynamics must be expressed implicitly and kinematic constraints are not enforced; we expect the prescribed kinematics to already obey the constraints.

Hermite-Simpson transcription
The Hermite-Simpson scheme transcribes the optimal control problem into a nonlinear program by approximating integrals using a Hermite interpolant and Simpson integration rule. As a third-order scheme, Hermite-Simpson transcription exhibits accuracy that improves eight-fold when halving the mesh interval.
We use a similar dimensionless time mesh as for the trapezoidal scheme, with n mesh intervals with durations h i . We also introduce collocation points at the midpoints of the mesh intervals, leading to a total of 2n + 1 time points at which we discretize the continuous variables: whereτ i denote mesh interval midpoints. For conciseness, we define the following functions: where hermite i () represents the Hermite interpolant, which yields the value of the midpoint of a segment in a Hermite spline, and simpson i () represents the Simpson integration rule. Again, F is a function for mesh interval i, and η represents any subset of continuous variables. The symbolsη i represent continuous variables evaluated at the mesh interval midpoints. Using the explicit multibody dynamics function fu defined previously, Hermite-Simpson 7/11 transcription results in the following nonlinear program: j (t, y, x, λ, p)) fq(q, u, γ, p)) i = 1, . . . , n q i = q i−1 + simpson i ( fq(q, u, γ, p)) i = 1, . . . , n fu(t, y, x, λ, p)) i = 1, . . . , n fu(t, y, x, λ, p)) i = 1, . . . , n The variablesȳ i (includingq i ,ū i ,z ex,i ,z im,i ),x i ,ζ i ,λ i , andγ i are associated with the midpoint of mesh interval i. The midpoint states are included as explicit optimization variablesȳ i since we are using the separated form of Hermite-Simpson transcription [8] Eq 12 includes the first and second time derivatives of the kinematic constraints: 0 =φ(t, y, x, λ, p) = G(q, p)u +Ġ(q, p)u = G(q, p) fu(t, y, x, λ, p) +Ġ(q, p)u.
Including these constraints reduces the index of the differential-algebraic equations from 3 to 1, which provides numerical benefits [6][7][8]. However, simply appending the kinematic constraints to the unconstrained optimal control problem would overconstrain the resulting non-linear program, so we use the method for handling kinematic constraints that was introduced by Posa and colleagues [9]. In addition to the Lagrange multiplier variables we introduced for trapezoidal transcription, we introduce velocity correction variables, γ, whose multiplication with the transpose of the kinematic constraint Jacobian, G(q, p) T γ, yields a correction to the generalized speeds that is perpendicular to the constraint manifold φ = 0. We update the function fq used in Eq 12 accordingly: Without this correction, the Hermite splines cannot simultaneously satisfy both the differential equations at the mesh interval midpoints and the kinematic constraints at the mesh points.
Since velocity-level kinematic constraints (Eq 13) are already enforced at the mesh points, this correction term only appears at the mesh interval midpoints (i.e., G(q i , p) Tγ i ); therefore, velocity corrections at the mesh points are zero and do not appear in Eq 12. Although Simbody supports non-holonomic and acceleration-level constraints, Moco only supports holonomic (position-level) constraints.
For implicit multibody dynamics, we remove the constraints involving fu and introduce generalized accelerations as algebraic variables υ to enforce multibody dynamics in "inverse dynamics" form:

Automated test suite
Moco contains an automated test suite with over 80 test cases. The test suite contains four major types of tests: analytical tests, interface tests, algorithmic tests, and regression tests.
In analytical tests, we ensure Moco produces the correct solutions for problems with known analytical solutions. These tests are commonly used to verify optimal control software and for performance benchmarking [15]. We solve two such problems: the "linear tangent steering" problem described in the main article and a boundary value problem involving a set of linear first-order differential equations presented by Kirk et al. [10] (Example 5.1-1, pg. 198). In addition, we include analytical tests using the SmoothSphereHalfSpaceForce contact model to ensure that the smooth model produces correct forces. For example, we test that a point mass containing a contact element with dissipation comes to rest with a normal force that matches the weight of the system after being dropped from a height.
Interface tests ensure that a change made by a user though the Moco interface is reflected within the software. Any values set by the user (using a "set" method) should be returned by Moco when queried (using a "get" method Algorithmic tests check that Moco solves direct collocation problems correctly. These tests check that problems using explicit and implicit dynamics produce the same trajectory, minimizing a cost term produces the expected behavior (e.g., minimum time problems hit the lower bound on final time), and that kinematic constraints are obeyed. For example, we check that adding a coordinate coupler constraint to a double pendulum problem produces coordinate trajectories that obey a linear relationship defined by the constraint. Algorithmic tests also encompass "self-consistency" tests. For example, using a model of a point mass hanging from a muscle, we check that tracking a reference trajectory obtained by a minimum time predictive problem produces the same controls as the original problem.
For complex problems for which there is no known exact solution, we use regression tests, which check that solutions to direct collocation problems do not regress from reference trajectories. These checks ensure that changes to the codebase do not produce any unintended effects. For example, we test that MocoInverse and MocoTrack always produce the same muscle activity given the same model and solver settings.

Moco's direct collocation solvers
Moco provides two solvers as subclasses of MocoSolver: MocoCasADiSolver uses the third-party CasADi library [11], and MocoTropterSolver uses a direct collocation solver we developed named Tropter. CasADi is an open-source package for algorithmic differentiation and is a bridge to nonlinear program solvers IPOPT [12], SNOPT [13], and others.
Gradient-based nonlinear program solvers require the gradient of the objective, the Jacobian of the constraints, and sometimes the Hessian of the Lagrangian [8]. To maximize computational efficiency, these derivatives are ideally computed exactly through either analytic expressions or algorithmic differentiation [11,14]. OpenSim's main distribution does not provide exact derivatives, so we use finite differences. CasADi is an ideal library for employing direct collocation, but two limitations led us to create Tropter: CasADi did not initially support finite differences, and CasADi's open-source license is more restrictive than OpenSim's. More recent versions of CasADi support finite differences and CasADi understands the structure of the nonlinear program objective and constraint functions, allowing for potentially more efficient finite difference calculations than with Tropter, which treats the nonlinear program objective and constraints as black-box functions [15]. If OpenSim provides exact derivatives in the future, we can exploit the algorithmic differentiation modes in Tropter and CasADi [16]. Those distributing Moco as a dependency of closed-source software may prefer distributing Moco without CasADi, as CasADi's "weak copyleft" GNU Lesser General Public License 3.0 places requirements on how CasADi is redistributed.