Fiber Velocity Outputs

[Russell Johnson]

I've been running Moco optimizations, and would like to output fiber velocity from each of the solutions. I've been following the suggested procedure of using study.analyze() to output all the fiber lengths and velocities, as suggested in a previous forum post. I've been finding some weird fiber velocity trajectories when I get them directly from the fiber_velocity command, but when I take the fiber length output and differentiate them, I get much cleaner and more sensible results.

Code: Select all

outputPaths100 = StdVectorString();
outputPaths100.add('.*fiber_velocity')
FiberVelocity100 = study.analyze(prevSolution_100p, outputPaths100);

The velocity trajectories match during certain parts of the gait cycle, but the dramatic peaks are inconsistent. There does seem to be a pattern between limbs, where corresponding muscles between limbs give the similar irregularities (offset by 50%). I'm wondering if I need to do some realizeVelocity or realizeAcceleration command, but I haven't been able to get the correct for that since analyze seems to need a MocoTrajectory but realize needs a StatesTrajectory.

Thanks,
Russell

Muscle Fiber Velocity Fig_v2.jpg (155.94 KiB) Viewed 943 times

[Ross Miller]

Hi Russell,

Interesting problem. I checked for this same thing on a recent simulation I ran, this was walking by a 3-D model with 84 muscles, constraint tolerance = 1e-04 and dt = 0.01 s. There are some small differences between the CC velocities calculated by study.analyze() vs. calculated by finite differencing the CC lengths, biggest is about 0.1 m/s:

Screen Shot 2020-11-12 at 10.38.49 AM.png (747.61 KiB) Viewed 921 times

This was a tracking simulation (guessing the "p" in your prevSolution_100p stands for "predictive"?), not sure if that's a major factor here. The code I used to get the CC lengths and velocities (101x84 matrices) is:

Code: Select all

outputPaths = StdVectorString();
outputPaths.add('.*fiber_velocity');
outputPaths.add('.*fiber_length');
outputTable = study.analyze(gaitTrackingSolution, outputPaths);
opensimMoco.writeTableToFile(outputTable,'MusOut.sto')
temp5 = importdata('MusOut.sto');
MusOut = temp5.data;
MusOutHeaders = temp5.colheaders;

j = 1;
for i = 1:length(MusOutHeaders)
    currentVarName = string(MusOutHeaders(i));
    if endsWith(currentVarName,'|fiber_velocity')
        ivec(j) = i;
        j = j + 1;
    end
end
MusFiberVelocities = MusOut(:,ivec);

j = 1;
for i = 1:length(MusOutHeaders)
    currentVarName = string(MusOutHeaders(i));
    if endsWith(currentVarName,'|fiber_length')
        ivec(j) = i;
        j = j + 1;
    end
end
MusFiberLengths = MusOut(:,ivec);

CC length isn't a state in Moco, so there's no constraint forcing Vcc - d(Lcc)/dt < ctol, but I'm still surprised you get values that differ by ~5 m/s let alone ~5000 m/s.

