Moco Inverse Clarification of Results

[Matthias Kindler]

Hello to all,

I already posted this in the OpenSim forum, but I wasn't aware that there also is this forum! Sorry for any potential confusion..

I'm trying to calculate muscle forces for a specific motion with MocoInverse. I based my code on the examples I found in the docs (https://simtk.org/api_docs/opensim/api_docs/).
Unfortunately, when evaluating the results, I noticed that the muscle activation was the same as the force, both in the range 0 to 1. So first, I figured that maybe I need to allow for higher max controls with dgf.setMaxControl(max_value) because I noticed that the default was set to 1, but that set my forces and activations to the range 0 to max_value. I am using DeGrooteFregley muscles.

My questions are..
1) How can the force and the activation be the same? Does this depend on the muscle model (DGF)?
2) Maybe I misunderstood the concept of activation and the generated force... What exactly does setMaxControl do? Is the control the same as the activation and the force as well?

Thank you!

[Nicholas Bianco]

Hi Matthias,

How are you computing muscle force? For typical muscles with optimal forces in the 100s or 1000s of Newtons, you should be producing larger muscle forces than 1.

setMaxControl() will set the bounds on the muscle excitation. Muscle activation dynamics will determine the activation value based on the excitation value. Muscle force will be the optimal force multiplied by the activation.

Best,
Nick

[Matthias Kindler]

Hi Nick,

thanks for your answer.

I built a lumbar spine model with 12 DGF muscle fascicles in this manner: (just as example)

Code: Select all

dgf = osim.DeGrooteFregly2016Muscle()
dgf.set_max_isometric_force(79.58)
dgf.set_optimal_fiber_length(0.1034)
dgf.set_tendon_slack_length(0.0363)
dgf.set_tendon_strain_at_one_norm_force(0.10) 
dgf.set_ignore_activation_dynamics(False) 
dgf.set_ignore_tendon_compliance(True) 
dgf.set_fiber_damping(0.01) # check value in original
dgf.set_tendon_compliance_dynamics_mode("implicit")
dgf.set_max_contraction_velocity(10) 
dgf.set_pennation_angle_at_optimal(0.18675022)
dgf.set_ignore_passive_fiber_force(True)
dgf.set_active_force_width_scale(1.5)
dgf.setMinControl(0.01)
dgf.setMaxControl(30)
dgf.set_default_activation(0.01) 
dgf.addNewPathPoint("PS_L5_TP_r-P1", b_lumbar5, 
                                    osim.Vec3(-0.00562001, 0.02913976, 0.02528807))
dgf.addNewPathPoint("PS_L5_TP_r-P2", b_pelvis, 
                                    osim.Vec3(-0.02380000, -0.0570000, 0.07590000))
                                   

I want to do a MocoInverse study to get the muscle forces for a specific kinematics table. This is roughly how I do it:

Code: Select all

 muscleDrivenModel = helpers.getMuscleDrivenModel()
 inverse = osim.MocoInverse()
 modelProcessor = osim.ModelProcessor(muscleDrivenModel)
 modelProcessor.append(osim.ModOpAddReserves(0.01))
 inverse.setModel(modelProcessor)
 tableProcessor = osim.TableProcessor("kinematics.sto")
 
 inverse.setKinematics(tableProcessor)
 inverse.set_initial_time(0)
 inverse.set_final_time(3.0)
 inverse.set_mesh_interval(0.06)
 inverse.set_convergence_tolerance(1e-4)
 inverse.set_constraint_tolerance(1e-4)
  
 inverse.set_kinematics_allow_extra_columns(True)
 inverse.set_minimize_sum_squared_activations(True)
 inverse.append_output_paths('.*normalized_fiber_length')
 inverse.append_output_paths('.*passive_force_multiplier')
 inverseSolution = inverse.solve()
 

If I understood you correctly, muscle excitation determines muscle activation which in turn multiplied by the optimal force will result in the muscle force...
How do I determine the correct maxControl, then?

Another thing, I'm confused about is the reserve.. I have set it to 0.01 but it always dominates more like in the hundreds (see report

July_11_1145_inverse_moco_muscle_report.pdf

(42.3 KiB) Downloaded 40 times

)... Am I implementing this correctly?

Thanks for your help!

-Matthias

[Nicholas Bianco]

Hi Matthias,

We typically bound muscle controls between [0, 1], so that at max activation, the muscle is producing its maximum isometric force. So I change the following:

Code: Select all

dgf.setMinControl(0.0)
dgf.setMaxControl(1.0)

You can optionally leave the minimum control to 0.01 if solving for muscle tendon dynamics, which has a singularity when muscle activation goes to zero.

The reserves you have in the model have strength 0.01, but they have no bounds. The plots in the report you attached is the value of the reserve actuator control value. The force produced by each reserve is the control value times the strength, so the reserve torques here are still relatively small. If you would like to bound the reserves to [-1, 1], you can do the following:

Code: Select all

 modelProcessor.append(osim.ModOpAddReserves(0.01, 1))

Best,
Nick