This example predicts the trajectory to swing a double pendulum from horizontal to vertical, and then tracks the resulting trajectory.
21 import opensim
as osim
24 This file performs the following problems using a
25 double pendulum model:
27 1. predict an optimal trajectory (and controls),
28 2. track the states from the optimal trajectory, and
29 3. track the marker trajectories from the optimal trajectory.
35 if os.getenv(
'OPENSIM_USE_VISUALIZER') ==
'0':
40 def createDoublePendulumModel():
42 model.setName(
"double_pendulum")
46 b0 = osim.Body(
"b0", 1, osim.Vec3(0), osim.Inertia(1))
48 b1 = osim.Body(
"b1", 1, osim.Vec3(0), osim.Inertia(1))
52 m0 = osim.Marker(
"m0", b0, osim.Vec3(0))
53 m1 = osim.Marker(
"m1", b1, osim.Vec3(0))
58 j0 = osim.PinJoint(
"j0", model.getGround(), osim.Vec3(0), osim.Vec3(0),
59 b0, osim.Vec3(-1, 0, 0), osim.Vec3(0))
60 q0 = j0.updCoordinate()
62 j1 = osim.PinJoint(
"j1",
63 b0, osim.Vec3(0), osim.Vec3(0), b1, osim.Vec3(-1, 0, 0), osim.Vec3(0))
64 q1 = j1.updCoordinate()
69 tau0 = osim.CoordinateActuator()
70 tau0.setCoordinate(j0.updCoordinate())
72 tau0.setOptimalForce(1)
73 model.addComponent(tau0)
75 tau1 = osim.CoordinateActuator()
76 tau1.setCoordinate(j1.updCoordinate())
78 tau1.setOptimalForce(1)
79 model.addComponent(tau1)
82 bodyGeometry = osim.Ellipsoid(0.5, 0.1, 0.1)
83 transform = osim.Transform(osim.Vec3(-0.5, 0, 0))
84 b0Center = osim.PhysicalOffsetFrame(
"b0_center", b0, transform)
85 b0.addComponent(b0Center)
86 b0Center.attachGeometry(bodyGeometry.clone())
87 b1Center = osim.PhysicalOffsetFrame(
"b1_center", b1, transform)
88 b1.addComponent(b1Center)
89 b1Center.attachGeometry(bodyGeometry.clone())
91 model.finalizeConnections()
92 model.printToXML(
"double_pendulum.osim")
96 def solvePrediction():
109 study = osim.MocoStudy()
110 study.setName(
"double_pendulum_predict")
112 problem = study.updProblem()
115 problem.setModel(createDoublePendulumModel())
118 problem.setTimeBounds(0, [0, 5])
122 problem.setStateInfo(
"/jointset/j0/q0/value", [-10, 10], 0)
123 problem.setStateInfo(
"/jointset/j0/q0/speed", [-50, 50], 0, 0)
124 problem.setStateInfo(
"/jointset/j1/q1/value", [-10, 10], 0)
125 problem.setStateInfo(
"/jointset/j1/q1/speed", [-50, 50], 0, 0)
126 problem.setControlInfo(
"/tau0", [-100, 100])
127 problem.setControlInfo(
"/tau1", [-100, 100])
131 ftCost = osim.MocoFinalTimeGoal()
132 ftCost.setWeight(0.001)
133 problem.addGoal(ftCost)
135 finalCost = osim.MocoMarkerFinalGoal()
136 finalCost.setName(
"final")
137 finalCost.setWeight(1000.0)
138 finalCost.setPointName(
"/markerset/m1")
139 finalCost.setReferenceLocation(osim.Vec3(0, 2, 0))
140 problem.addGoal(finalCost)
144 solver = study.initTropterSolver()
145 solver.set_num_mesh_intervals(100)
146 solver.set_verbosity(2)
147 solver.set_optim_solver(
"ipopt")
149 guess = solver.createGuess()
151 guess.setTime([0, 1])
152 guess.setState(
"/jointset/j0/q0/value", [0, -math.pi])
153 guess.setState(
"/jointset/j1/q1/value", [0, 2*math.pi])
154 guess.setState(
"/jointset/j0/q0/speed", [0, 0])
155 guess.setState(
"/jointset/j1/q1/speed", [0, 0])
156 guess.setControl(
"/tau0", [0, 0])
157 guess.setControl(
"/tau1", [0, 0])
158 guess.resampleWithNumTimes(10)
159 solver.setGuess(guess)
162 study.printToXML(
"examplePredictAndTrack_predict.omoco")
165 solution = study.solve()
166 solution.write(
"examplePredictAndTrack_predict_solution.sto")
169 study.visualize(solution)
173 def computeMarkersReference(predictedSolution):
174 model = createDoublePendulumModel()
176 states = predictedSolution.exportToStatesTable()
178 statesTraj = osim.StatesTrajectory.createFromStatesTable(model, states)
180 markerTrajectories = osim.TimeSeriesTableVec3()
181 markerTrajectories.setColumnLabels([
"/markerset/m0",
"/markerset/m1"])
183 for state
in statesTraj:
184 model.realizePosition(state)
185 m0 = model.getComponent(
"markerset/m0")
186 m1 = model.getComponent(
"markerset/m1")
187 markerTrajectories.appendRow(state.getTime(),
188 osim.RowVectorVec3([m0.getLocationInGround(state),
189 m1.getLocationInGround(state)]))
192 markerWeights = osim.SetMarkerWeights()
193 markerWeights.cloneAndAppend(osim.MarkerWeight(
"/markerset/m0", 1))
194 markerWeights.cloneAndAppend(osim.MarkerWeight(
"/markerset/m1", 5))
196 return osim.MarkersReference(markerTrajectories, markerWeights)
199 def solveStateTracking(stateRef):
201 study = osim.MocoStudy()
202 study.setName(
"double_pendulum_track")
204 problem = study.updProblem()
207 problem.setModel(createDoublePendulumModel())
213 finalTime = stateRef.getIndependentColumn()[-1]
214 problem.setTimeBounds(0, finalTime)
215 problem.setStateInfo(
"/jointset/j0/q0/value", [-10, 10], 0)
216 problem.setStateInfo(
"/jointset/j0/q0/speed", [-50, 50], 0)
217 problem.setStateInfo(
"/jointset/j1/q1/value", [-10, 10], 0)
218 problem.setStateInfo(
"/jointset/j1/q1/speed", [-50, 50], 0)
219 problem.setControlInfo(
"/tau0", [-150, 150])
220 problem.setControlInfo(
"/tau1", [-150, 150])
223 stateTracking = osim.MocoStateTrackingGoal()
224 stateTracking.setReference(osim.TableProcessor(stateRef))
225 problem.addGoal(stateTracking)
227 effort = osim.MocoControlGoal()
228 effort.setName(
"effort")
229 effort.setWeight(0.001)
233 solver = study.initTropterSolver()
234 solver.set_num_mesh_intervals(50)
235 solver.set_verbosity(2)
236 solver.set_optim_solver(
"ipopt")
237 solver.set_optim_jacobian_approximation(
"exact")
238 solver.set_optim_hessian_approximation(
"exact")
239 solver.set_exact_hessian_block_sparsity_mode(
"dense")
242 study.printToXML(
"examplePredictAndTrack_track_states.omoco")
245 solution = study.solve()
247 solution.write(
"examplePredictAndTrack_track_states_solution.sto")
250 study.visualize(solution)
254 def solveMarkerTracking(markersRef, guess):
256 study = osim.MocoStudy()
257 study.setName(
"double_pendulum_track")
259 problem = study.updProblem()
262 problem.setModel(createDoublePendulumModel())
268 finalTime = markersRef.getMarkerTable().getIndependentColumn()[-1]
269 problem.setTimeBounds(0, finalTime)
270 problem.setStateInfo(
"/jointset/j0/q0/value", [-10, 10], 0)
271 problem.setStateInfo(
"/jointset/j0/q0/speed", [-50, 50], 0)
272 problem.setStateInfo(
"/jointset/j1/q1/value", [-10, 10], 0)
273 problem.setStateInfo(
"/jointset/j1/q1/speed", [-50, 50], 0)
274 problem.setControlInfo(
"/tau0", [-100, 100])
275 problem.setControlInfo(
"/tau1", [-100, 100])
278 markerTracking = osim.MocoMarkerTrackingGoal()
279 markerTracking.setMarkersReference(markersRef)
280 problem.addGoal(markerTracking)
282 effort = osim.MocoControlGoal()
283 effort.setName(
"effort")
284 effort.setWeight(0.0001)
288 solver = study.initTropterSolver()
289 solver.set_num_mesh_intervals(50)
290 solver.set_verbosity(2)
291 solver.set_optim_solver(
"ipopt")
292 solver.set_optim_jacobian_approximation(
"exact")
293 solver.set_optim_hessian_approximation(
"exact")
294 solver.set_exact_hessian_block_sparsity_mode(
"dense")
296 solver.setGuess(guess)
299 study.printToXML(
"examplePredictAndTrack_track_markers.omoco")
302 solution = study.solve()
304 solution.write(
"examplePredictAndTrack_track_markers_solution.sto")
307 study.visualize(solution)
313 optimalTrajectory = solvePrediction()
315 markersRef = computeMarkersReference(optimalTrajectory)
317 trackedSolution = solveStateTracking(optimalTrajectory.exportToStatesTable())