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.
26 import opensim
as osim
28 def createHangingMuscleModel(ignore_activation_dynamics,
29 ignore_tendon_compliance):
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)
37 joint = osim.SliderJoint(
"joint", model.getGround(), body)
38 coord = joint.updCoordinate()
39 coord.setName(
"height")
40 model.addComponent(joint)
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)
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))
66 probe = osim.Umberger2010MuscleMetabolicsProbe()
67 probe.setName(
"metabolics")
68 probe.addMuscle(
"muscle", 0.5)
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)
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);
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)
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)
107 body.attachGeometry(osim.Sphere(0.05))
109 model.finalizeConnections()
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")
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])
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")
137 initial_equilibrium.setMode(
"cost")
138 initial_equilibrium.setWeight(0.001)
140 problem.addGoal(osim.MocoFinalTimeGoal())
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)
148 solution = study.solve()
149 osim.STOFileAdapter.write(solution.exportToStatesTable(),
150 "exampleHangingMuscle_states.sto")
151 osim.STOFileAdapter.write(solution.exportToControlsTable(),
152 "exampleHangingMuscle_controls.sto")
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")
166 analyze = osim.AnalyzeTool(
"exampleHangingMuscle_AnalyzeTool_setup.xml")
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()
182 if os.getenv(
'OPENSIM_USE_VISUALIZER') !=
'0':
187 print(
'Skipping plotting')
190 pl.plot(time, force * -velocity,
191 label=
'active_fiber_force * fiber_velocity', lw=4)
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')
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')
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')
208 shortening_rate = table_metabolics.getDependentColumn(
209 "shortening_rate_TOTAL").to_numpy()
210 pl.plot(time, shortening_rate, label=
'shortening rate')
212 basal_rate = table_metabolics.getDependentColumn(
213 "basal_rate_TOTAL").to_numpy()
214 pl.plot(time, basal_rate, label=
'basal rate')