API  4.3
For MATLAB, Python, Java, and C++ users
exampleSquatToStand_answers.py

This is an example that predicts a squat-to-stand movement and optimizes the stiffness of an assistive passive device.

1 
3 import opensim as osim
4 import exampleSquatToStand_helpers as helpers
5 import mocoPlotTrajectory as plot
6 import os
7 import numpy as np
8 torqueDrivenModel = helpers.getTorqueDrivenModel()
9 muscleDrivenModel = helpers.getMuscleDrivenModel()
10 
11 
13 study = osim.MocoStudy()
14 
15 # Part 1b: Initialize the problem and set the model.
16 problem = study.updProblem()
17 problem.setModel(torqueDrivenModel)
18 
19 # Part 1c: Set bounds on the problem.
20 #
21 # problem.setTimeBounds(initial_bounds, final_bounds)
22 # problem.setStateInfo(path, trajectory_bounds, inital_bounds, final_bounds)
23 #
24 # All *_bounds arguments can be set to a range, [lower upper], or to a
25 # single value (equal lower and upper bounds). Empty brackets, [], indicate
26 # using default bounds (if they exist). You may set multiple state infos at
27 # once using setStateInfoPattern():
28 #
29 # problem.setStateInfoPattern(pattern, trajectory_bounds, inital_bounds, ...
30 # final_bounds)
31 #
32 # This function supports regular expressions in the 'pattern' argument;
33 # use '.*' to match any substring of the state/control path
34 # For example, the following will set all coordinate value state infos:
35 #
36 # problem.setStateInfoPattern('/path/to/states/.*/value', ...)
37 
38 # Time bounds
39 problem.setTimeBounds(0, 1)
40 
41 # Position bounds: the model should start in a squat and finish
42 # standing up.
43 problem.setStateInfo('/jointset/hip_r/hip_flexion_r/value',
44  [-2, 0.5], -2, 0)
45 problem.setStateInfo('/jointset/knee_r/knee_angle_r/value',
46  [-2, 0], -2, 0)
47 problem.setStateInfo('/jointset/ankle_r/ankle_angle_r/value',
48  [-0.5, 0.7], -0.5, 0)
49 
50 # Velocity bounds: all model coordinates should start and end at rest.
51 problem.setStateInfoPattern('/jointset/.*/speed', [], 0, 0)
52 
53 # Part 1d: Add a MocoControlCost to the problem.
54 problem.addGoal(osim.MocoControlGoal('myeffort'))
55 
56 # Part 1e: Configure the solver.
57 solver = study.initCasADiSolver()
58 solver.set_num_mesh_intervals(25)
59 solver.set_optim_convergence_tolerance(1e-4)
60 solver.set_optim_constraint_tolerance(1e-4)
61 
62 if not os.path.isfile('predictSolution.sto'):
63  # Part 1f: Solve! Write the solution to file, and visualize.
64  predictSolution = study.solve()
65  predictSolution.write('predictSolution.sto')
66  study.visualize(predictSolution)
67 
68 
69 
73 tableProcessor = osim.TableProcessor('predictSolution.sto')
74 tableProcessor.append(osim.TabOpLowPassFilter(6))
75 
76 # Part 2b: Add a MocoStateTrackingCost to the problem using the states
77 # from the predictive problem (via the TableProcessor we just created).
78 # Enable the setAllowUnusedReferences() setting to ignore the controls in
79 # the predictive solution.
80 tracking = osim.MocoStateTrackingGoal()
81 tracking.setName('mytracking')
82 tracking.setReference(tableProcessor)
83 tracking.setAllowUnusedReferences(True)
84 problem.addGoal(tracking)
85 
86 # Part 2c: Reduce the control cost weight so it now acts as a regularization
87 # term.
88 problem.updGoal('myeffort').setWeight(0.001)
89 
90 # Part 2d: Set the initial guess using the predictive problem solution.
91 # Tighten convergence tolerance to ensure smooth controls.
92 solver.setGuessFile('predictSolution.sto')
93 solver.set_optim_convergence_tolerance(1e-6)
94 
95 if not os.path.isfile('trackingSolution.sto'):
96  # Part 2e: Solve! Write the solution to file, and visualize.
97  trackingSolution = study.solve()
98  trackingSolution.write('trackingSolution.sto')
99  study.visualize(trackingSolution)
100 
101 
102 
104 plot.mocoPlotTrajectory('predictSolution.sto', 'trackingSolution.sto',
105  'predict', 'track')
106 
107 
109 inverse = osim.MocoInverse()
110 
111 # Part 4a: Provide the model via a ModelProcessor. Similar to the TableProcessor,
112 # you can add operators to modify the base model.
113 modelProcessor = osim.ModelProcessor(muscleDrivenModel)
114 modelProcessor.append(osim.ModOpAddReserves(2))
115 inverse.setModel(modelProcessor)
116 
117 # Part 4b: Set the reference kinematics using the same TableProcessor we used
118 # in the tracking problem.
119 inverse.setKinematics(tableProcessor)
120 
121 # Part 4c: Set the time range, mesh interval, and convergence tolerance.
122 inverse.set_initial_time(0)
123 inverse.set_final_time(1)
124 inverse.set_mesh_interval(0.05)
125 inverse.set_convergence_tolerance(1e-4)
126 inverse.set_constraint_tolerance(1e-4)
127 
128 # Allow extra (unused) columns in the kinematics and minimize activations.
129 inverse.set_kinematics_allow_extra_columns(True)
130 inverse.set_minimize_sum_squared_activations(True)
131 
132 # Append additional outputs path for quantities that are calculated
133 # post-hoc using the inverse problem solution.
134 inverse.append_output_paths('.*normalized_fiber_length')
135 inverse.append_output_paths('.*passive_force_multiplier')
136 
137 # Part 4d: Solve! Write the MocoSolution to file.
138 inverseSolution = inverse.solve()
139 solution = inverseSolution.getMocoSolution()
140 solution.write('inverseSolution.sto')
141 
142 # Part 4e: Get the outputs we calculated from the inverse solution.
143 inverseOutputs = inverseSolution.getOutputs()
144 sto = osim.STOFileAdapter()
145 sto.write(inverseOutputs, 'muscleOutputs.sto')
146 
147 
150 device = osim.SpringGeneralizedForce('knee_angle_r')
151 device.setStiffness(50)
152 device.setRestLength(0)
153 device.setViscosity(0)
154 muscleDrivenModel.addForce(device)
155 
156 # Part 5b: Create a ModelProcessor similar to the previous one, using the same
157 # reserve actuator strength so we can compare muscle activity accurately.
158 modelProcessor = osim.ModelProcessor(muscleDrivenModel)
159 modelProcessor.append(osim.ModOpAddReserves(2))
160 inverse.setModel(modelProcessor)
161 
162 # Part 5c: Solve! Write solution.
163 inverseDeviceSolution = inverse.solve()
164 deviceSolution = inverseDeviceSolution.getMocoSolution()
165 deviceSolution.write('inverseDeviceSolution.sto')
166 
167 
168 print('Cost without device: ', solution.getObjective())
169 print('Cost with device: ', deviceSolution.getObjective())
170 
171 # This is a convenience function provided for you. See below for the
172 # implementation.
173 helpers.compareInverseSolutions(inverseSolution, inverseDeviceSolution)