Impulse-based dynamic simulation system

[Jan Bender]

At the moment I am writing my PhD thesis about impulse-based dynamic simulation. In contrast to Mirtich's work joint constraints are simulated by impulses as well as collisions and contacts. So if you are interested in my work, take a look at:

http://www.ibds.de.vu

[Erwin Coumans]

Thanks for the link, that looks interesting.

We present a new method for the dynamic simulation of linked rigid body systems

I haven't read the details, but it sounds very similar to what Bullet Physics Library is doing. How does this compare to Bullet? Bullet uses sequential impulses too, for contacts, collisions (no distinguishment made) and constraints. I would call most of such iterative approaches just relaxation, introduced in the rigid body simulation world by Bart Barenbrug and Kees van Overveld. Sequential impulse, as described by Erin Catto is another implementation of relaxation.

Could you tell us what is new about your approach, compared to say Bullet?

Thanks,
Erwin

[Jan Bender]

My approach works for joint and contact constraints as follows. A simulation step consists of these parts:

1) Determine the position of the joint or contact points after the next simulation step (in fact a preview of the future is computed).

2) Make an approximation of the velocity change that is needed to satisfy the constraints. Due to this approximation I get a linear equation for the impulse that can be easily solved. Here I use the same matrix K as in Bullet. Then the impulse is applied.

(R. Weinstein used a nonlinear equation to get the impulse in her paper. The rest of "her" algorithm is the same as the one I already published in the year 2003.)

3) Since I used an approximation and the constraints may depend on each other I continue iteratively with step 2) until all constraints are satisfied.

4) Now I can proceed in time and the constraints wil be satisfied after the time step due to the applied impulses.

5) After the time step I can determine impulses to satisfy all velocity constraints in a similar way.

This algorithm converges to the correct solution. The proof can be found in one of the papers on the homepage.

If you want to know more about this method, please feel free and ask me.

Jan

[Erwin Coumans]

My approach is very similar, so the difference must be in the details.

My approach works for joint and contact constraints as follows. A simulation step consists of these parts:

1) Determine the position of the joint or contact points after the next simulation step (in fact a preview of the future is computed).

2) Make an approximation of the velocity change that is needed to satisfy the constraints. Due to this approximation I get a linear equation for the impulse that can be easily solved. Here I use the same matrix K as in Bullet. Then the impulse is applied.

When do you perform collision detection? In between step 1 and 2?

My approach in Bullet also predicts an unconstraint preview of the future, for the object transforms. Then continuous collision detection can be used to calculate the time of impact, given a constant linear and angular velocity, which returns a fraction between [0..1] between current transform and this predicted future / next simulation step transform.

Also, I integrate/predict the unconstraint velocities at the beginning of the timestep.
At the moment, I use the current frames transform for the contact points and joints position, and use the unconstraint 'predicted' velocity of the next timestep.However, all information is available to use the next simulation steps contact point and joint positions. It would be intesting to compare both methods.

The interleaving and order of individual constraint corrections and measurements leave a lot of freedom and choice.

Have you considered how the order of solving constraints impacts the results? I found that an ordering of back and forward over the list is better then just going forward. Also I get better results by first solving all positional and velocity constraints, and then the friction constraints. This provides a good estimate for the total friction magniture for the current frame.

How do you determine the friction magnitude?
Do you use last frames friction normal impulse/force?

(R. Weinstein used a nonlinear equation to get the impulse in her paper. The rest of "her" algorithm is the same as the one I already published in the year 2003.)

3) Since I used an approximation and the constraints may depend on each other I continue iteratively with step 2) until all constraints are satisfied.

Step 3 is called relaxation.

4) Now I can proceed in time and the constraints wil be satisfied after the time step due to the applied impulses.

5) After the time step I can determine impulses to satisfy all velocity constraints in a similar way.

Do you consider continuous collision detection?
It can improve the 'prediction' of future positions in step 1. If you don't want the added complexity of subdividing the timestep, you can consider only using continuous collision detection between dynamic and static objects (ignoring the dynamic versus dynamic case).

Thanks,
Erwin

[Dirk Gregorius]

I found that an ordering of back and forward over the list is better then just going forward.

I tried this as well, but my experience was that it is not worth the effort. Can you give examples where this improved the convergence? Was this better for at least the majority of situations and can we see the improvement in the examples?

[Jan Bender]

