Page 1 of 3
Joint Reaction Analysis
Posted: Thu May 11, 2023 9:52 pm
by nvanden
Hi all,
I've been using Moco, working with simulating rehab in TFA models during stance phases. Previously I'd been looking into how different interventions affected changes in movement patterns, I've recently started looking at how the JRFs change (using code shared by Pasha). I'm wondering if anyone can explain what's goin on "under-the-hood" with this though, particularly 'analyzeSpatialVec'? I see it using the states and controls of the Moco trajectory to find the JRFs, but I'm curious how it does that. In previous OpenSim work I've used the analyze tool to find JRFs with results from a muscle optimization, some kinematics, and an external loads file. I wouldn't expect a 1:1 outcome between Moco and the other workflow but the differences over the same analysis timeframe are a bit alarming.
- HJRF.jpg (30.3 KiB) Viewed 2569 times
The JRF output from Moco seems to missing the loading response from early stance, even though the kinematic solution is pretty close to the experimental data that I use as the input kinematics,
- Kin.jpg (119.2 KiB) Viewed 2569 times
and the predicted GRFs from the contact spheres also look pretty good.
- GRF.jpg (40.21 KiB) Viewed 2569 times
I'm wondering if analyzeSpatialVec is doing some different type of analysis compared with the standard analyze tool and what could be going on to miss the loading during that phase of gait?
Re: Joint Reaction Analysis
Posted: Fri May 12, 2023 5:23 am
by rosshm
Hi Nicholas,
What do the hip muscle forces in early stance for the Moco solution look like?
I have not done hip but I get reasonable-looking knee contact forces in Moco solutions with the method here. Some small changes in syntax may be necessary, this code is from an early version of Moco:
pvb wrote: ↑Sat Mar 20, 2021 2:47 pm
Hi Nick,
Thank you so much, this works! 'flatten()' provides you with a TimeSeriesTable with 6 + 1 columns for each column of TimeSeriesTableSpatialVec. Essentially, this means that each element of a SpatialVec (3 Torque and 3 Force elements), gets unpacked to its "own column', and the extra column is the time column.
As a final question: the default behaviour of "analyzeSpatialVec" is to express the SpatialVecs in the "Ground" (i.e. an inertial) frame of reference, is this correct?
If so, to transform each spatialVec to the desired frame, I must implement something like the following (pseudo) code:
Code: Select all
state = model.realizeDynamics(statesTraj.get(0)) ; %statesTraj is the StatesTrajectory from a prediction solution
Torque_ground = Vec3(JR(1,1), JR(1,2), JR(1,3));
Force_ground = Vec3(JR(1,4), JR(1,5), JR(1,6));
%JR is an array obtained after flattening a TimeSeriesTableSpatialVec and converting
% it to a struct using osimTableToStruct
ground = model.getGround();
child_frame = model.getJointSet().get(0).getChildFrame; %change this to the desired frame for each joint
Torque_local = ground.expressVectorInAnotherFrame(state, Torque_ground, child_frame);
Force_local = ground.expressVectorInAnotherFrame(state, Force_ground, child_frame);
joint_reactions_local(1,1) = Torque_local.get(0);
joint_reactions_local(1,2) = Torque_local.get(1);
joint_reactions_local(1,3) = Torque_local.get(2);
joint_reactions_local(1,4) = Force_local.get(0);
joint_reactions_local(1,5) = Force_local.get(1);
joint_reactions_local(1,6) = Force_local.get(2);
Thank you so much, I think I'm all set to implement my project now!
Pasha
EDIT: different implementation to transform the vectors
Ross
Re: Joint Reaction Analysis
Posted: Fri May 12, 2023 12:28 pm
by nvanden
Hi Ross,
Some of the stance limb hip abductors and flexors are spiking around heel-strike, a couple are starting at a pretty high activation (~0.8) during the first frame that drops off right after.
The code you shared is what I've been using (or adapted from that same discussion chain). My question is more related to how that is doing the JRF analysis, are the GRF produced by the contact spheres affecting that at all? Or is the "weirdness" in the output JRF stemming from a lack of some initial activation bound in the muscles?
Thanks,
Nick
Re: Joint Reaction Analysis
Posted: Fri May 12, 2023 3:08 pm
by rosshm
There is not a lot of documentation on the technical side of JointReaction. You would probably have to find the CPP code for it to really figure out what it is doing. It reads and looks to me like it sums the net joint force with the actuator forces at the joint e.g. muscle forces.
Since your GRF look pretty good I think they or the contact spheres are probably not the source.
Does your simulation have a periodicity constraint? Without that, with at least some of the collocation schemes, the first timestep is unconstrained and you can get strange results there e.g. very-high activations. Can probably also use the various initial-state MocoGoals if periodicity is undesirable.
Ross
Re: Joint Reaction Analysis
Posted: Mon May 15, 2023 11:42 am
by nvanden
Thanks, I haven't been able to find much yet, but I'll keep looking into that.
I know I'd had some weirdness with my periodicity constraint previously, because I'm only running sims over the amputated limb stance phase. I'd tried to adapt the periodicity code from the 'example2Dwalker' since that was also over just a stance but had some issues with adapting to 3D as far as state names. I've been hard-coding my periodicity constraint (manually adding each goal pair) but could very well have done it wrong. I haven't really done anything with initial-state MocoGoals yet, could that be a good way to ensure the initial activations in addition to periodicity?
Nick
Re: Joint Reaction Analysis
Posted: Mon May 15, 2023 3:32 pm
by rosshm
If you are simulating the stance only (and not a full stride) then the periodicity goal will not be appropriate. I would look into the various goals with Initial in the name.
For limb loss since it is presumably asymmetric, the periodicity over even a full step is not applicable, would need to be a full stride.
Nick B can probably give a more specific recommendation. I have only ever done strides!
Ross
Re: Joint Reaction Analysis
Posted: Mon May 15, 2023 4:00 pm
by nbianco
Hi Nick,
If the initial activations are responsible for the large joint loads, you could try adding a MocoInitialActivationGoal to the problem, which constrains the initial activations to be equal to the initial excitations. This way, the initial activations will be appropriately penalized since they are effectively tied to the MocoControlGoal. (MocoInverse problems use this goal by default).
-Nick
Re: Joint Reaction Analysis
Posted: Mon May 15, 2023 9:10 pm
by nvanden
That's good to know, thank you both! I appreciate the help.
Nick
Re: Joint Reaction Analysis
Posted: Wed May 17, 2023 8:58 am
by rosshm
The "initial states" problem here can be a big problem. The way that the system dynamics are enforced as constraints, at least in most collocation schemes that I know of, are with finite difference calculations that involve the current state and the previous state, or something similar. But for the first timestep there is no previous timestep to use, so the initial state is "unconstrained" in this regard. You can specify a specific value or bounds on that value for those initial states, but even then it gets a little thorny. For example if doing a "walk with minimum effort" problem, the optimizer will generally pick the highest allowed value for the initial horizontal velocity, this is free movement for zero effort.
For motions that are, on average, periodic, like gait, this is where the periodicity constraint is helpful. There is still no previous timestep to constrain the initial states, but they are constrained to be equal to the final states, which do have a previous timestep to constrain their dynamics. Basically it is using the second-to-last timestep to inform the dynamics constraint of the initial timestep, in a roundabout way.
Nick, not sure if you have data (tracking targets) on a full stride, but if you do, you might get better results and simpler problems by simulating a full stride, and can avoid making questionable assumptions about the initial states e.g. is assuming the initial activations are equal to the initial excitations a good assumption? Not necessarily.
Ross
Re: Joint Reaction Analysis
Posted: Thu Jun 08, 2023 11:28 am
by nvanden
Hi Ross,
Sorry for the delayed response, went down a bit of a rabbit hole. I do have data from a full IK trial that I select windows of for the tracking problems, where I've been isolating to a stance phase previously. I'm wondering if I'm setting up my periodicity wrong or if I've got the wrong idea with expanding my time range? I've been hard-coding the periodicity goal (syntax for the lumbopelvic DOFs seems weird but I tried to follow this:
https://opensim-org.github.io/opensim-m ... ml#details). If I'm understanding the goal correctly this should be the way to define periodicity?
Code: Select all
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_tilt/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_tilt/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_list/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_list/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_rotation/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/groundPelvis/pelvis_rotation/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_flexion_l/value','/jointset/hip_l/hip_flexion_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_flexion_l/speed','/jointset/hip_l/hip_flexion_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_adduction_l/value','/jointset/hip_l/hip_adduction_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_adduction_l/speed','/jointset/hip_l/hip_adduction_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_rotation_l/value','/jointset/hip_l/hip_rotation_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_l/hip_rotation_l/speed','/jointset/hip_l/hip_rotation_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/knee_l/knee_angle_l/value','/jointset/knee_l/knee_angle_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/knee_l/knee_angle_l/speed','/jointset/knee_l/knee_angle_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/ankle_l/ankle_angle_l/value','/jointset/ankle_l/ankle_angle_l/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/ankle_l/ankle_angle_l/speed','/jointset/ankle_l/ankle_angle_l/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_flexion_r/value','/jointset/hip_r/hip_flexion_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_flexion_r/speed','/jointset/hip_r/hip_flexion_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_adduction_r/value','/jointset/hip_r/hip_adduction_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_adduction_r/speed','/jointset/hip_r/hip_adduction_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_rotation_r/value','/jointset/hip_r/hip_rotation_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/hip_r/hip_rotation_r/speed','/jointset/hip_r/hip_rotation_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/knee_r/knee_angle_r/value','/jointset/knee_r/knee_angle_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/knee_r/knee_angle_r/speed','/jointset/knee_r/knee_angle_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/ankle_r/ankle_angle_r/value','/jointset/ankle_r/ankle_angle_r/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/ankle_r/ankle_angle_r/speed','/jointset/ankle_r/ankle_angle_r/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_extension/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_extension/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_bending/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_bending/speed'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_rotation/value'));
periodicityGoal.addStatePair(MocoPeriodicityGoalPair('/jointset/lumbar/lumbar_rotation/speed'));
However, my output model is jumping backwards. When I comment out the lumbopelvic DOFs, I can move 'normally' but then am I going to have similar issues with not being constrained?
I've also changed the time window to look at a full gait cycle (HS --> HS) but I'm realizing that may not be analogous to a stride (in my head, HS --> contralateral TO)?
Thanks,
Nick