Skip to main content

Treatment Optimization


To utilize Treatment Optimization, you will need a valid GPOPS-II license and have it downloaded and set up with MATLAB. GPOPS-II is a powerful direct collocation optimal control solver that forms the core of Treatment Optimization.


For help with GPOPS-II, click here

GPOPS-II Setting Files

GPOPS-II offers two nonlinear problem (NLP) solvers: SNOPT and IPOPT. SNOPT utilizes a quasi-Newton SQP active set method, while IPOPT employs an interior-point method. The choice of solver and settings depends on convergence and computation time considerations. We suggest using IPOPT in first derivative mode as the NLP solver without mesh refinement for the tutorials and examples provided. A fixed mesh of 60 collocation points divided into 10 intervals over the entire gait cycle was used to reduce computation time.


For more information regarding GPOPS-II settings, click here


Treatment Optimization is an advanced tool that harnesses the power of neuromusculoskeletal modeling and optimal control techniques to revolutionize the field of biomechanics and human movement analysis. By combining the capabilities of OpenSim and GPOPS-II, Treatment Optimization offers a comprehensive solution for optimizing control strategies and generating predictions of human movement.

At its core, Treatment Optimization consists of three modules: Tracking Optimization, Verification Optimization, and Design Optimization. These modules work together to gradually refine and improve the predicted movement patterns.

The Tracking Optimization module focuses on finding the optimal control strategies, either torque-driven or synergy-driven, that closely match the experimental motion, moments, external loads (if applicable), and muscle activity (if applicable). By minimizing the discrepancy between the predicted and experimental data, Tracking Optimization aims to reproduce the experimental movement.

The Verification Optimization module ensures that the calibrated controllers, derived from Tracking Optimization, can reproduce the experimental data even when tracking those quantities is eliminated. This validation step enhances the reliability and robustness of the predictions generated by Treatment Optimization.

The Design Optimization module offers users the flexibility to test and explore various treatment strategies by customizing cost functions, path constraints, and terminal conditions. Design Optimization empowers researchers, clinicians, and practitioners to tailor the treatment to the specific needs and goals of the individual, paving the way for optimized outcomes and improved patient care.

Treatment Optimization enforces skeletal dynamic constraints and ensures consistency between muscle-produced moments and inverse dynamic moments. This important aspect helps to ensure that the predicted movements are dynamically consistent and align with the physiological principles governing human motion.


Treatment Optimization plays a crucial role in enhancing our understanding of neuromusculoskeletal dynamics and improving treatment strategies. By providing researchers and practitioners with a powerful tool to optimize control strategies and predict human movement, Treatment Optimization opens doors for personalized treatment plans and enhanced athletic performance. The ability to generate predictions easily using Treatment Optimization significantly contributes to advancements in personalized and tailored treatment plans, ultimately leading to improved outcomes and an enhanced quality of life for individuals.

Problem Formulation

When solving direct collocation optimal control problems involving human movement, implicit skeletal dynamics have been found to work better than explicit skeletal dynamics. However, GPOPS-II does not directly handle implicit dynamics and requires explicit dynamics for dynamic constraints.

To address this limitation of GPOPS-II, an implicit form of skeletal dynamics is introduced in the path constraints of GPOPS-II within Treatment Optimization. This formulation requires treating joint jerk as an additional control variable. The kinematic derivative relationships, where joint jerk is the first-time derivative of joint acceleration, joint acceleration is the first-time derivative of joint velocity, and joint velocity is the first-time derivative of joint position, are utilized to establish the required explicit dynamics. A joint jerk minimization term is included in the cost function to ensure that the joint jerk controls are unique. Allowing for the enforcement of implicit skeletal dynamics in GPOPS-II, even though the solver itself does not inherently handle implicit dynamics.

Common File Settings

Design Variable Bounds Terms

The search space for the both states and controls are specified under the "RCNLDesignVariableBoundsTerms> tag. Instead of assigning the upper and lower bounds for each state (ex. joint position, joint velocity, etc.), we assign a range multiple that applies for all coordinates. The range multiple works as follows for any joint and for both states and controls.

upperLimit(statePositions)=max(q1,2,3,...n)+jointPositionsMultiplerange(q1,2,3,...n)upperLimit(statePositions)= max(q_{1,2,3,...n}) + jointPositionsMultiple * range(q_{1,2,3,...n})
lowerLimit(statePositions)=min(q1,2,3,...n)jointPositionsMultiplerange(q1,2,3,...n)lowerLimit(statePositions)= min(q_{1,2,3,...n}) - jointPositionsMultiple * range(q_{1,2,3,...n})
Synergy driven controller
upperLimit(synergyActivations)=maxSynergyActivationsupperLimit(synergyActivations)= maxSynergyActivations
lowerLimit(synergyActivations)=0;lowerLimit(synergyActivations)= 0;
upperLimit(synergyWeights)=maxSynergyWeightsupperLimit(synergyWeights)= maxSynergyWeights
lowerLimit(synergyWeights)=0;lowerLimit(synergyWeights)= 0;
Torque driven controller
upperLimit(torqueController)=max(load1,2,3,...n)+torqueControlMultiplerange(load1,2,3,...n)upperLimit(torqueController)= max(load_{1,2,3,...n}) + torqueControlMultiple * range(load_{1,2,3,...n})
lowerLimit(torqueController)=min(load1,2,3,...n)torqueControlMultiplerange(load1,2,3,...n)lowerLimit(torqueController)= min(load_{1,2,3,...n}) - torqueControlMultiple * range(load_{1,2,3,...n})