Thelen2003Muscle::initMuscleState updated. Algorithm changes: bisecti…

[mjhmilla]

Fixes issue #<issue_number>

Thomas G's initialization problems (discussed in the meeting 3 March 2020)

Brief summary of changes

1. Algorithm change : bisection used for good initial conditions for the Newton method.
2. Algorithm change : A text book Newton method (not the previous version with the line search) polishes the root off.
3. Numerical change: Problems were being caused earlier when the Newton method encountered a zero gradient: either due to a slack tendon (this case was handled before) or when fv goes to zero (an un handled case). With Thelen's force-velocity curve fv can go to zero quite easily if activation is low. To prevent the gradients from going to zero (and the ensuing Newton method problems) I've constrained the tendon strain and fv to be 10 x larger than the solution tolerance. This is very small, but large enough to prevent these problems.

Testing I've completed

Ran all of the tests in release and debug mode.

Looking for feedback on...

1. It would be great if Thomas G tested this out on the case that was failing.
2. The method: the variable names do deviate from the camel_case convention. Shall I update them them?

CHANGELOG.md (choose one)

• no need to update because: this is deep in the bowels of the code.

---

This change is 

[chrisdembia]

I'm wary of making such big changes to the algorithm; is there a risk that we create other issues that didn't exist before? What about fixing only the force-velocity issue?

Reviewable status: 0 of 1 files reviewed, all discussions resolved

[mjhmilla]

Hey Chris,

I see above travis-ci and appveyor had some complaints: is this something I need to worry about? I looked through the details and didn't find anything obvious.

As for the big changes to the algorithm: in principle bisection+Newton should be more robust than the variable-step Newton method I had in there. This is actually what I've been using since my time at Stanford (and perhaps before I left). I originally used a hand-made variable step Newton method (e.g. this is something I made up!) in an effort to save a bit of computation. In hind-sight this was a bad decision. That's why I updated the root solving method.

If you wanted to just fix the force-velocity issue, I think that is straight forward: the "auto velocityFunc" needs to be updated. That's easy to do. In any case, it would be valuable to have the specific case where initMuscleState was failing: then I could verify (in debug mode) how it was failing, we can verify that the fix works, and include that as a test case. Without this new test case you actually cannot be sure that anything we've done is an improvement.

So do you have the test case?

[aseth1]

@mjhmilla @chrisdembia I am not fond of the idea of "fixing" the algorithm when it was not clear if this was even the issue. For example, I used Scott U's case and I was not seeing fv (if you mean the force-velocity multiplier) going to zero. Let me know if you want me to create a separate PR for a test case or add it in here. I should be able to do that this weekend.

[mjhmilla]

@aseth1 : Please do add the test case that was failing with the earlier version of initMuscleState in Thelen2003Muscle. I don't think we will reach a fix that everybody is happy with until this test case is in place.

[aseth1]

@mjhmilla see #2687 for the related test case. Feel free to merge and change with this PR as you see fit.

[mjhmilla]

@aseth1 Thanks! I'll try to sneak some work on this in the coming week.

[mjhmilla]

I have not forgotten about this, but my work on this will be delayed: I've been dropped into a project with paper deadline in just over 2 weeks. I'm currently working 7 days a week and I'm still not certain that we will make the deadline.