Does OpenSim report state equations?

[Sina Porsa]

Thankyou Ajay.

[Ton van den Bogert]

This thread has been quiet for a while, but I would like to follow up and advocate for better support of the type of work that Sina is doing. We had a collaborative R01 grant with Simbios during 2006-2009 and eventually had to give up on using Opensim with direct collocation. Sina has recently managed to make it work with Opensim, which is a great accomplishment but computation time, although orders of magnitude faster than shooting, is still slow and local optima are still a problem.

In Sina's original posting, he asked if the dynamics could be provided in implicit form:

f(x,u,dx/dt)=0

As I argued in my 2011 paper, this would have many advantages. Numerical conditioning would be improved, to the extent that you can even have a system with a non-invertible mass matrix. Muscles can have zero activation. Because of the improved numerical conditioning, local minima are probably less likely. In the implicit formulation, it is much easier to avoid singularities and discontinuities which wreak havoc on optimizations. And scaling of the optimization problem is much improved in when implicit dynamics is used. When dx/dt is solved from the dynamics, the magnitudes of these state derivatives can vary enormously. For instance, an ankle angular acceleration can get very large because the inertia of the foot is low. In the implicit formulation, f is the difference between the required joint moment for the given x and dx/dt, and the joint moment generated by muscles and other elements. For numerical scaling, this can easily be normalized to body weight or a maximum isometric joint moment, and then you can expect values between -1 and 1.

An API function that returns f should not be too difficult to make because all elements already exist in Opensim. The f for the multibody dynamics would be

f = M*qdd + C(q,qd) - R*F

where M is the mass matrix, R are the muscle moment arms, F are the muscle forces, and C is everything else (gravity, passive elements, coriolis, centrifugal). q,qd,qdd are the generalized coordinates and their velocities and accelerations. I know that Opensim projects the multibody dynamics onto the constraint manifold, and this could complicate things. Maybe we can at least use this for systems without kinematic constraints, i.e. completely free generalized coordinates. In the dynamic arm simulator project, we replaced kinematic constraints (scapula-thorax contact) by very stiff elastic contact and that worked fine and has the additional advantage that contact can be restricted to compressive force only.

The f for muscle contraction is simply the force imbalance between the series elasticity and the other muscle elements:

f = Fsee(Lmt - Lce) - (Fce(Lce, dLce/dt, activatin) + Fpee(x)) * cos(pennation)

For numerical scaling, this can easily be normalized to the maximum isometric force Fmax to get a dimensionless number near the -1 to +1 range. Optimizers such as IPOPT are very sensitive to scaling. In the conventional explicit dynamics, dLce/dt is solved and this will go to infinity when the muscle is near its eccentric maximum force.

We are keen on using Opensim with collocation approaches for optimal control, but so far, it has been necessary to write our own musculoskeletal dynamics code. Not everyone has the knowledge or the resources to do this. If Opensim had this capability, suddenly, modern optimal control techniques become available to the entire Opensim community and can be applied to any model that already has been developed in Opensim.

Ton van den Bogert

[Michael Sherman]

Hi, Ton. I like this idea a lot. OpenSim already has an inverse dynamics solver that comes very close to the implicit formulation you want. It uses the Simbody API for the residual form of the multibody dynamics , which is exactly as you described it for unconstrained systems. There is also an API for constrained systems, although how best to use it is an open question!

I think the only missing piece would be the residual form of the muscle dynamics, which should be an easy thing to add. Ajay Seth and Chris Dembia likely would have a better feel for what would be required; I'll see if I can get them to join in.

The Jacobian of these equations is another story, though. We don't have a good way to calculate these analytically, so couldn't provide the exact Jacobian required by Rosenbrock methods. We would need to approximate numerically and solve the discretized equations using something more Newton's method-ish.

Assuming the approximated derivatives are tolerable, I think what we would need to get this done is a grad student, postdoc, or visiting scholar who would drive it and test it in a real application. Any ideas?

Best regards,
Sherm

[Ton van den Bogert]

Sherm,

Thanks for pointing us to the residual form of the inverse multibody dynamics. That gets us halfway there and will already help a lot.