When do you perform collision detection? In between step 1 and 2?

Collision detection and resolution I perform before step 1. I made the assumption that collisions occur in an infinitesimal time intervall. So collisions can be resolved before the time step. Just contacts and joints have constraints that lasts for the whole time step.

It would be intesting to compare both methods.

Where can I find some detailed information about your method? Do you have a paper about it?

Have you considered how the order of solving constraints impacts the results?

I tried a few things but I didn't find an order that worked for all models. Which order do you take if you have a closed loop in the model?

How do you determine the friction magnitude?

I use the Coulomb friction law. So dynamic friction is computed depending on the normal impulse of the actual iteration. Static friction occurs, if the sum of all friction impulses and the sum of all normal impulses satisfy the condition defined by the law.

Do you use last frames friction normal impulse/force?

No, not yet. Dirk told me already about this and I really want to try this as soon as possible.

Do you consider continuous collision detection?

No, I use a discrete collision detection. In my thesis I mainly focus on joint constraints, collision response and contact handling. I didn't want to spend too much time on collision detection. So I haven't implemented an own one.

Jan

[Erwin Coumans]

When do you perform collision detection? In between step 1 and 2?

Collision detection and resolution I perform before step 1. I made the assumption that collisions occur in an infinitesimal time intervall. So collisions can be resolved before the time step. Just contacts and joints have constraints that lasts for the whole time step.

So you make a difference between colliding contact and contact points?
In Bullet I don't make any distinguisment. As soon as objects 'collide', the collision detection just adds contact points. Those are dealt with in the solver.

It would be intesting to compare both methods.

Where can I find some detailed information about your method? Do you have a paper about it?

Bullet's dynamics method is not novel, so no paper.
It is most recently described in Erin Catto's and Kenny Erlebens PhD (Projected Gauss Seidel), see below for links. It is really just relaxation using impulses, for contacts and other constraints.

Parts of the continuous collision detection are novel, and published in my draft paper here:
http://www.continuousphysics.com/Bullet ... ection.pdf
But you didn't express interest in CCD, so you can skip that

Have you considered how the order of solving constraints impacts the results?

I tried a few things but I didn't find an order that worked for all models. Which order do you take if you have a closed loop in the model?

You can solve in order first to last, and next iteration last until first. In practice I found it handles stacking a bit better. Also, I use warm starting, which means starting with previous-frames impulse as first estimate. This makes convergence a little better too.

How do you determine the friction magnitude?

I use the Coulomb friction law. So dynamic friction is computed depending on the normal impulse of the actual iteration.

During the actual iteration, the normal impulse is not really a valid approximation, especially when you perform clipping on the accumulated normal impulse. See Erin Catto's paper/slides on this clipping of the 'accumulated' impulse.

See the slides here
http://www.gphysics.com/files/GDC2006_ErinCatto.zip
This one is a bit older, still interesting.
http://www.gphysics.com/files/IterativeDynamics.pdf
Bullet implements this in 3D.

As we discussed on this forum, iteratively applying sequential impulses with clipping the accumulated impulse is mathematically equivalent to Projected Gauss Seidel. PGS is pretty well documented, I can recommend Kenny Erlebens PhD thesis or recent book: http://www.diku.dk/~kenny


I haven't had time to look deeper into your papers yet, but so far, the main diffference is that you are solving the positional contact constraint (penetration depth recovery?) in the future 'predicted' position.

What are the main novel ideas in your approach?

Erwin

[Erin Catto]

I arrived at my sequential impulse (SI) algorithm by studying the PGS algorithm in great detail. A coworker of mine (Gary Snethen) had the idea that the PGS algorithm was really just solving pairwise interactions and iterating over all interacting pairs. I extended this to the accumulated impulse concept and assembled the final algorithm. In its basic form, SI is mathematically equivalent to PGS.

The appeal of the SI algorithm is that it has the same performance and robustness of the PGS algorithm, but it is far less abstract. There are other benefits as well. It lets you solve in blocks, you can do n-way interactions, and you can improve the friction model.

The key to SI is coherence, as manifested in the accumulated impulse. This is essentially warm starting. You can also view it as amoritizing the solver cost over many time steps. PGS and SI require many iterations to converge. At 60Hz you effectively get a minimum of 60 iterations per second.

[Jan Bender]

[Erwin Coumans]

