Need ODE solver for C++ Drake

[psiorx]

Creating this issue to spawn a discussion and collect some notes on potential open source options for an ODE solver that could work in Drake.

The solver will be at the heart of most of the simulations so giving some thought as to whether we should build this from scratch or integrate an existing package is in order.

To start off, I found this:

http://headmyshoulder.github.io/odeint-v2/

In particular, it has a good implementation of an ode45 equivalent. Sounds like it will be included in Boost soon and I know we're trying to avoid dependencies on Boost but it's header-only and supposedly highly performant. Any thoughts?

[psiorx]

@avalenzu just mentioned this one out of LLNL:

https://computation.llnl.gov/casc/sundials/main.html

[psiorx]

Looks like odeint has already been integrated with Boost(header-only) and ships with the latest version. However, the version shipped with Boost doesn't seem to have built-in support for Eigen types.

One option is to do the computation on a std::vector or boost::array and wrap that in an Eigen map. It also appears that the latest version on GitHub does support Eigen types but it's unclear whether this will make it into the Boost distro.

[psiorx]

Update:

Got odeint working with Eigen types of different flavors (including AutoDiffScalar!). This could have potentially interesting implications (e.g. Shooting methods) if the gradients are maintained automatically throughout a simulation.

Using Eigen types as the underlying storage medium also seems to give a performance boost compared to std::vector especially for fixed size vectors.

Either way, the library is very flexible and fast and I'd like to move towards pulling this into Drake. Have we settled on an interface/contract for simulation within C++? In matlab, we're returning xtraj's but in the recent pendulum runLQR example, simulate returns nothing and simply calls output repeatedly.

[RussTedrake]

Excellent. I think we should move the trajectory classes over (note that PPTrajectory already exists in some form and could be reused.. not sure if it needs to be cleaned up) and continue to have simulation return trajectories.

[RussTedrake]

I commented on the pull request, but I now think odeint doesn't do everything we need. My two cents is that the main things that we need are:

• runge-kutta 4
• runge-kutta 4/5 + event detection (the variable step solvers will be almost useless without that)
• the trivial "integrator" for discrete systems
• support in rk4 and rk45 for mixed discrete/continuous systems

odeint doesn't do the event detection; i don't think. not sure about the mixed-discrete/continuous updates. neither is very hard. i'm thinking right now that it probably makes sense to just roll our own.

[RussTedrake]

Art kuo recommended sundials (which definitely does zero-crossings etc), looking at open modelica, and at guggenheim's code for taking bigger steps by using jacobians of the dynamics. He also is cleary under the impression that we could just use numerical recipes.

[sherm1]

I don't think you would be happy with the BDF multistep integrators in Sundials. Although they can work beautifully under the right circumstances, I find them too touchy, especially when used with discrete events (requiring restarts at low order). We have a custom one in Simbody written specifically for us by CVODE's author Radu Serban. When it works it's great but it has too many weaknesses for my tastes. A one-step RK method will be much more robust.

I do agree that a roll-your-own is likely to work best here if you want error control, hybrid discrete/continuous integration, and automatic differentiation. Simbody has an RK Merson that I like better than RK45, but either would work and not take a lot of code. Good event isolation is tricky though!

[tkoolen]

I'm not confident that using IntegratorType::RungeKutta2 is going to be sufficient for a decent Valkyrie simulation, so this is still a blocking issue for that.

[sherm1]

We're going to add a higher order variable step method to Simulator next; that should be sufficient for Valkyrie.

[sherm1]

I'll assume the variable-step RK3 we have had in Drake for ages addresses this issue. Other requests should get their own issues.