API  4.5
For C++ developers
exampleMocoInverse.py

This is an example using the MocoInverse tool with a complex model to prescribe walking. This example also shows how to track electromyography data.

1 # -------------------------------------------------------------------------- #
2 # OpenSim Moco: exampleMocoInverse.py #
3 # -------------------------------------------------------------------------- #
4 # Copyright (c) 2023 Stanford University and the Authors #
5 # #
6 # Author(s): Christopher Dembia, Nicholas Bianco #
7 # #
8 # Licensed under the Apache License, Version 2.0 (the "License") you may #
9 # not use this file except in compliance with the License. You may obtain a #
10 # copy of the License at http://www.apache.org/licenses/LICENSE-2.0 #
11 # #
12 # Unless required by applicable law or agreed to in writing, software #
13 # distributed under the License is distributed on an "AS IS" BASIS, #
14 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. #
15 # See the License for the specific language governing permissions and #
16 # limitations under the License. #
17 # -------------------------------------------------------------------------- #
18 
19 # This example shows how to use the MocoInverse tool to exactly prescribe a
20 # motion and estimate muscle behavior for walking. The first example does not
21 # rely on electromyography data, while the second example penalizes deviation
22 # from electromyography data for a subset of muscles.
23 #
24 # Both examples use the Python utility osim.report to automatically generate a
25 # PDF that includes the trajectories of all states and controls in the solution.
26 # This utility requires a Python environment with Matplotlib and NumPy installed.
27 #
28 # See the README.txt next to this file for more information.
29 
30 import opensim as osim
31 
32 def solveMocoInverse():
33 
34  # Construct the MocoInverse tool.
35  inverse = osim.MocoInverse()
36 
37  # Construct a ModelProcessor and set it on the tool. The default
38  # muscles in the model are replaced with optimization-friendly
39  # DeGrooteFregly2016Muscles, and adjustments are made to the default muscle
40  # parameters.
41  modelProcessor = osim.ModelProcessor('subject_walk_scaled.osim')
42  modelProcessor.append(osim.ModOpAddExternalLoads('grf_walk.xml'))
43  modelProcessor.append(osim.ModOpIgnoreTendonCompliance())
44  modelProcessor.append(osim.ModOpReplaceMusclesWithDeGrooteFregly2016())
45  # Only valid for DeGrooteFregly2016Muscles.
46  modelProcessor.append(osim.ModOpIgnorePassiveFiberForcesDGF())
47  # Only valid for DeGrooteFregly2016Muscles.
48  modelProcessor.append(osim.ModOpScaleActiveFiberForceCurveWidthDGF(1.5))
49  # Use a function-based representation for the muscle paths. This is
50  # recommended to speed up convergence, but if you would like to use
51  # the original GeometryPath muscle wrapping instead, simply comment out
52  # this line. To learn how to create a set of function-based paths for
53  # your model, see the example 'examplePolynomialPathFitter.py'.
54  modelProcessor.append(osim.ModOpReplacePathsWithFunctionBasedPaths(
55  "subject_walk_scaled_FunctionBasedPathSet.xml"))
56  modelProcessor.append(osim.ModOpAddReserves(1.0))
57  inverse.setModel(modelProcessor)
58 
59  # Construct a TableProcessor of the coordinate data and pass it to the
60  # inverse tool. TableProcessors can be used in the same way as
61  # ModelProcessors by appending TableOperators to modify the base table.
62  # A TableProcessor with no operators, as we have here, simply returns the
63  # base table.
64  inverse.setKinematics(osim.TableProcessor('coordinates.sto'))
65 
66  # Initial time, final time, and mesh interval.
67  inverse.set_initial_time(0.48)
68  inverse.set_final_time(1.61)
69  inverse.set_mesh_interval(0.02)
70 
71  # By default, Moco gives an error if the kinematics contains extra columns.
72  # Here, we tell Moco to allow (and ignore) those extra columns.
73  inverse.set_kinematics_allow_extra_columns(True)
74 
75  # Solve the problem and write the solution to a Storage file.
76  solution = inverse.solve()
77  solution.getMocoSolution().write('example3DWalking_MocoInverse_solution.sto')
78 
79  # Generate a PDF with plots for the solution trajectory.
80  model = modelProcessor.process()
81  report = osim.report.Report(model,
82  'example3DWalking_MocoInverse_solution.sto',
83  bilateral=True)
84  # The PDF is saved to the working directory.
85  report.generate()
86 
87 def solveMocoInverseWithEMG():
88 
89  # This initial block of code is identical to the code above.
90  inverse = osim.MocoInverse()
91  modelProcessor = osim.ModelProcessor('subject_walk_scaled.osim')
92  modelProcessor.append(osim.ModOpAddExternalLoads('grf_walk.xml'))
93  modelProcessor.append(osim.ModOpIgnoreTendonCompliance())
94  modelProcessor.append(osim.ModOpReplaceMusclesWithDeGrooteFregly2016())
95  modelProcessor.append(osim.ModOpIgnorePassiveFiberForcesDGF())
96  modelProcessor.append(osim.ModOpScaleActiveFiberForceCurveWidthDGF(1.5))
97  modelProcessor.append(osim.ModOpReplacePathsWithFunctionBasedPaths(
98  "subject_walk_scaled_FunctionBasedPathSet.xml"))
99  modelProcessor.append(osim.ModOpAddReserves(1.0))
100  inverse.setModel(modelProcessor)
101  inverse.setKinematics(osim.TableProcessor('coordinates.sto'))
102  inverse.set_initial_time(0.48)
103  inverse.set_final_time(1.61)
104  inverse.set_mesh_interval(0.02)
105  inverse.set_kinematics_allow_extra_columns(True)
106 
107  study = inverse.initialize()
108  problem = study.updProblem()
109 
110  # Add electromyography tracking.
111  emgTracking = osim.MocoControlTrackingGoal('emg_tracking')
112  emgTracking.setWeight(50.0)
113  # Each column in electromyography.sto is normalized so the maximum value in
114  # each column is 1.0.
115  controlsRef = osim.TimeSeriesTable('electromyography.sto')
116  # Scale the tracked muscle activity based on peak levels from
117  # "Gait Analysis: Normal and Pathological Function" by
118  # Perry and Burnfield, 2010 (digitized by Carmichael Ong).
119  soleus = controlsRef.updDependentColumn('soleus')
120  gasmed = controlsRef.updDependentColumn('gastrocnemius')
121  tibant = controlsRef.updDependentColumn('tibialis_anterior')
122  for t in range(0, controlsRef.getNumRows()):
123  soleus[t] = 0.77 * soleus[t]
124  gasmed[t] = 0.87 * gasmed[t]
125  tibant[t] = 0.37 * tibant[t]
126  emgTracking.setReference(osim.TableProcessor(controlsRef))
127  # Associate actuators in the model with columns in electromyography.sto.
128  emgTracking.setReferenceLabel('/forceset/soleus_r', 'soleus')
129  emgTracking.setReferenceLabel('/forceset/gasmed_r', 'gastrocnemius')
130  emgTracking.setReferenceLabel('/forceset/gaslat_r', 'gastrocnemius')
131  emgTracking.setReferenceLabel('/forceset/tibant_r', 'tibialis_anterior')
132  problem.addGoal(emgTracking)
133 
134  # Solve the problem and write the solution to a Storage file.
135  solution = study.solve()
136  solution.write('example3DWalking_MocoInverseWithEMG_solution.sto')
137 
138  # Write the reference data in a way that's easy to compare to the solution.
139  controlsRef.removeColumn('medial_hamstrings')
140  controlsRef.removeColumn('biceps_femoris')
141  controlsRef.removeColumn('vastus_lateralis')
142  controlsRef.removeColumn('vastus_medius')
143  controlsRef.removeColumn('rectus_femoris')
144  controlsRef.removeColumn('gluteus_maximus')
145  controlsRef.removeColumn('gluteus_medius')
146  controlsRef.setColumnLabels(['/forceset/soleus_r', '/forceset/gasmed_r',
147  '/forceset/tibant_r'])
148  controlsRef.appendColumn('/forceset/gaslat_r', gasmed)
149  osim.STOFileAdapter.write(controlsRef, 'controls_reference.sto')
150 
151  # Generate a report comparing MocoInverse solutions without and with EMG
152  # tracking.
153  model = modelProcessor.process()
154  output = 'example3DWalking_MocoInverseWithEMG_report.pdf'
155  ref_files = [
156  'controls_reference.sto',
157  'example3DWalking_MocoInverseWithEMG_solution.sto']
158  report = osim.report.Report(model,
159  'example3DWalking_MocoInverse_solution.sto',
160  output=output, bilateral=True,
161  ref_files=ref_files,
162  colors=['black', 'blue', 'red'])
163  # The PDF is saved to the working directory.
164  report.generate()
165 
166 
167 # Solve the basic muscle redundancy problem with MocoInverse.
168 solveMocoInverse()
169 
170 # This problem penalizes the deviation from electromyography data for a
171 # subset of muscles.
172 solveMocoInverseWithEMG()