Energy drift

[chudo19]

Good day. Please pardon for my bad english.
As part of my reasearch i have implemented particles with constraints physics as described by baraff&witkin docs.
But due to numerical integration error all constraints drift from proper positions. Doc says to use spring dumper feadback (Baugmante stabilization) to solve this problem. But i found that using such stabilization introduses huge energy drift. Taking simple example - "bead on a wire" it takes 5-7 second for radial velocity to go to infinity.
Doc says nothing about energy drift so i tryied different integration schemes. And i discovered that any of Euler ,MidPoint or RungeKutta6
introduse energy drift even without feedback step.
Only verlet integration , (Xn+1 = 2*Xn - Xn-1 + at^2 ) does not seem to introduce new energy to the system. At least at "bead on the wire " case.
But as i read through internet i see people often use simple Euler integration and no body admits any energy drift problems.
Can anybody explain, please, how to solve such issues?
Is there anybody with same problems as mine?

[chudo19]

i still have problem.

can anybody help ?

[mewert]

Energy drift is probably a problem with most solvers out there and it only really occurs due to angular error. It's usually safe to ignore. It's common to introduce artificial global damping as well, which is fine, kind of like air resistance. You don't usually have a situation where the system is free to move in such a way that the energy drift accumulates and becomes noticeable. Your rag-dolls arms will be intercepted by their head or a wall before it does a full 360 around it's shoulder joint.

If you calculate your constraint frame as if it was the end of your time-step, rather than the start, you may be able to get a pseudo implicit kind of effect, which may damp out your energy.

I've seen claims that a projected Gauss-seidel solver with an euler integrator is actually an implicit time-stepping scheme. Maybe someone can explain why energy drift is still a problem with such a solver?

[chudo19]

Energy drift is probably a problem with most solvers out there and it only really occurs due to angular error

i have tested Novodex engine with pendulum composed of 5-7 bodies and didnt notice any positive or negative drift after 3-5 minutes of simulation.
How can it be ? Altough i know they use euler as base integrator.

If you calculate your constraint frame as if it was the end of your time-step, rather than the start, you may be able to get a pseudo implicit kind of effect, which may damp out your energy.

Can you explain more on the topic. Please.

[ngbinh]

If you have contacts/collision between bodies, it's very hard to reserve energy. If they are all free, then the energy should be the same unless you have something wrong in your integrator.

[pspcontest]

If the drift is due to penetration depth recovery, there are solutions (besides using continuous collision detection of course):
Split the constraint solver impulse, and only use the penetration depth impulse for the integration, but not for the next frames velocity.

See Erin's comment on the ODE mailing list:

I didn't have time to check the patch, but I updated my 2D physics engine
Box2D. This engine uses sequential impulses, which are mathematically
equivalent to PGS. Here's how I did it:

- Each body has a bias velocity that is only used for correcting overlap.
- A separate bias impulse is accumulated for each contact.
- The bias velocity and impulse are initialized to zero each time step (no
warm starting).

The bias impulse and velocity are completely independent of the regular
impulses and velocity, so they could be solved in parallel. However, it was
convenient to solve them together and that should be more cache friendly.

You can get the code here: http://www.gphysics.com/files/Box2D.zip

