Energy conservation/extra sources of damping

mentisn0b commented 7 years ago

We are doing some physics simulations using DART. Even with joint damping set of zero, there appears to be another source of damping causing loss of energy over time.

Situation is: single sphere attached to world with ball joint. Set spring stiffness 0..2 to 1, and give the dofs 0..2 a nonzero initial position. The system loses energy quickly over time and oscillations decay rapidly.

Known sources of damping are:

dof->setDampingCoefficient()
ERP

The first was set at zero, the second did not appear to have an impact.

Are there other (hidden) sources of damping?

mxgrey commented 7 years ago

Numerical integration can do funny things to system stability. Explicit Euler method tends to make systems appear less stable (less damped) while implicit Euler method tends to make systems appear more stable (more damped), although there are exceptions to each of those statements.

Mathematically speaking, it's impossible to have perfect conservation while using numerical integration unless the order of the equations of motion is less than or equal to the order of the integration scheme (which will not be the case for an articulated dynamic system). DART uses numerical integration during simulation, although I'm not sure how best to characterize the scheme... it resembles something like explicit Euler or maybe leapfrog integration.

mxgrey commented 7 years ago

To at least reduce the issue that you're experiencing, you might want to try reducing the timestep that is used during simulation. If that doesn't help and the system still rapidly loses energy, then further investigation might be warranted.

mentisn0b commented 7 years ago

Changing the timestep does not appear to make a difference. As a minimal reproducing testcase, take the rigidChain example from the master branch, and replace the body of main() with:

[...]
// create and initialize the world
Eigen::Vector3d gravity(0.0, 0.0, 0.0);
myWorld->setGravity(gravity);
myWorld->setTimeStep(1./2000);

//myWorld->getConstraintSolver()->setLCPSolver(std::unique_ptr<dart::constraint::LCPSolver>(new dart::constraint::DantzigLCPSolver(myWorld->getTimeStep())));
myWorld->getConstraintSolver()->setLCPSolver(std::unique_ptr<dart::constraint::LCPSolver>(new dart::constraint::PGSLCPSolver(myWorld->getTimeStep())));

int dof =  myWorld->getSkeleton(0)->getNumDofs();
std::cout << "Skeleton has " << dof << " dofs\n";
Eigen::VectorXd initPose(dof);
for (int i = 0; i < dof; i++) {
myWorld->getSkeleton(0)->getDof(i)->setDampingCoefficient(0);
initPose[i] = dart::math::random(-0.5, 0.5);
myWorld->getSkeleton(0)->getDof(i)->setSpringStiffness(1.);
}
myWorld->getSkeleton(0)->setPositions(initPose);
[...]

Replace MyWindow::timeStepping() body so it only calls mWorld->step().

chain.skel is:

<?xml version="1.0" ?>
<skel version="1.0">
    <world name="world 1">
        <skeleton name="skeleton 1">
            <body name="dummy">
                <gravity>0</gravity>
                <transformation>0 0 0 0 0 0</transformation>
                <inertia>
                    <mass>.1</mass>
                    <offset>0 0 0</offset>
                </inertia>
                <collision_shape>
                    <geometry>
                        <sphere>
                            <radius>0.01</radius>
                        </sphere>
                    </geometry>
                </collision_shape>
                <visualization_shape>
                    <geometry>
                        <sphere>
                            <radius>0.01</radius>
                        </sphere>
                    </geometry>
                </visualization_shape>
            </body>
            <joint type="ball" name="joint 0">
                <parent>world</parent>
                <child>dummy</child>
                <transformation>0 0 0 0 0 0</transformation>
            </joint>
        </skeleton>
    </world>
</skel>

Results are the same for dt=1/2000 and 1/20000, fast damping to rest condition.

jslee02 commented 7 years ago

For clarification, did you use the original MyWindow class without modification? If so, this would be because of the adding damping here.

mentisn0b commented 7 years ago

Hi jslee02, we modified the original MyWindow class to turn off the damping introduced via setForce(). The body of the MyWindow::timeStepping() function contained only a call to mWorld->step(), and setForce() was never called.

mxgrey commented 7 years ago

I've tested the situation exactly as described by @mentisn0b. Here are my observations:

Starting from a randomized configuration with 1.0/2000 time step, 1.0 spring stiffness, no damping, and no gravity, the rigid chain will initially exhibit some skittish behavior where the links bounce around wildly. Within roughly 10 seconds of simulation time, this skittish behavior converges into a somewhat calm periodic swinging. Over time, this periodic swinging continues to converge, and all the links move in sync with each other.

