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!
Moco Inverse Clarification of Results
- Matthias Kindler
- Posts: 9
- Joined: Tue Dec 06, 2022 1:02 am
- Nicholas Bianco
- Posts: 1005
- Joined: Thu Oct 04, 2012 8:09 pm
Re: Moco Inverse Clarification of Results
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
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
- Posts: 9
- Joined: Tue Dec 06, 2022 1:02 am
Re: Moco Inverse Clarification of Results
Hi Nick,
thanks for your answer.
I built a lumbar spine model with 12 DGF muscle fascicles in this manner: (just as example)
I want to do a MocoInverse study to get the muscle forces for a specific kinematics table. This is roughly how I do it:
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 )... Am I implementing this correctly?
Thanks for your help!
-Matthias
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))
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()
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 )... Am I implementing this correctly?
Thanks for your help!
-Matthias
- Nicholas Bianco
- Posts: 1005
- Joined: Thu Oct 04, 2012 8:09 pm
Re: Moco Inverse Clarification of Results
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:
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:
Best,
Nick
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)
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))
Nick