I'm not so worried about the Jacobians yet, we can do that with finite differences. If you keep this in mind, some day we may have a dynamics engine that provides the Jacobians too. In our own code, we do this symbolically with Autolev. It works, but the expressions get terribly large. In the dynamic arm simulator, the equations of motions with derivatives are 13 MB of C expressions and take 3 days to compile when we turn on the compiler optimization. Execution time is probably also larger than necessary. I think it can be done much faster without symbolics and compilation, by propagating the derivatives (at run time) through all the steps that lead to Fr and Frstar in Kane's equations. All you need is the chain rule and some extra data structures. Sounds simple but probably isn't... BTW, you already have part of the Jacobian, the mass matrix which contains the derivatives of f with respect to the generalized accelerations. The finite differencing only needs to be done with respect to the generalized coordinates and velocities.

Looking forward to what Ajay and Chris say about the implicit form of muscle dynamics.

I have a graduate student (Brad Humphreys) who is working on this. He has already produced an optimal control solution for a 2-dof 8-muscle arm model with the Opensim API and direct collocation in Matlab. He could definitely drive this and test it on a real application. Brad works for NASA Glenn Research Center in Cleveland and they are very motivated to use Opensim for predictive simulations because they already have Opensim models which include detailed representations of exercise equipment that is used in the International Space Station.

Sina Porsa could also do it if he is interested.

Ton

[Christopher Dembia]

If Opensim had this capability, suddenly, modern optimal control techniques become available to the entire Opensim community and can be applied to any model that already has been developed in Opensim.

Ton, I care deeply about making direct collocation easy with OpenSim (and incorporating direct collocation within OpenSim) and so I am also very interested in the ability to easily get implicit dynamics out of OpenSim. This morning, Ajay, Sherm, and I briefly discussed how this might look. I think we should continue to discuss the design and implementation on GitHub (https://github.com/opensim-org/opensim-core/issues). Please keep me in the loop.

[Brian Umberger]

Hi All,

Ton and Ross Miller recently pointed me to this interesting discussion. My group is also intent on using collocation techniques to solve optimal control problems using OpenSim. As a first step in that direction, my former postdoc (Leng-Feng Lee) and I just published a paper that describes an approach for doing this using the Matlab scripting interface based on the explicit formulation and finite differences for the derivatives. The results are encouraging, but our current approach does not benefit from the potential advantages of the implicit formulation that Ton outlined in his 2011 paper. Here are the links to our paper and project page if you are interested:

https://peerj.com/articles/1638.pdf
https://simtk.org/home/directcolloc

I also just saw Sina's nice article on solving a vertical jumping problem using direct collocation vs. direct shooting. The examples we included in our paper were smaller in scale than the jumping problem that Sina solved. However, we are currently working on predictive walking simulations using our approach with a 2-D model and it takes about 4-5 hours to get a really nice result on a modest laptop. On the one hand this is amazing, as it is about 50x faster than direct shooting with coarse discretization of the controls, which yields a lower quality result (i.e., greater constraint violation). On the other hand, the performance could be even better, perhaps substantially so, if there was better support for this approach in OpenSim. Any steps that can be taken to enhance these capabilities in future releases of OpenSim will be greatly appreciated.

Brian

[Sina Porsa]

Hi All,
It’s exactly three years since when I started this thread and now it is converted to a very interesting discussion. I think we all do believe providing implicit formulation can highly help us to implement optimal control techniques.
Here I want to ask for a feature which can highly improve computational time, even with the current explicit formulation. As I mentioned in my paper, the sparsity of the Jacobian matrix cannot be fully exploited because of current version of the built in functions of OpenSim. OpenSim provides “computeStateVariableDerivatives” to calculate the f function (e.g to calculate x_dots). This built-in function accepts the current value of states (x) and controls (u) as input and calculates/reports the derivative of ALL of the states. (look at the for loop here:https://simtk.org/api_docs/opensim/api_ ... 3dc9c41d17)

If this function could calculate and report the x-dot values only for the states of interest, then the sparsity of the jacobiam matrix can be fully exploited and it will improve the computational times a lot. 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). In this case we just need to calculate the x_dot for the activation of the same muscle, but the current version of OpenSim calculates x_dot for all of the states.
I think it should be fairly easy to address this problem and it will be greatly appreciated by all of the researchers who are trying to use collocation with OpenSim.
Sina

[Ross Miller]

Hi all,

It is great to see this method becoming more popular, I think it has a lot of potential.

I don't have any experience with them personally so I can't speak to their robustness or ease-of-implementation, but there are well-established methods for numerically estimating the sparsity of a system's Jacobian (e.g. the Curtis et al. (1974) reference in Brian's PeerJ paper). In theory I think the software would only ever need to form this matrix once per model. It would only need to be re-computed if the model's dynamics are changing between iterations and some optimizers may not even allow that (I think IPOPT does not).

