EMG Tracking missing states

[Molly Shepherd]

Hi all,

I have recently been incorporating EMG tracking into my MocoInverse problem. I am having no issues actually getting a tracking solution. However, I have noticed that the EMG tracking solution does not contain any of the jointset states (either value or speed). This is an issue as I am trying to use this EMG tracking solution to run a JR analysis.

The same issue appears to have been mentioned on this thread:
viewtopicPhpbb.php?f=1815&t=14650&p=42757&start=0&view=
however, there did not appear to be a solution.

Any help would be great, thanks!
Molly

Full code for moco inverse and emg tracking:

Code: Select all

import org.opensim.modeling.*;
% Construct the MocoInverse tool.
inverse = MocoInverse();

% Construct a ModelProcessor and set it on the tool. The default
% muscles in the model are replaced with optimization-friendly
% DeGrooteFregly2016Muscles, and adjustments are made to the default muscle
% parameters.
modelProcessor = ModelProcessor('DDH_0Run01_postRRA.osim');
modelProcessor.append(ModOpAddExternalLoads('DDH_0Run01_RRA_GRF.xml'));

modelProcessor.append(ModOpReplaceMusclesWithDeGrooteFregly2016());
% Only valid for DeGrooteFregly2016Muscles.
modelProcessor.append(ModOpIgnorePassiveFiberForcesDGF());
% modelProcessor.append(ModOpIgnoreTendonCompliance());
modelProcessor.append(ModOpUseImplicitTendonComplianceDynamicsDGF());
% modelProcessor.append(ModOpTendonComplianceDynamicsModeDGF('opimplicit'));
% Only valid for DeGrooteFregly2016Muscles.
modelProcessor.append(ModOpScaleActiveFiberForceCurveWidthDGF(1.5));
modelProcessor.append(ModOpAddReserves(10, 100,true));

jointNames = StdVectorString();
jointNames.add('mtp_l');
jointNames.add('mtp_r');

modelProcessor.append(ModOpReplaceJointsWithWelds(jointNames));
inverse.setModel(modelProcessor);



% Construct a TableProcessor of the coordinate data and pass it to the
% inverse tool. TableProcessors can be used in the same way as
% ModelProcessors by appending TableOperators to modify the base table.
% A TableProcessor with no operators, as we have here, simply returns the
% base table.
inverse.setKinematics(TableProcessor('DDH_SS_postRRA_states.sto'));
% inverse.append_output_paths('.*fiber_length');
% Initial time, final time, and mesh interval.
inverse.set_initial_time(1.29);
inverse.set_final_time(2.01);
inverse.set_mesh_interval(0.02);

% By default, Moco gives an error if the kinematics contains extra columns.
% Here, we tell Moco to allow (and ignore) those extra columns.
inverse.set_kinematics_allow_extra_columns(true);
inverse.set_minimize_sum_squared_activations(true);
% constrains muscle residuals/error
inverse.set_constraint_tolerance(1e-05);
% when change from an iteration is smaller than tolerance, iterations terminated
inverse.set_convergence_tolerance(1e-04);

% solution = inverse.solve();
% solution.getMocoSolution().write('MocoInverse_0Run01_solution.sto');

% get MocoStudy
study = inverse.initialize();
problem = study.updProblem();

EMG = importdata('EMG.mat');
emgReference = osimTableFromStruct(EMG.DDH);

% Part 3b: Create a MocoControlTrackingGoal, set its weight, and provide
% the EMG data as the tracking reference. We also need to specify the 
% reference labels for the four muscles whose EMG we will track.
tracking = MocoControlTrackingGoal('emg_tracking');
tracking.setWeight(5);
tracking.setReference(TableProcessor(emgReference));
tracking.setReferenceLabel('/forceset/bflh140_r', 'RBFLH');
tracking.setReferenceLabel('/forceset/ercspn_r', 'RES');
tracking.setReferenceLabel('/forceset/gasmed_r', 'RGASMED');
tracking.setReferenceLabel('/forceset/glmax2_r', 'RGMAX');
tracking.setReferenceLabel('/forceset/glmed2_r', 'RGMED');
tracking.setReferenceLabel('/forceset/recfem_r', 'RRECFEM');
tracking.setReferenceLabel('/forceset/tfl_r', 'RTFL');
tracking.setReferenceLabel('/forceset/vaslat140_r', 'RVASLAT');


% Part 3c: The EMG signals in the tracking are all normalized to have
% a maximum value of 1, but the magnitudes of the excitations from the 
% effort minimization solution suggest that these signals should be
% rescaled. Use addScaleFactor() to add a MocoParameter to the problem that
% will scale the reference data for the muscles in the tracking cost.


% Part 3d: Add the tracking goal to the problem.
problem.addGoal(tracking)

% Part 3e: Update the MocoCasADiSolver with the updated MocoProblem using 
% resetProblem().
solver = MocoCasADiSolver.safeDownCast(study.updSolver());
solver.resetProblem(problem);

% Part 3f: Tell MocoCasADiSolver that the MocoParameters we added to the 
% problem via addScaleFactor() above do not require initSystem() calls on
% the model. This provides a large speed-up.
solver.set_parameters_require_initsystem(false);

% Solve the problem and write the solution to a Storage file.
solution = study.solve();
solution.write('trackingSolution.sto');

[Aaron Fox]

Hi Molly,

We ran into this problem recently as well trying to use a MocoInverse solution for muscle and JRF analyses. Our solution was to generate a states table from the RRA/IK results that we were feeding into the MocoInverse problem, and then there is a function that allows you to insert a states trajectory into your solution

Code: Select all

solution.insertStatesTrajectory(kinematicsTable);

There's a bit of fiddling around to get everything in the correct format, but once you understand the process it can be done in a few lines of code.

Nick or some other Moco users however might be able to give a more efficient response.

Aaron

[Molly Shepherd]

Thanks Aaron,

I will definitely try manipulating the problem this way.

Molly

[Nicholas Bianco]

Hi Molly,

Aaron's approach should work, but another alternative is to obatain the model from the problem configured by MocoInverse. It will include the PrescribedMotion constraints which eliminate the coordinate values and speeds from the State. Therefore, when you use this model for any analyses, it will be consistent with the trajectory from the EMG tracking problem.

Code: Select all

study = inverse.initialize();
problem = study.updProblem();
model = problem.getModel();

Best,
Nick