Blips in simulated muscle variables

[aymanhab]

Reported by Ton on the user forum at
https://simtk.org/plugins/phpBB/viewtopicPhpbb.php?f=91&t=15001&p=0&start=0&view=&sid=05adfb8d846e5cd0e9f3b37ed5d3e37d
Preliminary investigation suggests this is related to the muscle model rather than faulty reporting. Using rigid-tendon muscles shows no such blips and increasing damping to 10 or more eliminates the blips. The properties of the muscle model have little documentation of default values or acceptable ranges of values. This could be a bug or just picking bad values to initialize the muscle model in this example. The same issue is reported for Millard and Fregly muscle models in case that helps.

[aymanhab]

@tkuchida Any idea what maybe wrong here or how to fix/adjust the example so that these curves are smooth? Thanks

[tkuchida]

@aymanhab I don't have much time to dig into this; however, the Simbody kerplosion 💥 suggests something more serious than a non-smooth muscle force (even if that is, indeed, an issue). Looking at build_and_simulate_simple_arm.m, I noticed that the moments and products of inertia are zero (Inertia(0) on lines 38 and 39). I would first rule this out as the underlying cause (a zero inertia could be causing infinite or at least high accelerations). The example in README.md was updated in #3186 so the inertia matrix makes sense… perhaps try that model instead? 🐧

[aymanhab]

Thanks for the quick response @tkuchida This was created using the latest README which has the corrected inertias. I made a PR to include the TableReporter to demonstrate/reproduce.
The plot is shown below:

[aymanhab]

@mjhmilla Is there a chance you can explain what's wrong with the setup? Thanks

[mjhmilla]

I'm just about to board a plane for a holiday, so I'm out for next couple of weeks. If these are the facts:

• blips present on Millard & Fregly elastic tendon muscle models
• blips absent on a rigid tendon model

then my guess is that something in the model could be putting noise at the position level. Since an elastic tendon muscle model has fiber length as a state, any noise at the position level will result in large force transients when the tendon is stretched. In contrast, a rigid tendon model force output will not be greatly affected by small amounts of position level noise. If I'm correct, then in both cases there will be small discontinuous changes in path length that show up in all simulations at the position level of affected muscles (with wrapping cylinders, perhaps). But what could be causing this type of error?

@tkuchida mentioned that there are some bodies with zero inertia, this could be a place to look. Ton mentioned that there are wrapping cylinders used in the model, this could be another place to look. My suspicion is that the muscle wrapping routines are periodically not meeting the desired tolerance and are injecting position errors into the muscle path length. Because path velocity is usually solved separately, these defects might be large at the position level and very small at the velocity level. In contrast, if this is a problem of small masses, the defects would be continuous and would show up in the path length and path velocity outputs.

If I had the time, I would look at the path length and velocity profiles of the affected muscles near the occurrences of the blips. After that I'd either update the inertia and/or wrapping (perhaps just eliminating the cylinder wrapping and associated muscles) to confirm the problem. If this is still unresolved when I get back from holiday I'll take a look at it.

[aymanhab]

Thanks @mjhmilla Some clarifications:

1. Model has no wrapping (README file of this repo).
2. Mass is non-zero and inertia values are set from latest code updated by @tkuchida so non-zero.
3. Changing muscle model to rigidTendon shows smooth states, length and no blips or noise whatsoever AFAICT
4. Changing damping affects the blips significantly e.g. he're the curve for damping = 10
   

I understand you're on vacation, just wanted to keep everybody in the loop regarding status. Thank you

[pepbos]

Hi! @aseth1 , @adamkewley and I have been investigating these Blips. It appears that the muscles have all the traits of a numerically stiff problem.

Numerical stiffness is a nasty problem, which lacks a formal definition. It is actually directly defined in relation to explicit integrators. However, it is not necessarily an indication that there is something wrong with either the integrator, nor the model. Some differential equations are unfortunately numerically stiff, in which case perfectly fine explicit integrators are not performing well.

You should be able to remove the blips by increasing the accuracy: A trademark of numerical stiffness is that this will hardly affect the integrator stepsize. Increasing the fiber-damping will also work, because it makes the problem less stiff, but you will also get a different response.

We are still further investigating this issue :)

[mjhmilla]

Thank you very much @pepbos , @adamkewley , and @aseth1 for digging into this in more detail. There's a good chance that numerical stiffness is the cause: muscle modes are numerically stiff. From here it would be great to confirm that numerical stiffness is the source of the this problem: do the blips go away if one of OpenSim's implicit integrators is used (cpodes)? I see that cpodes is still in Simbody, so this test should be possible.

If numerical stiffness is confirmed, the next step is to confirm the source:

• If pennation angles are approaching 90 degrees this will cause stiffness in the model, but this can be fixed. I ran out of time when I wrote the first version, but I wish I had been able to implement a suggestion of Ton's: to make the state variable lce*cos(alpha) rather than lce. This removes the considerable nonlinearities that the pennation model injects into d/dt lce which can reliably produce a numerically stiff problem. That said, its hard for me to imagine that this is the issue for a bicep.
• If instead the multipliers of the curves are rapidly changing, then things are a bit more difficult. I've just looked at the source code to confirm that the Bezier curves are evaluated to a very tight tolerance: EPS*100. As such, if the multipliers are rapidly moving around, then its likely because of the shape of the curve. At this point there are only two options: switch to an implicit integrator (or use a smaller maximum step size) or change the shape of the offending curve.

If numerical stiffness is not confirmed:

• That would be interesting. At this point in time I don't see any obvious reasons for the noisy output given the problem.

I'm drowning in deadlines right now. I hope I can steal a little time to run this myself once I get a submission out the door in the next week or so.

[mjhmilla]

Just to add, from @aymanhab first message the statement appears:

"Using rigid-tendon muscles shows no such blips and increasing damping to 10 or more eliminates the blips. The properties of the muscle model have little documentation of default values or acceptable ranges of values. "

The damping parameter should be small as possible without introducing numerical problems during integration. Why should the damping parameter be small? If you look at the equilibrium equation, it should be clear: the damping element will act as a drag on the CE and slows it down. When the damping parameter is small (say 0.1, or 0.01) the simulated force-velocity curve of the model is very close to force-velocity Bezier curve. If a large damping value is used (like 10, which is ENORMOUS!) the CE won't be anywhere close to following the force-velocity Bezier curve: it will be much slower.

This also points to another numerical improvement that can be made to the model: make a new force-velocity curve that compensates for the damping term, at least at maximum activation. This would allow larger damping coefficients to be used without killing the accuracy of the model. The caveat is that this compensation would only be perfect at maximum activation, unless the numerical damping term is also made to scale with activation.

[adamkewley]

Just as an update @mjhmilla , we (@aseth1 , @pepbos , and me) are currently looking into this in parallel with other changes (e.g. perf-optimizing the underlying bezier curve implementation: #3481).

Apologies if we're a bit slow on responses (I tend to scour GitHub in very sporatic waves) but, from what I understand, we are currently doing the basic background work (e.g. reading your original paper, but also figuring out why models like Rajagopal2015 use a damping parameter of 0.01, rather than 0.1).

It might be an idea for us to shoot you an email or at least set up some kind of IM system (teams/slack), or a call, once we start to have concrete ideas about what we think might be a good idea!

[mjhmilla]

Dear @adamkewley , @aseth1, @pepbos,

That sounds great. The paper isn't going to tell you much about the Bezier curve implementation, for that you're going to need to look at SegmentedQuinticBezierToolkit's doxygen and associated code. I'd be quite happy to chat muscle code with you when you're ready. In regards to the two issues you described. I assume that the priorities are #1 reduce the chatter, #2 make the curves faster.

1. I'm quite sure we can remove the chattering Ton observed using these two changes:
a. Adjust the force-velocity curve to permit a larger (constant) damping parameter to be used without degrading accuracy of the model when replicating maximum activation shortening experiments. It would take some work to set this up, but this shouldn't be more than a few days work.
b. If 1a works, we could proceed to decide how to handle shortening during sub-maximal activation: as 1a is described the model will follow the force velocity curve at maximum activation, but at sub-maximal activation the damping parameter will be too big, and the contractile element will be slower.
c. Unrelated: the force velocity curve constructor should fit the curviness parameter such that the Bezier curve follows Hill's hyperbola as close as possible. I've written this code elsewhere, but haven't yet moved it into OpenSim.
d. I also suspect giving the tendon a damping model would help as well. It turns out that tendons are lightly damped - which is a fact I've had to include in a model I've made that's under review.

2. Curve performance improvements:
a. Quick change: I strongly suspect that calcU can be improved to reduce simulation time (near line 788 of SegmentedQuinticBezierToolkit.cpp). This is the only function that uses iteration in the Bezier spline code, and it uses a really expensive hint to get close to the solution: a large SimTK::Spine of the curve. This is expensive because the spline has on the order of 100's of points, and waiting for memory to be fetched is time consuming. The spline should be replaced with the bisection method to get close to the solution, and then Newton's method to polish it off. We'd need to benchmark this pretty rigourously to make sure the memory vs. computation trade-off is worth it. I suspect it will be.
b. Medium-time-consuming change: As far as I can tell, nobody (?) makes use of the C2 continuity of these curves. The C2 quintic Bezier curves can be replaced with C1 cubic curves, as I have done with an implementation in LS-DYNA. With cubic curves there is no need to use numerical iteration, and so these are much faster. The consequence is that subdivision is required to make the cubic curves actually follow the path that you want, but this is pretty trivial.
c. Time consuming change: replace the entire set of curves with a series solution of hyperbolic functions. I've found that I can describe the derivative function of all of the curves using a series solution of tanh functions. Lucky for us tanh has a symbolic integral. The result from this conversion would be a new curve library that is continuous, requires no iteration (beyond what's done to evaluate hyperbolic functions) and yet rich enough to faithfully follow whatever curve you'd like. The downside is that iteration would be needed to construct these things, but evaluation would be cheap ... under the assumption that evaluating hyperbolic functions are cheap.