So you make a difference between colliding contact and contact points?
In Bullet I don't make any distinguisment. As soon as objects 'collide', the collision detection just adds contact points. Those are dealt with in the solver.

Yes. You use the same algorithm for solving collisions and contacts? Do you use any prediction for the position of the contact points? With a prediction I get very nice results since static friction can stop bodies completely. With other algorithms I made the experience that bodies on an inclined plane never stop moving totally.

It is most recently described in Erin Catto's and Kenny Erlebens PhD (Projected Gauss Seidel), see below for links. It is really just relaxation using impulses, for contacts and other constraints.

Erin Catto's and Kenny Erleben's works are already on my Must-read-this-stack.

But you didn't express interest in CCD, so you can skip that

Interest yes. Time no. I just visited a presentation of Stephane Redon about this topic at the Eurographics conference. It is very interesting but my main goal is to finish with writing my thesis. (I hope the first version will be finished in two weeks.)

I haven't had time to look deeper into your papers yet, but so far, the main diffference is that you are solving the positional contact constraint (penetration depth recovery?) in the future 'predicted' position.

In fact I use the predicted positions to resolve the contact constraints exactly and for the simulation of static friction.

What are the main novel ideas in your approach?

First of all the main idea is to solve all kind of constraints by impulses. Okay, now it is not really new anymore because many people do this. But when our research with this started nobody used impulses to make an exact simulation of joint constraints (I don't know any work about this). At my institute we make research on this method since the year 2000. In fact we started with mass point systems and continued to extend the algorithm for rigid bodies. The algorithm has the main advantages that it can handle closed loops without special treatment and no kind of extra stabilization like Baumgarte is needed. The newest research was done on collision and contact resolution which I have presented a few month ago on the CASA conference. Now I have again a paper in review. In this paper I describe a fast method on how all impulses of a complex system can be computed at once. Also we did some research on methods of higher order. Alltogether we wanted to have an algorithm that allows an exact simulation and that is fast and easy to implement. Many of the physics engines for games have the goal to deliver just plausible results and to be very fast. Since we can define how exact the simulation should be, we can also make plausible simulations which are very fast.

I think in the next time I really must read the thesis of Erleben and the papers of Erin Catto. This could be very helpful to understand the main differences in the approaches


Jan

[Jan Bender]

[Erwin Coumans]

During the actual iteration, the normal impulse is not really a valid approximation, especially when you perform clipping on the accumulated normal impulse. See Erin Catto's paper/slides on this clipping of the 'accumulated' impulse.

Could you tell me a reason for this? In Erin Catto's work I just found something about his simplified friction model.

Jan

[Dirk Gregorius]

During the actual iteration, the normal impulse is not really a valid approximation, especially when you perform clipping on the accumulated normal impulse. See Erin Catto's paper/slides on this clipping of the 'accumulated' impulse.

Could you be a little bit more precise, this is interessting . If I understand you correctly you mean that the difference impulse (as opposed to the accumulated impulse) is not a valid approximation, right? But if I understand Erin's work correctly the friction model is based on the accumulated impulse which is (hopefully) converging against the correct solution and therefore is indeed a good approximation. You might be correct regarding the first impulse model you used in Bullet. Here the normal impulses were of course not a good approximation. Did I miss something here???

Regarding constraint ordering:
The GS scheme needs that the system matrix is diagonal dominant in order to converge. The means the diagonal element in a row must be greater than the sum of all other elements in this row. An example

100 1 1
1 100 1
1 1 100

This is a nice diagonal dominant matrix and GS will find a nice solution for this. If you now swap the first and second row

1 100 1
100 1 1
1 1 100

This matrix is not diagonal dominant and the GS will not find a solution. That means that you can have equivalent linear systems, but one time you find a solution and the other time not. I think this is what they try to catch in the ODE when they arbitrarily mix the constraints in each iteration. Again this is Voodo for me. On the other hand I would like to know if there is a constraint order (e.g. joint, non-penetration, friction) that has the best changes of being diagonal dominant. Ideas?

-Dirk

[Erwin Coumans]

[Erwin Coumans]

During the actual iteration, the normal impulse is not really a valid approximation, especially when you perform static friction clipping on the incremental normal impulse. See Erin Catto's paper/slides on this clipping of the 'accumulated' impulse.

Could you tell me a reason for this? In Erin Catto's work I just found something about his simplified friction model.

Jan

I'll try to clarify this statement below.

If I understand you correctly you mean that the difference impulse (as opposed to the accumulated impulse) is not a valid approximation, right?

Indeed, I amended the line into "especially when you perform static friction clipping on the incremental normal impulse".

Basically, during the iterations, the estimated incremental normal impulse can go fairly wild, and it gets corrected iteratively. The correction impulses mean that sometimes you get negative normal impulses.

I'll try to give a simple example where things might fail, unless you perform friction at the end of all iterations (where the normal impulse is properly known):

Say we use 2 iterations in the solver, and we have 2 contacts shared between 2 rigidbodies, one static and one dynamic (so a cube resting on another cube on its edge for example).

Now, the friction model might clip at a certain value MaxStaticFriction, say 1.0

first iteration, we traverse the contact points, and calculate the normal impulse, which happens to be 1.5 and 0.5 respectively for each contact point.
Next iteration, the incremental normal impulse is -0.5 and 0.5, so the accumulated normal impulse is 1.0 and 1.0 for both contacts.

If you perform friction at the end of the 2 iterations, no clipping is done, and all is fine. If you would interleave the friction with each iteration, you would clip the friction for the first contact from 1.5 to 1.0. Then next iteration you perform either -0.5 or 0.0 (depend wether you allow negative normal impulses to apply friction). This is obviously a wrong and different result.

Therefore, you can either use the last frames normal impulse magniture, or perform the friction iterations after all 'normal impulse' iterations are done. I do the latter in Bullet. There is also an issue with this, but not a big issue for me: it means that the interaction between friction and contact is not perfect. Such interaction happens only between the frames, not between the iterations in my approach. Problem with using the previous frames' normal magnitude is that you need to approximate it for new contacts (where there is no 'previous' frame). Doing such approximation can become voodoo art

[Erwin Coumans]

Regarding constraint ordering:
The GS scheme needs that the system matrix is diagonal dominant in order to converge. ....

Take a stack of spheres. If you traverse the contact constraints from top to bottom you get worse convergence then bottom to top, because the effect from the ground (infinitely heavy object) is not propagated at all during the first iteration. If you go bottom up, the 'ground' effect is nicely propagated.

If you randomize the constraint ordering, like ODE does, this effect is better propagated, but you loose simulation reproducability (unless you have a reproducable random number generator), and debugging becomes a nightmare.

THerefore, just simply traversing first -> last constraint and then last->first constraint can improve in practice. So it is a middleway between not re-ordering at all, and randomizing. I think Kenny Erleben does analysis to do the ordering, but that fails in closed-loop situations I think. My simplified forward/backwards interleaving doesn't fail in closed loops.

[Dirk Gregorius]

Hmmm, Erin actually says it is the other way around

Here's one thing that can improve convergence. I'm not sure if it has been mentioned before. Consider a vertical stack of boxes. Your intuition may tell you that solving the PGS according to the order of the stack may help. This is true. It turns out that solving from the top downward to ground effectively gets you a free iteration versus the opposite order.

http://continuousphysics.com/Bullet/php ... ration#374

Regarding friction:
I apply the friction at the manifold center as we discussed a while ago . And of course I have all normal impulses then. What I don't see is the problem with wild impulses. I compulte the friction using the accumulated impulse, not the intermediate "Delta Impulses". I don't see how this should get wild, since the accumulated converges (hopefully) against the correct solution each iteration as the friction does then as well. Right? I tried the last frame stuff and it was much worse then what I do now. I must admit that from the very beginning I always computed all normal impulses and then the friction constraints. Maybe that saved me some issues and I don't see the problem here

We need more discussions of this kind on your forum

[Erin Catto]

I actually did the sphere stacking test awhile ago with PGS. Top down gives you approximately a one iteration advantage over bottom up. I also found this counter intuitive. You might think of it as getting the full weight to ground quicker.

With regard to going back-and-forth, this is called Symmetric Successive Over Relaxation (SSOR). I also thought this sounded good. Apparently it has provably worse convergence than regular SOR (of which GS is a subset). See Matrix Computations by Golub.

Randomization of constraints kills the cache. I saw solver performance double on the PS2 by turning off randomization.

Dirk, convergence of GS is determined by the eigenvalues. These should be immutable under row and column swaps. I think ordering can only affect convergence in small ways (like 1 iteration).