Does OpenSim report state equations?

[Ton van den Bogert]

Small correction to my previous post: I mistakenly referred to Ajay rather than Chris.

[Christopher Dembia]

Thanks for all of that great advice, Ton. I'm also curious about time spent in OpenSim versus IPOPT. Based on some stuff I've done, I'd guess OpenSim is the bottleneck.

As for analytical versus numerical derivatives, getting Simbody to produce analytical derivatives is definitely an enormous undertaking. I'll continue to talk with Sherm about its feasibility, but I agree that in the mean time, our focus should be on finite differences. Perhaps a hybrid is possible in the interim, with analytical derivatives for muscles and finite differences otherwise?

As for muscle models, my preference is to create a new muscle model in OpenSim that is specifically tuned to play well with direct collocation rather build upon the ones that already exist in OpenSim. Once we have a muscle model in OpenSim that is tuned for direct collocation, we can examine making the existing muscle models work better with direct collocation.

[Michael Sherman]

As an example, assume that we are using finite difference method to estimate the Jacobian matrix and assume we want to calculate the derivative of x_dot with respect to one of the controls (e.g. excitation of a particular muscle).

Hi, Sina. I don't understand your example. If the control causes the muscle to generate forces those will in general affect many or all of the accelerations in x_dot. Those are coupled so can't be calculated individually. Or were you referring only to the Jacobian of the muscle force w.r.t. the control?

Regards,
Sherm

[Sina Porsa]

If the control causes the muscle to generate forces those will in general affect many or all of the accelerations in x_dot. Those are coupled so can't be calculated individually. Or were you referring only to the Jacobian of the muscle force w.r.t. the control?

It's because of the way that Thelen2003 muscle calculates muscle force (i'm not sure about other muscle models), in thelen 2003 model, the muscle force is not directly a function of excitation or even activation. This model calculates the muscle force (tendon force) based on generalized coordinates and muscle fiber lenghts. knowing the joint angles, total length on musculotendon actuator can be calculated. as the fiber length is one of the states, the length of the tendon can be calculated as : Lt = Lmt-Lm. Thelen2003 uses this tendon length to calculate to force.
Sina

[Ross Miller]

If someone could get this working before next winter so I don't have to renew my Autolev license that would be great =)

My computer programming ability when it comes to creating general-purpose applications in environments like C++ is minimal, so I'm not sure how much direct help I would be in expanding the existing capabilities of OpenSim, but I'm happy to contribute in any other ways that are useful or helpful.

Re: bottlenecks, in my experience with DC, the time spent in the optimizer itself is the bottleneck (I think from estimating the "Hessian"), but that is using an optimizer that's very good at taking advantage of sparsity information (IPOPT) and with a fully symbolic Jacobian (from Autolev). I've not tried it with a symbolic Hessian (maybe Ton has?) just because it's a hassle to form it and get it Matlab-ready unless the model is very simple, but that could possibly provide another dramatic speed-up.

I did not have good results when using fmincon with DC. It worked eventually and gave reasonable results but was very slow. It also required the full version of Matlab, the student version had a memory limit that was too restrictive (a nine-DoF / 18-muscle gait model would not run on a 51-node grid). I could never figure out how to get it to take advantage of sparsity in the nonlinear constraints. This was 2012 or so, I've used IPOPT exclusively since then and it works very well (and is free). If fmincon is different now then this could be an attractive option. Getting IPOPT running on a specific machine/OS/compiler combo can be tricky.

Ross

[Ton van den Bogert]

IPOPT estimates the Hessian of the Lagrangian, but that's essentially "free" if you use the "limited-memory" option. This uses BFGS updates that only require information from successive evaluations of the objective gradient and the constraint jacobian. Finite differences for the Hessian is probably not going to be feasible for systems with many state variables. Let's first get efficient first derivatives...

Analytical Hessian requires second derivatives of the system dynamics, and I fear that Autolev would choke on that for models of moderate size. And it would require NxNxN matrices which Autolev does not do, so some data structure trickery is needed. Not something I want to try soon...

SNOPT, on the other hand, only requires the hessian of the objective function which is sometimes quite easy to compute. For instance if the objective only contains controls and tracking of state variables. The Hessian of the objective, in that case, is actually constant and I have sometimes used that advantage. We used to have TOMLAB/SNOPT but it's quite expensive. We then switched to the version from Stanford Business Software but it is harder to use and not much luck yet.

In my experience, time spent in the model and in the optimizer is similar, with a setup similar to Ross's, but it can vary depending on how hard the optimization problem is. With Opensim, I doubt that IPOPT will be the bottleneck.

Thanks for sharing your experience with fmincon. I have tried it every 5 years or so since the optimization toolbox first appeared in the early 1990s and was always disappointed. Maybe it's better now.

[Brad Humphreys]

