Coordinate Actuator Outputs

[Francine Leung]

I am writing code based on this section of the example3DWalking file. I want to see the torque values generated by the coordinate actuators that are driving the model (in place of the muscles).

Code: Select all

function torqueDrivenMarkerTracking()

import org.opensim.modeling.*;

% Create and name an instance of the MocoTrack tool.
track = MocoTrack();
track.setName("torque_driven_marker_tracking");

% Construct a ModelProcessor and add it to the tool. ModelProcessors
% accept a base model and allow you to easily modify the model by appending
% ModelOperators. Operations are performed in the order that they are
% appended to the model.
% Create the base Model by passing in the model file.
modelProcessor = ModelProcessor("subject_walk_armless.osim");
% Add ground reaction external loads in lieu of a ground-contact model.
modelProcessor.append(ModOpAddExternalLoads("grf_walk.xml"));
% Remove all the muscles in the model's ForceSet.
modelProcessor.append(ModOpRemoveMuscles());
% Add CoordinateActuators to the model degrees-of-freedom. This
% ignores the pelvis coordinates which already have residual 
% CoordinateActuators.
modelProcessor.append(ModOpAddReserves(250));
track.setModel(modelProcessor);

% Use this convenience function to set the MocoTrack markers reference
% directly from a TRC file. By default, the markers data is filtered at
% 6 Hz and if in millimeters, converted to meters.
track.setMarkersReferenceFromTRC("marker_trajectories.trc");

% There is marker data in the 'marker_trajectories.trc' associated with
% model markers that no longer exists (i.e. markers on the arms). Set this
% flag to avoid an exception from being thrown.
track.set_allow_unused_references(true);

% Increase the global marker tracking weight, which is the weight
% associated with the internal MocoMarkerTrackingGoal term.
track.set_markers_global_tracking_weight(10);

% Increase the tracking weights for individual markers in the data set 
% placed on bony landmarks compared to markers located on soft tissue.
markerWeights = MocoWeightSet();
markerWeights.cloneAndAppend(MocoWeight("R.ASIS", 20));
markerWeights.cloneAndAppend(MocoWeight("L.ASIS", 20));
markerWeights.cloneAndAppend(MocoWeight("R.PSIS", 20));
markerWeights.cloneAndAppend(MocoWeight("L.PSIS", 20));
markerWeights.cloneAndAppend(MocoWeight("R.Knee", 10));
markerWeights.cloneAndAppend(MocoWeight("R.Ankle", 10));
markerWeights.cloneAndAppend(MocoWeight("R.Heel", 10));
markerWeights.cloneAndAppend(MocoWeight("R.MT5", 5));
markerWeights.cloneAndAppend(MocoWeight("R.Toe", 2));
markerWeights.cloneAndAppend(MocoWeight("L.Knee", 10));
markerWeights.cloneAndAppend(MocoWeight("L.Ankle", 10));
markerWeights.cloneAndAppend(MocoWeight("L.Heel", 10));
markerWeights.cloneAndAppend(MocoWeight("L.MT5", 5));
markerWeights.cloneAndAppend(MocoWeight("L.Toe", 2));
track.set_markers_weight_set(markerWeights);

% Initial time, final time, and mesh interval. The number of mesh points
% used to discretize the problem is computed internally using these values.
track.set_initial_time(0.81);
track.set_final_time(1.65);
track.set_mesh_interval(0.05);

I have tried using Force Reporter and creating an output table using analyze(). This only output the forces generated by the reserves in the original model, not the modified one. How can I see the forces and torques generated by the actuators created with ModOpAddReserves()? Thanks.

Code: Select all

% Add ForceReporter to the model
model = modelProcessor.process();
forceReporter = ForceReporter(model);
model.updAnalysisSet().adoptAndAppend(forceReporter);

% Solve! Use study.solve() to skip visualizing.
study = track.initialize(); % turn MocoTrack into a MocoStudy
solution = study.solve();
study.visualize(solution);
solution.write('MocoTrack_markertracking_solution.sto');

% Write the solution 
solution.unseal(); % need this to be able to write the solution - otherwise it says 'This trajectory is sealed, to force you to acknowledge the solver failed; call unseal() to gain access'
outputs = StdVectorString(); 
outputs.add('.*power');
outputs.add('.*speed');
% '.*' outputs everything available
% e.g. '.*length' matches anything that ends with length, to output everything associated with addbrev: outputs.add('/forceset/addbrev_r.*');
% actuation is the force/torque value in N/Nm.
% speed is in rad/s for rotational coordinates

% Fill a table with the selected outputs and create STO file
outputTable = study.analyze(solution, outputs);
STOFileAdapter().write(outputTable,'forces_and_torques.sto'); % write outputs (whatever is stated above e.g. actuation)
controlsTable = solution.exportToControlsTable();
STOFileAdapter().write(controlsTable, 'controls.sto');

% Now run a forward dynamics simulation using RKM with known controls
% Reinitialize the system to account for the added ForceReporter
state = model.initSystem();

% Use a TableProcessor to set the controls during the simulation
controlSet = TableProcessor(controlsTable);

%% Does not work after this line
% Set the controls from the controls table to the model
model.setControls(state, controlSet);

% Create a Manager to run the simulation
manager = Manager(model);
manager.setInitialTime(0.81);
manager.setFinalTime(1.65);

% Integrate the system from initial to final time
manager.integrate(state);

% Get the forces table from the ForceReporter
forcesTable = forceReporter.getForcesTable();

% Write the forces and torques to an STO file
STOFileAdapter.write(forcesTable, 'internal_joint_forces.sto');
disp('Internal joint forces written to internal_joint_forces.sto')

[Nicholas Bianco]

Hi Francine,

The torque value of each CoordinateActuator is simply the control value of the actuator in the solution multiplied by the optimal force of the actuator.

Code: Select all

modelProcessor.append(ModOpAddReserves(250));

This ModelOperator adds CoordinateActuators with optimal force of 250 N-m.

Hope that helps!
-Nick