Issue related to DeGrooteFregly2016Muscle computeEquilibrium

[Mohammadreza Rezaie]

Hi, I'm updating coordinates value and speed as well as muscles activation at each time frame for further analysis related to muscles. Everything works well with Millard2012EquilibriumMuscle. But when I replace the muscles by DeGrooteFregly2016Muscle, it raises the following error at a specific time and a specific muscle:

muscle.computeEquilibrium(state)
File F:\Python311\Lib\site-packages\opensim\simulation.py:31345 in computeEquilibrium
return _simulation.Muscle_computeEquilibrium(self, s)
RuntimeError: std::exception in 'void OpenSim::Muscle::computeEquilibrium(SimTK::State &) const': Function has same sign at bounds of 0.0 and 5.0.
Thrown at CommonUtilities.cpp:164 in solveBisection().

I did a sensitivity analysis and found that in a specific pose, activations less than 0.2 may raise that error.

Any idea why this is happening and how I can fix it? Thank you in advance.

[Nicholas Bianco]

Hi Mohammadreza,

You are setting the activation for all muscles? I know from experience that this muscle is less stable in explicit dynamics mode if you don't have a good guess for the state variables (activation and normalized tendon force). If you have a solve Moco problem in the form of a MocoTrajectory, you could plot tendon force using that. Otherwise, I would recommend the Millard muscle for this application (unless you absolutely need the DGF muscle for some reason).

Best,
Nick

[Mohammadreza Rezaie]

Thanks Nick for your response.

Yes, I'm updating coordinates pose and velocity and muscles activation in a loop (for calculating joint contact force), and I'm not using Moco. I got the activations from Muscle Redundancy Solver (in Python) and since it uses De Groote et al. 2016 muscle mathematical curves, I really need DeGrooteFregly2016Muscle.

One more question (out of the scope of this post): If I use joint reaction tool with muscle force file, does the model of muscles (Millard or DeGroote) matter?

Thanks for your help.

[Nicholas Bianco]

Hi Mohammadreza,

You should be able to extract the tendon forces directly from KULeuven Muscle Redundancy Solver if you're already using that tool. I would recommend that approach so you don't have to worry about converting the outputs from their tool into OpenSim.

-Nick

[Mohammadreza Rezaie]

Hi Nick, thanks for your response. Yes, I have both activation and fiber force (which is in equilibrium with tendon force), and eventually I need to use OpenSim Joint Reaction tool.

Please correct me, in addition to setActivation, I need to update DeGrooteFregly2016Muscle tendon force via setNormalizedTendonForce as well. But I'm getting the same error. Should I change anything else in this muscle?

One more question (out of the scope of this post): If I use joint reaction tool with muscle force file, does the model of muscles (Millard or DeGroote) matter?

I also tested this option in the GUI. It works without error when Solve for equilibrium for actuator states is unchecked. When it is checked, I see the error but it doesn't interrupt the tool and the outputs are identical (with and without Solve for equilibrium ...). Please let me know your thoughts on this as well.

Thank you again for your help.

[Nicholas Bianco]

Hi Mohammadreza,

The OpenSim DeGrooteFregly2016Muscle uses activation and normalized tendon force as a state (not fiber force or fiber length), so you might need to resolve your problem in Muscle Redundancy Solver to get a compatiable states trajectory.

Once you do that, you might consider using the SimulationUtilities::analyze() function to compute the joint reaction forces. It could give you more flexible in creating the states trajectory and computing the contact forces over the Joint Reaction tool (as long as you're certain you've specified the states trajectory correctly).

-Nick

[Mohammadreza Rezaie]

Hi Nick,

I appreciate your help. SimulationUtilities::analyze seems fantastic and I'm going to use it. Thanks for introducing.

I created two TimeSeriesTable(s) as follows:

1. state variables containing coordinates value and speed, and muscles activation and normalized_tendon_force in the same order as in the model.getStateVariableNames()
2. control variables containing muscle activation in the same order as in the createControlNamesFromModel(model)

(I will add actuators as well after I could get it working). This is my analyze snippet:

Code: Select all

reporter = osim.TableReporterSpatialVec()

outputVariables = osim.StdVectorString()
for joint in model.getJointSet():
    outputVariables.append(joint.getAbsolutePathString()+'|reaction_on_child')
    reporter.addToReport(joint.getOutput('reaction_on_child'))

model.addComponent(reporter)
model.finalizeConnections()
model.initSystem()

osim.analyze(model, stateVariables, controlVariables, outputVariables)

Without reporter, it returns an empty table and I get this warning:
Warning: No outputs were connected to 'reporter' of type TableReporter__double_. You can connect outputs by calling addToReport().
Therefore, I added a reporter component to the model. If I'm not mistaken, I think I should use TableReporterSpatialVec() because joint.getOutput('reaction_on_child').getTypeName() return 'SpatialVec'. But the code raises this error:
AttributeError: 'TableReporterSpatialVec' object has no attribute 'addToReport'

Could you please guide me on this as well? Is my API implementation correct?

Thanks again,
-Mohammadreza

[Nicholas Bianco]

Hi Mohammadreza,

In scripting, you can use the function 'analyzeSpatialVec()' to compute SpatialVec quantities from Outputs. It will create a reporter internally return a table of type TimeSeriesTableSpatialVec, no need to attach a reporter yourself:

Code: Select all

outputVariables = osim.StdVectorString()
for joint in model.getJointSet():
    outputVariables.append(joint.getAbsolutePathString()+'|reaction_on_child')
outputTable = osim.analyzeSpatialVec(model, stateVariables, controlVariables, outputVariables)

-Nick

[Mohammadreza Rezaie]

Hi Nick, thanks again for your help.

1. The append method for StdVectorString didn't get it working (the output table had no column). I used push_back('.*reaction_on_child') instead and the code worked without error. But after flattening and writing the table, there are nans for all columns. Any idea why this is happening?

2. I created state table based on model.getStateVariableNames(). Is it correct?

3. I created control table based on osim.createControlNamesFromModel(model) and put muscle activation for each. Is it correct?

4. Doesn't this type of analysis require external load file (ground reaction force)? I added the file as a component to the model. But there are still nans.

Thank you.

[Nicholas Bianco]

Hi Mohammadreza,

1. Try .add() if .append() didn't work. I'm not sure if that's why you're getting NaNs, but hopefullly that fixes things.

2. That should work, just make sure that the names are the full path to each state variable.

3. That should also work, but note that activations are states, not controls. Muscle excitations are controls.

4. The model should exactly match the model you used in your original simulations. The Muscle Redundancy Solver computes joint moments and uses those to formulate the OCP, but for an equivalent in OpenSim/Moco, yes, you will need to apply the external loads to the model. You can use the ModelOperator ModOpAddExternalLoads with a ModelProcessor to do this.

Double check that there are not locked coordinates in your model, that is often the cause of NaNs in OpenSim simulations/analyses.

-Nick