After more than six minutes of simulation time, the converged periodic swinging behavior seems to remain steady with no qualitative signs that energy is being lost from the system (although it may oscillate out of the camera's view plane, so you may need to rotate the camera with left-click to see the full extent of the oscillation).

This behavior seems sensible to me. Intuitively, I would fully expect the multi-body interactions to produce effects that appear to "damp" the relative motions of the bodies by syncing their motion into something periodic. I think energy is being redistributed through the system, but not lost from it. If my intuition is incorrect, then please let me know.

If you are certain that you are seeing the system come to a complete rest within some amount of simulation time that's less than ~6 minutes, then that means we are observing very different outcomes, in which case we'll need to verify that we run the same test code and then perhaps consider what platforms we're each using (or other confounding factors).

mxgrey commented 7 years ago

Testing with a time step of 1.0/200, I observed the system coming close to rest after 3 minutes of simulation time, but a time step of 1.0/2000 only had very minor energy loss throughout 12 minutes of simulation time. (To be clear, "simulation time" here refers to N*dt where N is the number of simulation steps and dt is the size of the time step). This suggests to me that any damping that is being experienced is likely due to typical numerical integration error, unless you're observing significantly different results than I am.

Edit: Also, a time step of 1.0/1000 (which is the default time step) had a noticeable decay in energy after 6-7 minutes of simulation time and started coming close to rest after 15 minutes of simulation time.

mentisn0b commented 7 years ago

mxgrey, thanks for the comprehensive response. Apparently I am getting different results. Just to be sure we're working on the same code, I did a fresh clone from the master branch.

For me, the rigidChain demo, with the damping via setForce() in the time stepping function disabled, it blows up after about 1 minute of world time (dt = 1/2000).

I print the energy at each step via mWorld->getSkeleton(0)->computeKineticEnergy() and computePotentialEnergy(). Strangely, potential energy is negative (approx. -1.8) in the initial timestep. Eventually it goes to positive infinity. The negative energy is apparently a separate issue (should an issue be filed?)

As the most simple check for non-physical damping phenomena: a sphere with diameter = 1 and mass = 0.01 bouncing back and forth around origin (translation joint, stiffness = 100) loses >3/4th of its energy after 30 sec of world time (dt = 1/2000). For dt=1/20000, the loss is only about 1/3rd. The loss apparently gets worse the smaller the mass is compared to 1.

mxgrey commented 7 years ago

Would you be able to fork the DART repo and link to your code? That should guarantee that we're looking at the same thing.

Negative potential energy is not inherently problematic in the presence of gravity, because if the "zero" point gravitational potential energy is above the CoMs of the system, then the PE can correctly end up being negative. But if there is no gravity and all the potential energy is from joint spring stiffness, then that's definitely non-physical.

mentisn0b commented 7 years ago

Please see:

Case 1: blows up. Might be acceptable behaviour because it only occurs after a long number of timesteps and the system does not have explicit damping turned on. However, seems to show different results on your machine.

https://github.com/dartsim/dart/compare/master...mentisn0b:rigidChain-without-damping?expand=1

Case 2: too much damping. The chain is reduced to 1 single segment for simplicity. Gravity is turned off and spring stiffness is set to 1. The system comes to rest very quickly, even though it should neither lose nor gain energy.

https://github.com/dartsim/dart/compare/master...mentisn0b:rigidChain-too-much-damping?expand=1

You are right about the negative potential energy; I was probably looking at a situation where the CoMs were below zero on the gravitational axis.

Thanks!

mxgrey commented 7 years ago

I’ve finally gotten around to testing the second case that you’ve linked, and I definitely see what you’re describing now, so we don’t have to worry about platform issues.

There are a few observations I’d like to make about test case 2 right off the bat: The mass-to-stiffness relationship is rather lopsided. There is a considerable amount of spring stiffness being applied to a tiny amount of mass. You can probably see how it oscillates at an extreme frequency right off the bat. These wild values result in quite a bit of numerical integration error, because the time step is not able to keep up with the extreme changes in acceleration. Increasing the body’s mass to 20 increases the time before coming to rest by a few orders of magnitude.

Furthermore, applying only gravity instead of spring stiffness shows the system to be much less damped. This again points to numerical integration error, because gravity produces a lower-order dynamic system (assuming constant acceleration due to gravity) than spring stiffness (assuming behavior according to Hooke’s Law). Higher order systems will exhibit larger amounts of numerical integration error (which will usually present itself as either excessive damping or instability) when the numerical integration scheme is lower-order than the system that it’s approximating.

All the evidence that I’m seeing points to numerical integration error. If your particular use case requires both (1) systems that have very high joint stiffness compared to their other dynamic properties, and (2) simulations that accurately capture the energy conservation of the system’s passive dynamics, then it may be necessary to alter DART’s numerical integration scheme to accommodate that. I think we’ve had discussions about making the integration scheme modular before, but I’m not aware of any ongoing work on that matter. The discussion for it would be beyond the scope of this particular issue, but if it’s something that you need, then we can begin a new issue on the topic of our numerical integration framework.

mentisn0b commented 7 years ago

Hi Michael, thanks very much for looking into this. For our purposes we do need to simulate systems with a high stiffness-to-mass ratio. As decreasing the timestep does not seem to be effective against the numerical integration error (a well-known phenomenon, I suppose, for forward Euler), it makes sense to select a different integration method. If you could point us in the right direction for how to get started on this, it would be much appreciated. I do not see how we can select a different integrator (e.g. those under dart/integrators); is it instead the case that forward Euler is more or less hard-wired into DART, because of e.g. all the integratePositions() and integrateVelocities() functions in all of the different types of joints?

mxgrey commented 7 years ago

There wouldn't be any need to make changes to integratePositions() or integrateVelocities(). I would recommend using a customized version of the World::step() command. You'll want to use a high-order Runge-Kutta method like RK4 instead of DART's native simulation step (which looks vaguely like leapfrog to me). You can either fork the DART repo and make changes to World::step() directly, or you could remove the line where you call World::step() in your application and replace it with your custom integration scheme (I would recommend the latter).

Truthfully, I'm not exactly sure how to implement Runge-Kutta when there are collisions involved. The discrete dynamics changes will likely complicate the integration procedure. If you don't need to worry about collisions for your simulations, then you'll have a much easier time implementing Runge-Kutta.

jslee02 commented 6 years ago

Closing due to lack of a response. Please re-open if you have further trouble.

dartsim / dart

Energy conservation/extra sources of damping #846