[pepbos]

@mjhmilla hi,

Just to reply on how we verified that the muscle is very likely numerically stiff:

1. The result is noisy despite a smooth solution (the response is smooth for very tight accuracy, like 1e-10),
2. Many small steps despite smooth solution,
3. Integrator time step size is stability limited:
   3.1 changing the accuracy hardly affects the step size,
   3.2 changing the fiberdamping does affect the step-size (and the dynamics, and the numerical stiffness)
4. Implicit integrator (Sundials CVODE) gives smooth result, and is able to take much larger time-steps.

Some of the solutions we are currently exploring for this problem:
a. Effect of fiberdamping,
b. Choice of state (fiber length, or tendon force)
c. Implicit vs explicit integrators
Each of these appear to have a significant impact on the response.

The fiberdamping has a major effect on the numerical stiffness. We would love to investigate the range of permissable fiberdampings, and how to adjust the curves. As a first step we got our hands on the original data used in 2013, but it would be great if we can sent you an email or something once we have a plan of attack.

Curve performance improvements is something we can also pursue. Point 2a makes a lot of sense.
Interestingly, we are also experimenting with a specific implicit integrator that would make use of C2 continuity! But this is still pending results.

[mjhmilla]

Hello @pepbos,

You've checked all of the numerical boxes for numerically stiff problem. The last time I talked to Sherm about this very issue he mentioned that the numerical integration experts just use 4: if an implicit integrator is much faster than an explicit integrator the problem is numerically stiff.

That reminds me: what is the ratio of the tendon slack length to optimal fiber length? Numerical stiffness increases as tendonSlackLength/optimalFiberLength becomes small. According to the build_and_simulate_arm_modified.txt file from Ton's email the optimal fiber length is 60 cm long, and the tendon is 55 cm long. So ...

1. When tendonSlackLength/optimalFiberLength is below 1, there's not much of a loss in accuracy to just switch to a rigid tendon. This is already an option in the model.
2. If the tendonSlackLength/optimalFiberLength is above 1, the tendon damping model I mentioned may also be helpful.
3. The pennation model can also make for a numerically stiff problem. But as I see from the file, there is no pennation.

Ok onto the specifics of your message:

The fiberdamping has a major effect on the numerical stiffness. We would love to investigate the range of permissable fiberdampings, and how to adjust the curves. As a first step we got our hands on the original data used in 2013, but it would be great if we can sent you an email or something once we have a plan of attack.

I could send you a Matlab prototype to generate the force-velocity curve that I'm thinking of, given a specific fiber damping value. To really test this properly, we'll need to simulate a series of ramp shortening experiments to make sure that the force-velocity curve produced when replicating a force-velocity experiment actually matches what we want.

Curve performance improvements is something we can also pursue. Point 2a makes a lot of sense.
Interestingly, we are also experimenting with a specific implicit integrator that would make use of C2 continuity! But this is still pending results.

In that case either option 2a or 2c would be the best long term solution.

[pepbos]

Hi @mjhmilla ,

The ratio of tendonSlackLength to optimalFiberLength is indeed just below one (~0.92) in this case, and switching to a rigid tendon seems to give a similar response, and greatly reduces the numerical stiffness.

However, looking at your paper from 2013 ("Flexing computational muscle"), I concluded that it was found that a rigid-tendon gives a larger modeling error, and would be a less accurate representation of a real muscle. In that paper, as well as others, I gathered that the muscle dynamics behave as a numerically stiff system. I therefore assumed that the rigid-tendon muscle is sort of checked off the list as a solution, because it no longer matched the experiments accurately.

1. Do I understand correctly that certain muscles are just as accurately represented using a rigid-tendon? (due to their specific parameters differing from those used in the experiment)
2. Does it make sense to define when a rigid tendon should be used? Perhaps this could then be used to take some kind of action (give warning, or switch model etc.).

Still, I assume that some muscles will require a non-rigid-tendon model, for which the numerical stiffness poses a challenge. In that regard, getting back to your other points:

If the tendonSlackLength/optimalFiberLength is above 1, the tendon damping model I mentioned may also be helpful.

Would it be possible for us to access this damped tendon model? It would be interesting to see how big the effect is on the numerical stiffness and response.

The pennation model can also make for a numerically stiff problem. But as I see from the file, there is no pennation.

The pennation angle dynamics that I see being used is given by $\dot{\alpha} = - v^M\sin\alpha / \left( l^M \cos \alpha\right) $. Not sure if this is the pennation model you refer to?

I could send you a Matlab prototype to generate the force-velocity curve that I'm thinking of, given a specific fiber damping value. To really test this properly, we'll need to simulate a series of ramp shortening experiments to make sure that the force-velocity curve produced when replicating a force-velocity experiment actually matches what we want.

The effect of fiber damping seems an interesting avenue to explore: We would be very interested in a matlab script outlining the curve generation!