Predicted GRF and Inverse Dynamics

[Aaron Fox]

Hi Moco Team and Users,

I have a query that those who have used contact spheres to predict GRFs may be able to answer. The predicted GRFs that stem from contact spheres in, for example, tracking or predictive sims in Moco only estimate force values and not position (or CoP) values. To perform something like inverse dynamics or similar the force values would need to appropriately be applied to a body. Is there a specific way that this could be done either within or outside of Moco/OpenSim? Something that's coming to mind at the moment would be to 'centre' it within the positions of the contact spheres on the foot (i.e. average medial/lateral and anterior/posterior position of the spheres), with this averaging weighted by the amount of force being produced by the different spheres (i.e. a sphere producing no force contributes no weight to the positional average).

Has anyone else developed similar or different methods?

Aaron

[Nicholas Bianco]

Hi Aaron,

When reporting the reaction forces and moments from SmoothSphereHalfSpaceForce ground contact components via getRecordValues(), the reaction forces and moments are expressed in ground about the origin of the ground frame. You should be able to compute the center of pressure like you would from force plate forces and moments(http://www.kwon3d.com/theory/grf/cop.html), where the "true" origin is just (0, 0, 0).

-Nick

[Aaron Fox]

Hi Nick,

Thanks for the heads up - looks relatively easy enough to calculate the CoP from the data output.

Aaron

[Ross Miller]

I think you could also use study.analyze() or study.analyzeSpatialVec() to get the resultant joint moments of a solution, although I'm speculating since I haven't tried this myself, and I'm not sure if it would give 3-D moments or just the moments for the model's DoF.

Ross

[Russell Johnson]

Hi All,

I am wondering if anyone had been able to get the study.analyze() method to work to calculate joint moments, as Ross suggested was maybe possible in the above post. I have been trying to calculate the joint moments with Moco solutions. I have a pipeline to calculate moments after calculating centers of pressure, but my ID results don't quite align with what I would expect. I would like to try this other method to check how well it does or does not match with my COP-based calculations. I have been able to calculate muscle-based outputs using study.analyze() (e.g., fiber length or velocity) but haven't quite figured out the syntax to dig into the joint moments.

Thanks,
Russell

[Ross Miller]

Hi Russell,

You can use analyzeSpatialVec to compute what Moco calls the "joint reactions" (forces and moments). There is some sample Matlab code in this thread:

viewtopicPhpbb.php?f=1815&t=13087&p=0&s ... 9795935455

I get odd results using it for moments. Some of them e.g. knee adduction moment look right, but others have strange magnitudes and shapes. I think I may be misunderstanding what a "reaction moment" is, or expressing them in the wrong reference frame, or [third mistake here].

When I compute joint moments by summing the moments from muscle forces and other actuators, I get normal-looking joint moments (example attached for hip/knee/ankle moments from a tracking simulation of walking with a 3-D model). The code I use for this is:

Code: Select all

muscles = model.getMuscles;
coordSet = model.getCoordinateSet();
statesTraj= gaitTrackingSolution.exportToStatesTrajectory(model); % Get state trajectory to loop through it


for i = 1:statesTraj.getSize
    state = statesTraj.get(i-1);
    model.realizeDynamics(state);
    
    % Muscle moment arms
    for j = 1:muscles.getSize
        mus = muscles.get(j-1);
        coord = coordSet.get('ankle_angle_r');
        MusMomentArms_ankle_angle_r(i,j) = mus.computeMomentArm(state, coord);
        coord = coordSet.get('knee_angle_r');
        MusMomentArms_knee_angle_r(i,j) = mus.computeMomentArm(state, coord);
        coord = coordSet.get('hip_flexion_r');
        MusMomentArms_hip_flexion_r(i,j) = mus.computeMomentArm(state, coord);
    end
    
    % Passive joint moments
    force = forceSet.get('PassiveAnkleDamping_r');
    EBCforce = ExpressionBasedCoordinateForce.safeDownCast(force);
    JMpas_ankle_angle_r(i) = EBCforce.calcExpressionForce(state);
    force = forceSet.get('PassiveKneeDamping_r');
    EBCforce = ExpressionBasedCoordinateForce.safeDownCast(force);
    JMpas_knee_angle_r(i) = EBCforce.calcExpressionForce(state);
    force = forceSet.get('PassiveHipDampingX_r');
    EBCforce = ExpressionBasedCoordinateForce.safeDownCast(force);
    JMpas_hip_flexion_r(i) = EBCforce.calcExpressionForce(state);
end

% Joint moments: sum of moments from muscle forces and passive joint moments
JM_ankle_angle_r = sum(MusMomentArms_ankle_angle_r.*MusTendonForces,2) + JMpas_ankle_angle_r';
JM_knee_angle_r = sum(MusMomentArms_knee_angle_r.*MusTendonForces,2) + JMpas_knee_angle_r';
JM_hip_flexion_r = sum(MusMomentArms_hip_flexion_r.*MusTendonForces,2) + JMpas_hip_flexion_r';

Before this you'd have to extract the tendon forces (I call them "MusTendonForces" here) using study.analyze.

Let me know if any of this doesn't work.

Ross

JointMoments.png (64.97 KiB) Viewed 1709 times

[Russell Johnson]

Hi Ross,

Thanks for linking to that thread for the reaction moments - I was able to calculate those as described, but agree that the joint moments there are a bit unexpected and I'm not sure why - wondering if anyone might have some insight as to which part I am missing.

I also used your suggestion to calculate joint moments based on muscle forces x moment arms plus the passive joint moments - these worked perfectly. And these seem to be a lot more realistic for expected joint moments than what I was calculating with my COP-based code.

For others that might reproduce this, make sure to also include a call to get the forceset:

Code: Select all

muscles = model.getMuscles;
coordSet = model.getCoordinateSet();
forceSet = model.getForceSet();
statesTraj= gaitTrackingSolution.exportToStatesTrajectory(model); % Get state trajectory to loop through it
time = gaitTrackingSolution.getTimeMat();

Thanks again - this was incredibly helpful!
Russell

[Pasha van Bijlert]

Hi all,

It is my understanding that the reaction moments and forces are the moments and forces that the joints produce to stop the segments from flying apart (i.e. constraint forces & moments, if you were to set up the EoM manually). These then do not have a direct correspondence with joint moments (i.e. net moment caused by summing inertial forces and rxF for each muscle).

So for a pin joint that can rotate in the Z plane, I would expect no joint reaction Moment in that plane. Muscle forces can affect the joint reaction forces, but are not the only factors that influence them, so I don't think there is a 1 to 1 correspondence. Does this interpretation explain the unexpected results? I suppose these are more general OpenSim musings than Moco-specific.

Best,
Pasha

[Carlos Gonçalves]

Sorry for jumping late in the discussion, but something that I used a lot to compare RRA results to CMC is use the Muscle Analysis inside the "analysis" task. This post helped me a lot to make it work inside Python https://github.com/opensim-org/opensim-core/issues/2076

I always set the tool to calculate the moments (setComputeMoments(true)) and sum the muscle moments for the desired joint.

It seems dump, but since the moment arm is zero for muscles that don't act on the desired joint. Even summing them all, you get the correct value.

Hope it helps.

Best regards.