Based on the timing of the spikes, I would guess this might have something to do with how <DeGrooteFregly2016Muscle> calculates Vcc when activation is ~zero or when Lsec = ~Lslack, something like that? Rather than d(Lcc)/dt, you could try calculating Vcc analytically from the output states and the muscle model and see if that matches the study.analyze() result better (I assume that's what study.analyze() is doing). DeGroote et al. (2016) has the full algorithm for the muscle model although I think some small changes were made in the Moco implementation, they might be detailed in the Moco documentation.

Ross

[Nicholas Bianco]

Hi Russell,

Thanks for pointing this out, it would be good to get to the bottom of this issue. Are you generating these results with tendon compliance enabled? In Moco, normalized tendon force is the state variable when tendon compliance is enabled, and both fiber length and fiber velocity are computed based on the normalized tendon force state. (This is done using the inverse of the tendon force-length and fiber force-velocity curves).

I don't think it's a realization issue, analyze() realizes to the Velocity stage. As Ross suggests, it could be helpful to plot activation and normalized tendon length alongside fiber length and fiber velocity to see if there's a correlation with times of low force production. You can also plot tendon force, either via analyze(), or the state variable directly (if tendon compliance is enabled).

Best,
-Nick

[Russell Johnson]

Hi Nick and Ross,

Thanks for the feedback. I had also noticed that, even for the times where the fiber velocity output seems reasonable, there are some differences between that and the derivative of fiber length (+/- 0.1 m/s), especially noticeable around the peaks in velocity. Obviously for now, I have been more focused on the big peaks to begin with. As Ross guessed, these results are from a "purely predictive" optimization, walking at 1.00 m/s, minimizing controls cubed.

I've plotted the muscle lengths below for the same result as my first post in this thread (this time only the right side), with the shaded gray regions the areas where I got the big spikes in velocity. One of my thoughts was for fiber lengths that were very small or large, there could be difficulty calculating velocity, but there doesn't seem to be a consistent pattern with the length data. Except for one instance in the Vastus Int, each muscle is lengthening at the time of the spikes, but maybe that is coincidental.

Muscle Fiber Length.jpg (107.45 KiB) Viewed 874 times

The plot for normalized tendon force and muscle activation is here, with the same shaded regions. The activation is very low in most of the regions where the fiber velocity spikes, except for the TA and EHL. It overlaps a bit with the ramping up of activation. Let me know if you'd like to see anything else visualized.

Muscle Activation and Tendon Force.jpg (140.01 KiB) Viewed 874 times

I assume the fiber damping is one potential culprit, in my model the fiber damping is set to default of 0.01. I have not tried different values here yet.

One follow up question: would the discrepancy in fiber velocity output here have a negative effect on a metabolic cost calculation, since it needs fiber velocity as an input to calculate fiber work? I have been calculating this using the states and controls file, and my values are actually pretty reasonable. I assume that given these large peaks in velocity, if this error was a factor then it would be pretty obvious in the metabolic cost output. But perhaps it is because I'm doing these calculations outside of Moco that is the big difference?

[Ross Miller]

Here are CC velocities from study.analyze() and finite differences for a predictive simulation (this one here: https://twitter.com/rosshm16/status/1320793141352550412). There are still some small disagreements but not the large spikes from Russell's original post.

Screen Shot 2020-11-12 at 5.43.09 PM.png (715.84 KiB) Viewed 850 times

I wouldn't worry overly much about differences on the order of 0.1 m/s, I think the finite differencing errors can be about that large assuming ~100-Hz data, but it would be good to get to the bottom of the very large differences.

I'd be surprised if velocities that big produced realistic metabolic costs, but hard to say, will depend on the model.

Russell, it looks like all the big Vcc spikes except TA and EHL are happening when activation is nearly zero. Do you have the lower bound on activation set to literally zero? If so and if the optimizer is actually picking zero or something very close to zero, that could be causing problems when computing Vcc. Or similarly, if the lower bound is zero, I think IPOPT can actually then converge with very small negative activations (within constraint tolerance) which could also cause weirdness when study.analyze() computes Vcc. I've been using 0.001 as the lower bound on activation and excitation.

Ross

[Nicholas Bianco]

Hi Russell,

Thanks for the quick response and super clear plots. Low activation seems a likely culprit. I agree with Ross suggestion on modifying the lower bound on activation. Fiber velocity is calculated differently depending on whether or not tendon compliance is enabled, and if enabled, which mode (implicit or explicit). What settings are you using?

-Nick

[Russell Johnson]

Hi Nick and Ross,

My above results were with a minimum activation bound set at 0.001, which seems to line up with what Ross has done. When looking through the activations, it looks like the solution is staying within the minimum bounds and not producing muscle activations less than 0.001 or negative values. Out of curiosity I ran a new optimization last night with the minimum activation bound at 0.01 and here are my muscle fiber velocities below. The results seem to be more consistent with Ross's figures comparing the derivative of length and velocities, but still a few smaller peaks here and there.

Muscle Fiber Velocity Fig_MusActBound0p01.jpg (173.3 KiB) Viewed 799 times

I have been using enabled tendon compliance in implicit mode.

For my metabolic cost results, I have plotted the COT across a few speeds here from the results with the minimum muscle activation bounds set to 0.001 (my original settings). The model is 2-D so I suspect the 1.50 m/s condition might be a little high due to the constraints of 2-D, but otherwise I think it's pretty decent. But again, I'm calculating this outside of Moco using the Umberger probe, so maybe the probe is calculating velocity in a different manner that doesn't result in those big spikes seen with the analyze() function?

metabolic_COT_speeds.jpg (15.72 KiB) Viewed 799 times

When calculating the metabolic cost with the increased muscle activation bound for the 1.00 m/s speed, the COT was ~4.1 J/kg/m which is a little higher than the result shown above, but consistent with the fact that the minimum activation bound is now 10x greater.

I could be off, but perhaps another potential culprit is my set_implicit_auxiliary_derivatives_weight setting, which is currently set to 0.001? I noticed Ross's most recent push of his UMocoD has that setting to 0.00001. Could that play a critical role here?

Thanks for your help with this.
Russell

[Nicholas Bianco]

Hi Russell,

Thanks for the update, this is helpful. Activation bounds around 0.01 are pretty typical; I'm actually a little surprised the COT changed as much as it did. I'm also surprised that fiber velocity is so sensitive to activation considering that you're using the implicit form of tendon compliance (which should avoid low activation singularities).

The derivative of normalized tendon force, which is minimized when set_minimize_implicit_auxiliary_derivatives is enabled, should have a larger impact on fiber velocity. You could be right: increasing the auxiliary derivatives weight should reduce the tendon velocity, which could mean the fiber has to take up a larger portion of the whole muscle-tendon velocity. I'd be interested to see how the fiber velocity results change when changing that weight.

-Nick

[Ross Miller]

My model in 3-D also has 84 muscles vs. I think 20 in Russell's model in 2-D, and Moco sums the dF/dt term over muscles, so the weight on my dF/dt term is not dramatically lower, although it's still about an order of magnitude lower [(0.001*20)/(0.00001*84) = ~24x lower].

The weight on dF/dt I think will definitely affect the fiber velocities, but I wouldn't expect it to cause differences in study.analyze() vs. d(Lcc)/dt that differ by 4-6 orders of magnitude.

When my simulations converge, I'll often get warnings like this from Moco:

Code: Select all

DeGrooteFregly2016Muscle::computeForce, muscle semimem_l force < 0 at time = 0.113333
DeGrooteFregly2016Muscle::computeForce, muscle semimem_l force < 0 at time = 0.116667
DeGrooteFregly2016Muscle::computeForce, muscle bflh_r force < 0 at time = 0.593333
DeGrooteFregly2016Muscle::computeForce, muscle semimem_r force < 0 at time = 0.613333
DeGrooteFregly2016Muscle::computeForce, muscle semimem_r force < 0 at time = 0.616667

I think it's from the fiber damping as Russell suggested, producing very small negative forces. They tend to happen around toe-off as you can see from the timing here (this was a walking simulation with Tstride = 1.0 s). I've not checked to see if the timing of those negative forces coincides with wonky-looking Vcc's, but that might be worth cross-referencing.

Ross

[Nicholas Bianco]

Russell's suggestion is correct, you'll occasionally see negative muscle forces due to damping when activation is low and fiber velocity is high.

In our experience, these negative forces are usually very low (~1/1000 of max isometric force).