API  4.2
For MATLAB, Python, Java, and C++ users
exampleHangingMuscle.py

This example optimizes the motion of a point mass actuated by a single muscle and shows you how to use the AnalyzeTool and a metabolic probe with a MocoSolution.

1 # -------------------------------------------------------------------------- #
2 # OpenSim Moco: exampleHangingMuscle.py #
3 # -------------------------------------------------------------------------- #
4 # Copyright (c) 2020 Stanford University and the Authors #
5 # #
6 # Author(s): Christopher Dembia #
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 # This example includes a point mass hanging by a muscle (+x is downward),
19 # and shows how to use MocoStudy with a model that includes a muscle.
20 # Additionally, this example shows how to use OpenSim's Analyses with a
21 # MocoSolution.
22 # The trajectory optimization problem is to lift the point mass by a small
23 # distance in minimum time.
24 
25 import os
26 import opensim as osim
27 
28 def createHangingMuscleModel(ignore_activation_dynamics,
29  ignore_tendon_compliance):
30  model = osim.Model()
31  model.setName("hanging_muscle")
32  model.set_gravity(osim.Vec3(9.81, 0, 0))
33  body = osim.Body("body", 0.5, osim.Vec3(0), osim.Inertia(1))
34  model.addComponent(body)
35 
36  # Allows translation along x.
37  joint = osim.SliderJoint("joint", model.getGround(), body)
38  coord = joint.updCoordinate()
39  coord.setName("height")
40  model.addComponent(joint)
41 
42  # The point mass is supported by a muscle.
43  # The DeGrooteFregly2016Muscle is the only muscle model in OpenSim that
44  # has been tested with Moco.
45  actu = osim.DeGrooteFregly2016Muscle()
46  actu.setName("muscle")
47  actu.set_max_isometric_force(30.0)
48  actu.set_optimal_fiber_length(0.10)
49  actu.set_tendon_slack_length(0.05)
50  actu.set_tendon_strain_at_one_norm_force(0.10)
51  actu.set_ignore_activation_dynamics(ignore_activation_dynamics)
52  actu.set_ignore_tendon_compliance(ignore_tendon_compliance)
53  actu.set_fiber_damping(0.01)
54  # The DeGrooteFregly2016Muscle is the only muscle model in OpenSim that
55  # can express its tendon compliance dynamics using an implicit
56  # differential equation.
57  actu.set_tendon_compliance_dynamics_mode("implicit")
58  actu.set_max_contraction_velocity(10)
59  actu.set_pennation_angle_at_optimal(0.10)
60  actu.addNewPathPoint("origin", model.updGround(), osim.Vec3(0))
61  actu.addNewPathPoint("insertion", body, osim.Vec3(0))
62  model.addForce(actu)
63 
64  # Add metabolics probes: one for the total metabolic rate,
65  # and one for each term in the metabolics model.
66  probe = osim.Umberger2010MuscleMetabolicsProbe()
67  probe.setName("metabolics")
68  probe.addMuscle("muscle", 0.5)
69  model.addProbe(probe)
70 
71  probe = osim.Umberger2010MuscleMetabolicsProbe()
72  probe.setName("activation_maintenance_rate")
73  probe.set_activation_maintenance_rate_on(True)
74  probe.set_shortening_rate_on(False)
75  probe.set_basal_rate_on(False)
76  probe.set_mechanical_work_rate_on(False)
77  probe.addMuscle("muscle", 0.5)
78  model.addProbe(probe)
79 
80  probe = osim.Umberger2010MuscleMetabolicsProbe()
81  probe.setName("shortening_rate")
82  probe.set_activation_maintenance_rate_on(False)
83  probe.set_shortening_rate_on(True)
84  probe.set_basal_rate_on(False)
85  probe.set_mechanical_work_rate_on(False)
86  probe.addMuscle("muscle", 0.5)
87  model.addProbe(probe);
88 
89  probe = osim.Umberger2010MuscleMetabolicsProbe()
90  probe.setName("basal_rate")
91  probe.set_activation_maintenance_rate_on(False)
92  probe.set_shortening_rate_on(False)
93  probe.set_basal_rate_on(True)
94  probe.set_mechanical_work_rate_on(False)
95  probe.addMuscle("muscle", 0.5)
96  model.addProbe(probe)
97 
98  probe = osim.Umberger2010MuscleMetabolicsProbe()
99  probe.setName("mechanical_work_rate")
100  probe.set_activation_maintenance_rate_on(False)
101  probe.set_shortening_rate_on(False)
102  probe.set_basal_rate_on(False)
103  probe.set_mechanical_work_rate_on(True)
104  probe.addMuscle("muscle", 0.5)
105  model.addProbe(probe)
106 
107  body.attachGeometry(osim.Sphere(0.05))
108 
109  model.finalizeConnections()
110 
111  return model
112 
113 ignore_activation_dynamics = False
114 ignore_tendon_compliance = False
115 model = createHangingMuscleModel(ignore_activation_dynamics,
116  ignore_tendon_compliance)
117 model.printToXML("hanging_muscle.osim")
118 
119 study = osim.MocoStudy()
120 problem = study.updProblem()
121 problem.setModelAsCopy(model)
122 problem.setTimeBounds(0, [0.05, 1.0])
123 problem.setStateInfo("/joint/height/value", [0.14, 0.16], 0.15, 0.14)
124 problem.setStateInfo("/joint/height/speed", [-1, 1], 0, 0)
125 problem.setControlInfo("/forceset/muscle", [0.01, 1])
126 
127 # Initial state constraints/costs.
128 if not ignore_activation_dynamics:
129  initial_activation = osim.MocoInitialActivationGoal()
130  problem.addGoal(initial_activation)
131  initial_activation.setName("initial_activation")
132 if not ignore_tendon_compliance:
133  initial_equilibrium = osim.MocoInitialVelocityEquilibriumDGFGoal()
134  problem.addGoal(initial_equilibrium)
135  initial_equilibrium.setName("initial_velocity_equilibrium")
136  # The problem converges in fewer iterations when this goal is in cost mode.
137  initial_equilibrium.setMode("cost")
138  initial_equilibrium.setWeight(0.001)
139 
140 problem.addGoal(osim.MocoFinalTimeGoal())
141 
142 solver = study.initCasADiSolver()
143 solver.set_num_mesh_intervals(25)
144 solver.set_multibody_dynamics_mode("implicit")
145 solver.set_optim_convergence_tolerance(1e-4)
146 solver.set_optim_constraint_tolerance(1e-4)
147 
148 solution = study.solve()
149 osim.STOFileAdapter.write(solution.exportToStatesTable(),
150  "exampleHangingMuscle_states.sto")
151 osim.STOFileAdapter.write(solution.exportToControlsTable(),
152  "exampleHangingMuscle_controls.sto")
153 
154 # Conduct an analysis using MuscleAnalysis and ProbeReporter.
155 # Create an AnalyzeTool setup file.
156 analyze = osim.AnalyzeTool()
157 analyze.setName("analyze")
158 analyze.setModelFilename("hanging_muscle.osim")
159 analyze.setStatesFileName("exampleHangingMuscle_states.sto")
160 analyze.updAnalysisSet().cloneAndAppend(osim.MuscleAnalysis())
161 analyze.updAnalysisSet().cloneAndAppend(osim.ProbeReporter())
162 analyze.updControllerSet().cloneAndAppend(
163  osim.PrescribedController("exampleHangingMuscle_controls.sto"))
164 analyze.printToXML("exampleHangingMuscle_AnalyzeTool_setup.xml")
165 # Run the analysis.
166 analyze = osim.AnalyzeTool("exampleHangingMuscle_AnalyzeTool_setup.xml")
167 analyze.run()
168 
169 table_force = osim.TimeSeriesTable(
170  "analyze_MuscleAnalysis_ActiveFiberForce.sto")
171 table_velocity = osim.TimeSeriesTable(
172  "analyze_MuscleAnalysis_FiberVelocity.sto")
173 time = table_force.getIndependentColumn()
174 force = table_force.getDependentColumn("muscle").to_numpy()
175 velocity = table_velocity.getDependentColumn("muscle").to_numpy()
176 
177 # Plot the terms of the metabolics model, and compare the metabolics model's
178 # mechanical work rate to the mechanical work rate computed using the
179 # MuscleAnalysis.
180 plot = False
181 # The following environment variable is set during automated testing.
182 if os.getenv('OPENSIM_USE_VISUALIZER') != '0':
183  try:
184  import pylab as pl
185  plot = True
186  except:
187  print('Skipping plotting')
188 
189 if plot:
190  pl.plot(time, force * -velocity,
191  label='active_fiber_force * fiber_velocity', lw=4)
192 
193  table_metabolics = osim.TimeSeriesTable("analyze_ProbeReporter_probes.sto")
194  time = table_metabolics.getIndependentColumn()
195  metabolics_total_rate = table_metabolics.getDependentColumn(
196  "metabolics_TOTAL").to_numpy()
197  pl.plot(time, metabolics_total_rate, label='total metabolic rate')
198 
199  mech_work_rate = table_metabolics.getDependentColumn(
200  "mechanical_work_rate_TOTAL").to_numpy()
201  pl.plot(time, mech_work_rate, label='mechanical work rate')
202 
203  activation_maintenance_rate = table_metabolics.getDependentColumn(
204  "activation_maintenance_rate_TOTAL").to_numpy()
205  pl.plot(time, activation_maintenance_rate,
206  label='activation maintenance rate')
207 
208  shortening_rate = table_metabolics.getDependentColumn(
209  "shortening_rate_TOTAL").to_numpy()
210  pl.plot(time, shortening_rate, label='shortening rate')
211 
212  basal_rate = table_metabolics.getDependentColumn(
213  "basal_rate_TOTAL").to_numpy()
214  pl.plot(time, basal_rate, label='basal rate')
215  pl.legend()
216  pl.show()
217