Low Muscle Activation

[Molly Shepherd]

Hi All!

I have been working on improving my Moco Inverse solution for walking/running trials, however I continue to get very low muscle activations (primarily < 0.2). I am using implicit tendon dynamics, and have lowered my reserves as much as possible while still allowing the model to run (100). In addition I have set constraint and convergence tolerances to 1e-03 and 1e-02 respectively based on information I have found on the forum. My residuals are all small, on the order or 10^-5. I did a comparison of the running activations to experimental EMG data collected in lab. The comparison normalized both the MocoInverse activations and the experimental EMG activations peak values to 1 which caused them to align pretty well. However, the walking activations do not match EMG data as well and are quite a bit smaller. I have attached the MocoInverse Walking and Running solutions.

Code is listed below:

Code: Select all

% function solveMocoInverse()
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('DDH09_SS_0Wlk01_postRRA.osim');
modelProcessor.append(ModOpAddExternalLoads('DDH09_0Wlk01_RRA_GRF.xml'));

modelProcessor.append(ModOpReplaceMusclesWithDeGrooteFregly2016());
% Only valid for DeGrooteFregly2016Muscles.
modelProcessor.append(ModOpIgnorePassiveFiberForcesDGF());
% modelProcessor.append(ModOpIgnoreTendonCompliance());
modelProcessor.append(ModOpUseImplicitTendonComplianceDynamicsDGF());
% Only valid for DeGrooteFregly2016Muscles.
modelProcessor.append(ModOpScaleActiveFiberForceCurveWidthDGF(1.5));
modelProcessor.append(ModOpAddReserves(100, 1));
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('DDH09_SS_0Wlk01_states.sto'));
% inverse.append_output_paths('.*fiber_length');
% Initial time, final time, and mesh interval.
inverse.set_initial_time(3.38);
inverse.set_final_time(4.349);
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);
inverse.set_constraint_tolerance(1e-03);
inverse.set_convergence_tolerance(1e-02);

% Solve the problem and write the solution to a Storage file.
solution = inverse.solve();
solution.getMocoSolution().write('MocoInverse_0Wlk01_solution.sto');

% Generate a report with plots for the solution trajectory.
model = modelProcessor.process();
report = osimMocoTrajectoryReport(model, ...
        'MocoInverse_0Wlk01_solution.sto', 'bilateral', true);
% The report is saved to the working directory.
reportFilepath = report.generate();
open(reportFilepath);

Thanks!
Molly

[Aaron Fox]

Hi Molly,

It looks like you're setting the optimal force of the reserves to be 100, correct? If that's the case then I would actually view this as being pretty high - so it's likely that the model is preferentially using these reserves over the muscles (i.e. the high optimal force is at such a low cost because they don't need to be activated much to use).

If I were setting reserve torques in a walking problem like this, I would probably set the optimal forces in reserves around 1 or 2 (and probably no higher than 5). With such a low optimal force, setting the activation bounds to 1 might mean they max out - so you might consider increasing the max control signal they can reach. A balance I often use is to have a low optimal force (i.e. 1 or 2), but allow an infinite control signal for these actuators - this means they can produce as much torque as necessary, yet it comes at a high cost in the optimisation due to the low optimal force.

Hope that makes sense.

Aaron

[Molly Shepherd]

Hi Aaron,

Thanks for the reply. That makes sense as to why the activations are low, unfortunately if I try to lower the reserves at all (even to 90) the problem either takes many hours to run or never converges at all. I'm not sure if there's a way to fix that issue? If it is possible, I may try lowering reserves on a joint by joint basis to see if that helps.

I will also try lowering the optimal force and increasing the activation bounds.

Thanks!

Molly

[Edward Syrett]

Hi Molly,
I had a similar issue, but I found the problem was with the application of my GRF's. I initially had the GRF's positioned at the force plate origin, but I think that was causing large torques on the model, which needed the reserve actuators to counter. I switched my data to where the GRF's are applied at the COP (if you are using the c3dexport Matlab function, that's option 1 rather than 0), and my reserve actuators decreased significantly. I don't know if that helps at all!

Dan

[Aaron Fox]

Hi Molly,

If you change the reserve torques control limits to infinity - then theoretically it should be able to achieve a solution, just with a higher cost function. This might cause issues with your problem converging though due to the difficulty in finding that optimal solution that satisfies your convergence tolerance.

As another suggested on the forum here there may be some things to fix with your experimental data. Have you ran it through anything like RRA to check what the residuals and required torques are like?

Aaron

[Molly Shepherd]

Hi Aaron and Dan,

Yes I have run my data through RRA and my kinematics and residuals all look good so I am not concerned with my GRFs or experimental data.

I can get the problem to converge with the reserve torques set to 100, just not with anything lower. I am trying to output the reserves utilized by each joint to see if I can systematically lower them on a joint by joint basis.

Thanks!
Molly