You can switch between the two methods by using the macro
BIAS_PRESERVES_MOMENTUM in Arbiter.cpp. With the split velocity I was able
to use a large bias factor (0.. This means that overlap is reduced 80% per
step. Such a large bias factor leads to serious bouciness with the merged
velocity.

I also tried using a split velocity for joints. The results were not nearly
so nice, especially for my suspension bridge test. I think that particular
arrangement is more stable with some compliance in the joints.

Erin

[mewert]

You say your implementation is based on the Barraff/Witkin docs. So I assume you are enforcing your constraints as an application of a force? i.e. dv/dt ( acceleration ). Novodex, etc. skip acceleration and apply impulses, rather than forces i.e. dx/dt

If you split your impulses, as was mentioned, it appears that this is, in fact, an implicit time-stepping scheme. The penetration recovery impulse is an implicit update for the previous time-step. It's just being applied a frame later. An implicit update looks at the state of the system at the end of the time-step when calculating a solution. So to understand how the penetration recovery works as an implicit update, you go back a step which was looking forwards to the current step... confused? So the result should be a loss of energy in you system. Try scribbling diagrams of what happens to your bead on a wire example after over a couple of time-steps to see what I mean. This is in 2D, though maybe it doesn't apply in 3D so well, I'm pretty sure I recall my rag-doll shoulder constraint spiraling out of control with large-time steps. But that was awhile ago, maybe I'm thinking of what happened before I fixed a bug...

If you are applying forces ( rather than impulses ) and you take the constraint directions from what they would be at the end of the time-step, that would be like converting it to an implicit update ( but not as strong for a complicated system, so it might not work ).

Anyways, switch to applying impulses rather than forces, ( the baumgart thing seems to also be an implicit update ), if that doesn't work, do the impulse splitting technique.
If you are still gaining energy, check your constraint calculations, something is probably off.

[gee]

If the drift is due to penetration depth recovery, there are solutions (besides using continuous collision detection of course):
Split the constraint solver impulse, and only use the penetration depth impulse for the integration, but not for the next frames velocity.

See Erin's comment on the ODE mailing list:

I didn't have time to check the patch, but I updated my 2D physics engine
Box2D. This engine uses sequential impulses, which are mathematically
equivalent to PGS. Here's how I did it:

- Each body has a bias velocity that is only used for correcting overlap.
- A separate bias impulse is accumulated for each contact.
- The bias velocity and impulse are initialized to zero each time step (no
warm starting).

The bias impulse and velocity are completely independent of the regular
impulses and velocity, so they could be solved in parallel. However, it was
convenient to solve them together and that should be more cache friendly.

You can get the code here: http://www.gphysics.com/files/Box2D.zip

You can switch between the two methods by using the macro
BIAS_PRESERVES_MOMENTUM in Arbiter.cpp. With the split velocity I was able
to use a large bias factor (0.. This means that overlap is reduced 80% per
step. Such a large bias factor leads to serious bouciness with the merged
velocity.

I also tried using a split velocity for joints. The results were not nearly
so nice, especially for my suspension bridge test. I think that particular
arrangement is more stable with some compliance in the joints.

Erin

Hello,

I'm sorry but I cannot see the difference between doing that and doing another LCP at the position update step. Why doing it at the velocity one ?

My main problem today is that error correction applied on the velocity completely brake the contact caching thing (and of course it's not stable at all with stacking). I have other problems but which can be solved with all I can read here and there.

By the way I'm writing a simulator based mostly on kenny's thesis, and I would like to thank him, and all of you for every thing I've learned through your papers and this forum.

[Erin Catto]

Yes, the impulse splitting can be done in separate loops, but it may be more efficient to combine the loops. Only velocity impulses should be cached.

[gee]

Hello,

thanks for the reply but I cannot see how you do that. I've tried to look at box2D, but I'm not getting it completely.

If I do a 'normal' MLCP to resolve collisions I will only get one impulse array for my contacts right ? then if I want to have impulses with and without the baumgart error term I need to do this twice no ? I think this is what I'm thinking wrong, but as I'm not the doing doing the LCP I might miss some knowledge here, and my colleague does not seem to know another way as well.

By the way, on monday I tried a position correction only, and stacking was not that better also, and to do a good one I had to recalculate the constraints which is doing twice the job per iteration, for a not so good result (this without contact caching).

Thanks

[Dirk Gregorius]

You can think of it like this. Instead of solving:

J*W*JT * lambda = -C(x) / dt - J*( v(t) + W * f_ext * dt )

you solve the two MLCPs individually:

J*W*JT * lambda1 = -C(x) / dt
J*W*JT * lambda2 = -J*( v(t) + W * f_ext * dt )

Instead of superpositioning both results, you bybass the velocity for the position correction (Baumgarte term).

I think this is little different from the original post-stabilization idea like e.g. in the Cline thesis, where you really compute a second MLCP *after* the position update - with a new Jacobian and delta velocity. The idea of Erin is that the position correction should not effect the velocity state (the persistent momentum) of the rigid bodies. Note that this works actually not well for all scenarios. I think Erin mentioned in this forum that he ran into problems with the brige example. For contacts it seems to be a nice idea and let's contacts appear much more rigid...

-Dirk

[gee]

Hello,

I'll try that later and see how this works. But this is still solving 2 LCPs as I thought, not just one

Also another way would be to compute
J*W*JT * lambda = -C(x) / dt - J*( v(t) + W * f_ext * dt )

and then compute :
J*W*JT * lambda = -J*( v(t) + W * f_ext * dt )
using the old lambda values for warm start no ? As the errror should be small the lambda values should be close to new ones I think.

I will try this for contact, and maybe keeping it as the first equation for joints (well I will try both, but as you said it does not seem to be always a good idea).

Thank you for your help.

[Dirk Gregorius]

Basically you are right, you can see it as two LCPs if you want. The advantage of what I suggested is:

1) The Jacobian is the same for both LCP, so when implemented correctly the second LCP comes for free
2) Since you now know the "penalty" impulse to catch up the position error you can simple find the displacement from this without adding to the velocity.

dv = P * m => dx/dt = P * m => dx = P * m * dt

So you apply this displacement directly and dont add it to the velocity

HTH,
-Dirk

PS: Please note that in your formula here you are missing the scaling constant for the Baumgarte term - this:

J*W*JT * lambda = -C(x) / dt - J*( v(t) + W * f_ext * dt )

must be for example

J*W*JT * lambda = -0.1 * C(x) / dt - J*( v(t) + W * f_ext * dt )

See Erin's GDC 2005 presentation for a derivation on how to choose values here...

[gee]

Hello,

1) True.
2) Yes, that I know, I just thought from Erin answer that there was a way with only 1 LCP (trying to resolve 2 systems at the same time did not seem like a good idea so I wondered).
3) And yes I know, I was just using something simple. By the way using a baumgarte correction for now, I have to use a really small coefficient, while with joints 0.8 is fine most of the time. With contacts I have something around 1e-6 (well that is for around 200 balls in contact, for 3-6 balls a higher value is ok).

Thanks

[Dirk Gregorius]

I made the same observation (e.g. when implementing cloth using sequentiel impulses). There seems to be a dependency between the number of constraints and the Baumgarte coefficent. The formula I used was:

tau ~ 1 / Sqrt( n )

where n is the number of constraints

Actually I read through a lot of papers recently in order to find better stabilization. The problem with Baumgarte stabilization is that it combats the drift, but not necessarily on the best trajectory. I think for contacts it is really ok, but for several joints building an articulated structure it might be a problem (see e.g. papers of Uri Asher on DAE). I am currently interessted in high quality joints (e.g. to couple animation and physics), but iterative solvers don't seem to be good enough. In the newest paper of R. Weinstein she actually moves from here iterative solver approach to a direct least square solver that minimes the residum...

-Dirk