Maybe this should be a separate thread, but I wanted to add a +1 to Sina's suggestion on computeStateVariableDerivatives. I have read his and Brian's great papers and believe our formulations are similar but I wanted to share mine. My formulation (and constraint jacobian) has as "inputs": joint, angles, muscle fiber length, and activation. "Outputs" are: joint accelerations and muscle velocity from computeStateVariableDerivatives (and compared to finite differences). Recently I have tried to understand how/if excitation plays a roll in computeStateVariableDerivatives. A low overhead API function would be great and obviously an implicit capability would be wonderful.

Brad

[Brian Umberger]

Hi All,

Sorry for not responding sooner. I thought if there was more activity in this thread I would be alerted by SimTK via email (maybe that only happens if you initiate the thread), thus I didn't realize until yesterday how much further the conversation had progressed.

Chris: I'm happy to help out as best I can. I'm not a C++ programmer so my ability to make low-level contributions will be severely limited. I also don't have first-hand experience with the implicit approach like Ton and Ross (and maybe others). That said, I will be excited to see this move forward and will help where I can.

Ton: Here is the timing output from IPOPT for the lower-limb positioning simulations in our recent paper. This should answer your question about where the time is spent. As you can see, most of it is spent in OpenSim getting the information needed to compute derivative via finite differences.

25 Node case (825 unknowns, 615 constraints)
Total CPU secs in IPOPT (w/o function evaluations) = 5.477
Total CPU secs in NLP function evaluations = 1158.557

200 Node case (6600 unknowns, 4815 constraints)
Total CPU secs in IPOPT (w/o function evaluations) = 27.341
Total CPU secs in NLP function evaluations = 7800.382

In our current walking simulations, the effect even is greater if we track GRFs, because of the additional overhead associated with getting the magnitudes of the GRFs from OpenSim. For predictive walking simulations the ratios are similar to above.

Regarding efficiently calculating sparse Jacobians via finite differences, I agree with Ton that does not necessarily need to be an OpenSim-side issue. There are several algorithms for doing this effectively. We started with the Curtis algorithm cited in our paper, but most recently we have been using the built-in numjac.m function in Matlab. We determine the full sparsity pattern of the Jacobian (ahead of time) by evaluating the dense Jacobian for a variety of values of X, and then flag elements that are zero in all cases. This sparsity pattern is provided to IPOPT at run-time. Because of the small number of non-zero elements and the nature of the sparsity pattern for these types of problems, numjac is able to greatly reduce the number of calls to OpenSim required to evaluate the full, sparse Jacobian. In a typical example, a constraints Jacobian with 5000 columns can be determined with only about 50 calls to OpenSim. That is a glass-half-full, glass-half-empty situation: it would be 5000+1 calls in the dense, finite difference case, but only 1 call in the analytical case. For the time being, we are happy with the speed-up this gives us over direct shooting.

I will close by echoing what other have said that it is great to see the level of interest in this approach both within and outside of the core OpenSim team. I'll also check back on this thread more frequently.

Brian

[Ton van den Bogert]

Brian,

At the top of the page is a "subscribe topic" button to get e-mail notifications. When it has a red x next to it, that does not mean that you are already subscribed... I made that mistake before. Click the button until it says "unsubscribe topic".

Thanks for posting those timing results.

Here are my timings for a 50-node case (2d model, half gait cycle, minimize squared controls, random initial guess).

Code: Select all

Number of nonzeros in equality constraint Jacobian...:    32767
Total number of variables............................:     4182
Total number of equality constraints.................:     3382
Number of objective function evaluations             = 1408
Number of objective gradient evaluations             = 768
Number of equality constraint evaluations            = 1408
Number of equality constraint Jacobian evaluations   = 831
Total CPU secs in IPOPT (w/o function evaluations)   =    193.659
Total CPU secs in NLP function evaluations           =     44.090

So for me right now, IPOPT dominates the computation time. This is with implicit dynamics and analytical Jacobians. This indicates what may be possible with Opensim when it has that capability.

Ton

[Ross Miller]

A couple related comments that I've found through trial and error have a big effect on the CPU time in IPOPT/DC:

- IPOPT takes a long time to converge (lots of iterations) with constraints that include any function that decays to an asymptote (I'm not sure of the technical term for these functions, sorry; sigmoid?). I discovered this from using Hatze's force-velocity equation, which is nice because it's one single equation regardless of which FV "limb" you are on, but uses hyperbolic tangent and decays to F=0 as V->-Inf (and to F=Fecc and V->Inf). I think this could probably be minimized with careful/thoughtful setup of the problem and convergence criteria, but it's easier to just use functions that don't have these characteristics (e.g. have an FV curve that can go below zero or above Fecc; I think Ton mentioned this earlier).

- I'm not sure if this is strictly necessary, but it definitely helps if the optimizer can always "feel" the generalized forces being applied to the model even when the modeler might intuitively want these forces to be zero. E.g. for a deformable discrete-element contact model, it helps if the force is still a function of the contact point's kinematics even when contact is inactive. This could also appear in the passive joint stiffness (e.g. if it's literally zero across a certain chunk of the ROM then becomes non-zero) or the SEC force-extension relationship.

Ross