Muscle equilibrium bug thelen2003

[aseth1]

Test case for issue #2686

Brief summary of changes

A test case was added that is a stripped down (one muscle) model and conditions that reproduce the muscle equilibrium solve failure encountered in #2686.

Testing I've completed

This test case fails as the activation approaches the activation value that caused the original failure.

Looking for feedback on...

@mjhmilla please use this test case to investigate and verify that #2684 is the correct solution

CHANGELOG.md (choose one)

• no need to update because this is a test case for existing functionality

---

This change is 

[chrisdembia]

Approaches 0.08

[aseth1]

Thank you for correcting the typo.

[aseth1]

Thank you @chrisdembia for the review and approval.

[mjhmilla]

I wrote way too much earlier - a shorter version is below. I am not issuing another pull request until either @chrisdembia, @carmichaelong or @aseth responds and approves the idea of the minimal fix, or asks for something else: I do not have a lot of time these days to work on code that gets thrown out.

Problem

The initialization method is using Newton's method and the Thelen2003Muscle has C1 and C2 disconinuties: the test case Ajay made happens to have a its solution at lceN=1 and so that C1 discontinuity is causing the Newton method to jump from branch to branch.

Minimal Fix

I can add a bisection method ahead of the Newton method and return the better of the two solutions. In the worst case the user will get an initialization with a tolerance of fiso/2^15, in the best case the Newton method will polish it up to 1e-8 N or better.

Expectations: fixing bugs

There are two other problems I noticed while I was looking at the initialization method: if fv goes to zero the Newton method fails, and there is a failure case with weird Newton + line-search method that is in there. I see now that I cannot just fix things, so, in the future if I see something that really should be fixed should I build a test case first?

In principle the minimal fix will cover these problems: in the worse case the rougher bisection-tolerance solution will be returned when a much better Newton-tolerance could have been. That's not so bad.

Notes: test case

Before Ajay's test case gets pulled into master the CMake has to be updated to copy the osim file over to the build directory: on my machine I had to do this manually, which isn't going to work for the average user. The alternative is that I hand craft a test case that has an initialization right at the sweet spot of lceN = 1.

Related Future Work

Is there interest in an implicit version of the 2012 Equilibrium model? If yes, this is something I could start to work on. It should require some modest changes to the curves (so that the gradients don't go to zero) and one additional function.

[chrisdembia]

The dev team discussed what to do about these changes and we would be a lot more comfortable with larger changes to the equilibrium algorithm if there's a more thorough test across activation levels and MTU lengths. We have started prototyping such a test for the DeGrooteFregly2016Muscle in Moco:

https://github.com/opensim-org/opensim-moco/blob/4e77683b6c8e1dbaa5937bc994722f3de92e3f05/Moco/Tests/testDeGrooteFregly2016Muscle.cpp#L1176-L1247

Such a test should be straightforward to implement and can reduce uncertainty when changing the algorithm. I'm happy to work on the test with you.

I see now that I cannot just fix things, so, in the future if I see something that really should be fixed should I build a test case first?

We have a pull request review process outlined here. It would be great to build a test case first and discuss the proposed changes with the development team before opening a pull request. Naturally, this depends on the scope of the changes.

Related Future Work

That's a great idea! I created a separate issue #2711 to discuss this.

[mjhmilla]

Thanks for the response - I'll set aside some time this week to work on:

• a more generic version of the initialization test to go in testMuscles
• the minimal fix to the Thelen2003Muscle code

[chrisdembia]

Sounds great :)