Something like this (automatic determination of sparsity matrix, output to Matlab) would be a really helpful addition to the OpenSim API.

Ross

[Christopher Dembia]

Brian, a few of us in our lab have read your recent paper. I really appreciate that you and Leng-Feng Lee shared code. I've also seen Sina's recent paper. Both papers are very exciting.

This thread has been exciting for me. I want OpenSim to be a platform for direct collocation. Not just for direct collocation to be possible with OpenSim, but for it to be easy and fast with OpenSim. I'd like direct collocation to be first-class functionality of OpenSim. I plan on taking many of the steps necessary to realize this.

This task is big, though. Ton, Sina, Brian, Ross: would you be interested in contributing to OpenSim in order to realize this? I'd like us to come together to discuss and then implement the features that will satisfy everyone's needs.

Sina:
The state derivatives for the multibody system (the generalized accelerations) must all be computed at the same time, but you can still ask for them one-by-one by going to individual coordinates. For the muscle state derivatives, you could also go to the muscles one-by-one and ask them for their derivatives. However, I think the index sets that Betts mentions in his textbook is a better way to deal with calculating sparse finite differences (by formulating the perturbation in a clever way). Even better yet, I think OpenSim should have a routines to calculate numerical sparse Jacobians for you. I think this is more desirable than a single method that allows you to get individual state derivatives, since such a method still might require calculating most of the state derivatives at once anyway.

[Ton van den Bogert]

It's great to have your support for this, Chris.

I personally won't be able to do coding and testing, but I am more than happy to give advice, since I have done this before. My student Brad Humphreys (who has been following this thread) may be interested.

Ideally, Opensim would provide analytical derivatives, yes. That is going to be a lot of work. Muscle state equations are not that hard to differentiate. Muscle force will have derivatives with respect to joint angles and fiber length (via series elasticity), and these derivatives would propagate into joint torques and finally the generalized accelerations. Differentiating the multibody dynamics would require input from Sherm and I am not sure if that is even feasible to add to Simbody. This differentiation is a lot easier when multibody dynamics is written in implicit form (essentially the inverse dynamic form, with generalized force being explicit). But it will still be an enormous undertaking, from what I understand of Simbody.

For the time being we can probably live with finite differences and see what we can do to make that work as well as possible.

For optimizers, it is important that the dynamics equations (which become constraints in the collocation scheme) are twice differentiable so IPOPT can estimate the Hessian of the Lagrangian. This would require some attention in the muscle models. Hard saturations etc. must be removed. For instance, I currently allow my force-velocity curve to go negative when shortening faster than Vmax, to avoid the sudden saturation at zero. The series elasticity will not transmit negative forces, so it seems to work, but may need some serious review before putting that in an Opensim muscle model.

Sina: like Ajay, I don't think you need to do muscle dynamics one at a time. The muscles are not coupled to each other (only to the skeleton), so you can apply the same finite difference to all of them at the same time, ask Opensim for the perturbed state derivatives, and then compute the entire Jacobian matrix for muscle dynamics (which is diagonal). This saves a lot of time, because there are so many muscles. It's one Opensim function call rather than many. Jacobian of the multibody dynamics must be evaluated as a dense block (as Ajay said), but in the implicit formulation, there is no coupling between different branches of the multibody tree. Some time saving is possible. The Jacobian df/dx has a checkerboard pattern in the implicit multibody dynamics. You can apply a finite difference to four angles at the same time, for instance two hip angles and two shoulder angles, and get four columns of the Jacobian with one function call. This is another advantage of the implicit formulation, in addition to better numerical conditioning when the mass matrix is near-singular.

The methods that Ross refers to would automatically detect those opportunities for quick Jacobian evaluation. In our early work we used SNOPT which has that capability built in. I am not sure it is the job of Opensim to provide that information, it is more an optimizer issue.

Further time savings are possible by doing the finite differences on multiple processor cores. So there are plenty of easy opportunities for speeding up the Jacobian calculations without doing analytical derivatives in Opensim.

Implicit multibody dynamics may be relatively easy to add to Opensim, and that would also be an opportunity to provide an implicit muscle model that is twice differentiable, so we have full implicit system dynamics available.

I am curious how much time Brian’s optimizations spend in Opensim vs. IPOPT (it may be in his paper...). This may help decide what the priorities should be for Opensim. If IPOPT takes most of the time, finite differences are fine, but we may want to make sure that the problem is numerically well conditioned and differentiable to make it easier to solve. If Opensim takes most of the time, the finite differences are still a bottleneck.

Ton