## Equivalence sequential impulses & Projected Gauss Seidel

Issues with the forum and all license/patent related discussion
Erin Catto
Posts: 324
Joined: Fri Jul 01, 2005 5:29 am
Location: Irvine
Contact:

### Equivalence sequential impulses & Projected Gauss Seidel

So ODE got the PGS optimizations from Croteam. I wonder where Croteam got it from.

If you study the mathematics of PGS, you can find that they are equivalent to the sequential impulse algorithm of Mirtich. In the sequential impulse approach, the sparsity and auxilliary variables are implicit. The impulse method has low/high limits for friction so that's covered too.
Antonio Martini
Posts: 135
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London
Erin Catto wrote: If you study the mathematics of PGS, you can find that they are equivalent to the sequential impulse algorithm of Mirtich. In the sequential impulse approach, the sparsity and auxilliary variables are implicit. The impulse method has low/high limits for friction so that's covered too.
please correct me if im wrong, in PGS the contact constraints are solved at the same istant of time, in Mirtich after an impulse has been applied, bodies are numerically integrated, collision detection is performed again etc. So the system is stepped in time while solved. Are you may be referring to the propagation method mentioned in Baraff's and in the more recent Kenny Erleben's book or in code the Bullet solver?

http://www.cs.cmu.edu/~baraff/papers/sig89.pdf

if they are the same it would be interesting to see whats the equivalent of general preconditioning, warm starting and overrelaxation.
Last edited by Antonio Martini on Thu Jan 05, 2006 4:44 pm, edited 1 time in total.
Erin Catto
Posts: 324
Joined: Fri Jul 01, 2005 5:29 am
Location: Irvine
Contact:
Right, I forgot that Mirtich moves time forward after each impulse.

But Nonconvex Rigid Bodies with Stacking http://graphics.stanford.edu/papers/rigid_bodies-sig03/ uses sequential impulses which is equivalent to PGS, and that paper is from 2003.
kenny
Posts: 34
Joined: Thu Jun 30, 2005 8:49 am
Location: Denmark
Contact:
There is a subtle difference in the approaches... It dependens on how the right-hand side vector b in the complementarity formulation is set up.

I believe Projected Gauss Seidel as done by Moreau and Jeans Sweeping Process (blocked non-smooth gauss seidel) is exactly the same as the ``propagation method'', this is because the right-hand vector in the sub-problem is re-evaluated each time... this is the same as would happen in a propagation method.

However, in the non-blocked approach used by nearly everybody, one computes the b-vector at a single instant in time. Therefore one can observe a difference in simulation results, if you don't believe me take a look at these comparisons that I did sometime ago

http://www.diku.dk/~kenny/visionday05.pdf

Of course the Guendelmann et al. approach is not allowed to run until convergence, but neither is the PGS solver:-)

/Kenny
Dirk Gregorius
Posts: 875
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA
Guendelman - Mirtich:
The major difference - at least in my understanding - is that Guendelman doesn't need any velocity threshold since he uses the new integration scheme:
- Candidates
- Collision
- Velocity Update
- Contact
- Position Update
What I also like is, that with this method you don't need the exact TOI.

Impulse Based Simulation vs. Velocity Constraints (LCP - here PGS)
There is a big difference between them and I don't agree here. In the impulse based formulation you seek a sequence of impulses such that the relative normal velocity at each contact point becomes zero. Here you solve one contact point after the other re-evolving the velocity after each contact point as opposed to the PGS where you solve simultan. The iterative nature of the PGS makes them similar, but the idea is totally different.

Bullet - Guendelman
A very good idea of Guendelman (and probably a difference to Bullet) is that instead of applying the full impulse at each contact point he uses negative coefficients of restitution which physically is equivalent to slowing down the object instead of suddenly stop it. It also gives you a better sampling. What also makes a difference is using the current position for the impulse, but still using candidates for finding contacts. I also apply the impulses at the contact with highest penetrating velocity first, what is the same to Guendelman who applies the impulse at the contact with the deepest penetration. Since I don't refresh my contact information I use this heuristic. See:

K. Erleben book: Chapter 5 "Impulse Based Dynamics"
Chatterjee: http://ruina.tam.cornell.edu/research/t ... d_body.pdf
Guendelman: PowerPoint- Contact processing ( Improving accuracy or so )

On average impulses:
In order to avoid the awful tilting which happens through sequential applying impulses (s. Kenny) I use an averaged impulse in the collision phase. Instead of using the arithmetic mean of the contact points I use a weighted sum scaled by the relative normal velocity - only considering points with penetrating velocity. John Schultz at GameDev.net uses something similar. Instead of using the relative normal velocity he scales by the normal impulse (length) at each contact point.
Both schemes give good result for well defined stacking, but have problems for twisted and offset stacks. Also note that you totally miss the friction when the upper box only rotates around the contact normal. The averaging will cancel out the "Torsion Momentum". This is work in progress here...

@Kenny:

(1) Novodex uses Block LCP (at least for the joints). I have a paper of the ETHZ where Adam is so nice to explain quite some details. IIRC they have removed the paper, so if you can't find it I can send it Erwin and he maybe makes it accessable to a bigger audience...

(2) When relaxing the contacts in the Guendelman approach instead of correcting the whole relative velocity at each contact point directly Guendelman starts with smaller impulses (s. above). Do you think that this could be carried over to the LCP formulation.

It can be shown that a SOR paramater must lie between 0 < w_sor < 2 for convergence:

- How do I coose this parameter ? The ODE uses 1.3 - why?
- Wouldn't make under-relaxation make sense here?
- What about making the relaxation parameter be a function of the iteration count, e.g. w1 = 1.0, w2 = 1.0, w3 = 0.9, ... Or maybe ascending, e.g. w1 = 0.1, w2 = 0.2, w3 = 0.3 . Is there something like "adaptive" SUR/SOR around? I think of two opposing constraints here, where once cancels out the other. When you use smaller w parameter later in the solver, you maybe could smear the error a bit this way...?

Regards,

-Dirk

HTH,
-Dirk
Erin Catto
Posts: 324
Joined: Fri Jul 01, 2005 5:29 am
Location: Irvine
Contact:
Impulse Based Simulation vs. Velocity Constraints (LCP - here PGS)
There is a big difference between them and I don't agree here. In the impulse based formulation you seek a sequence of impulses such that the relative normal velocity at each contact point becomes zero. Here you solve one contact point after the other re-evolving the velocity after each contact point as opposed to the PGS where you solve simultan. The iterative nature of the PGS makes them similar, but the idea is totally different.
If you study PGS, you can discover that what it does is solve one constraint at a time and update the net force on the two bodies (assuming only two-body constraints). Since PGS works at the velocity level, this is equivalent to applying an impulse and updating the current velocity of the bodies. This is the mathematical similarity.

Now PGS and sequential impulses may use different impulses based on different friction models, position error resolution, and so forth. But the basic idea of solving one constraint at a time and propagating the result is identical.

If you consider a frictionless system with inelastic contacts, the two algorithms give identical results assuming collision and position updates are handled identically. Of course, the order of collision processing and position updates are independent of how velocities are updated.
Dirk Gregorius
Posts: 875
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA
I see - you mean since you update the lambda (impulse) the velocity change is implicit. Thanks for pointing this out!

@Erin and Kenny:
Did you try different relaxation methods, like SOR or SUR.

And maybe another question. Would it be possible to solve the contact situation the "classical" way, while solving the joints in a block wise fashion? Here is a quote from a Novodex presentation at the ETHZ:
? Because of strong interdependency of
Jacobian rows, pure GS has problems.
? Use blocked GS for joints.
How could this be achieved?

Then two questions on setting up the Jacobian:

1) Would I gain an advantage when setting up the first friction direction parallel to the tangential direction of the tangential velocity at the contact point instead of choosing two arbitrary perpendicular vectors in the contact plane?

2) For long lever arms the rotation term gets predominant. Are there more clever ways to setup the joint Jacobians such that they are more PGS friendly? I mean are there some consideration on choosing the basis for the joint constraint? When one body is fixed I see that I should choose one basis vector in the direction of the lever arm, but how should I choose the joint basis for two dynamic bodies?
? Choice of Basis for constraint
? Rotation term predominant:
dv = n*j/m; dw = j*inv(I)*r x n, r being lever arm.
? If both objects are movable then we can?t do #2

-Dirk
kenny
Posts: 34
Joined: Thu Jun 30, 2005 8:49 am
Location: Denmark
Contact:
DonDickieD wrote:Did you try different relaxation methods, like SOR or SUR.

Yes, but for my purpose it was not worth the trouble.
DonDickieD wrote: And maybe another question. Would it be possible to solve the contact situation the "classical" way, while solving the joints in a block wise fashion?
I am not sure if I entirely understand your question/problem?

/Kenny
Erin Catto
Posts: 324
Joined: Fri Jul 01, 2005 5:29 am
Location: Irvine
Contact:
Here's a paper that makes a similar claim about the Gauss-Seidel nature of sequential impulses:

http://www.ep.liu.se/ecp/010/004/ecp01004.pdf

I find it interesting that with PGS it is viewed as a clever optimization to use an accumulator, but with sequential impulses the optimization is built-in.
Dirk Gregorius
Posts: 875
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA
Thanks for the link Erin. I found the paper very interessting.

I recently looked at the Doom3 SDK source and they seem to have an impulse based formulation as well as a LCP solver. I have not enough knowledge or experience with LCP solvers to recognize the kind of the solver. Do they use a GS or similar method or is this already a more advanced solver that is superior to the GS in the context of computer games?

-Dirk
Erwin Coumans
Posts: 4235
Joined: Sun Jun 26, 2005 6:43 pm
Location: California, USA
Contact:
I recently looked at the Doom3 SDK source and they seem to have an impulse based formulation as well as a LCP solver. I have not enough knowledge or experience with LCP solvers to recognize the kind of the solver. Do they use a GS or similar method or is this already a more advanced solver that is superior to the GS in the context of computer games?
A few years ago, Jan Paul van Waveren (idSoftware's Doom3 physics programmer) told me he used impulse-based method for contacts and an LCP solver for ragdolls (non-contact joints). Probably he uses a direct /pivoting LCP method.

I'll ask if he can explain here on the forum.

It is typical to mix a direct LCP solver with iterative method. Especially high speed Ragdolls don't behave well on impact, their limbs detach. Typical scenario is a high speed car against a pedestrian, or a character standing next to a barrel (explosion). I've seen this in a game I worked on in the past, and we used Featherstone for the ragdolls, and impulse based for contacts.

Erwin
Dirk Gregorius
Posts: 875
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA
A few years ago, Jan Paul van Waveren (idSoftware's Doom3 physics programmer) told me he used impulse-based method for contacts and an LCP solver for ragdolls (non-contact joints). Probably he uses a direct /pivoting LCP method.

I'll ask if he can explain here on the forum.
Good idea - this would be very interessting.
It is typical to mix a direct LCP solver with iterative method. Especially high speed Ragdolls don't behave well on impact, their limbs detach. Typical scenario is a high speed car against a pedestrian, or a character standing next to a barrel (explosion). I've seen this in a game I worked on in the past, and we used Featherstone for the ragdolls, and impulse based for contacts.
Oh yes, we have exactly this problem. When a dragon takes a knight into his mouth and throws him around the ragdoll is always torn apart

What exactly do you mean by "mixing"? Are the non-articulated bodies handles in an iterative solver, while the articulated bodies are solved by some kind of LCP solver (PGS or whatever).

Regarding Featherstone:
How is the collison impulse propagated through the body? While maybe not accurate, the impulse could maybe only be applied to colliding body of the hierarchy. While physical not correct, this could maybe look good enough in a game.

Here is a paper I once spotted, but never read actually - silly me . Maybe interesting in this context:

http://www.gdconf.com/conference/archiv ... ngelis.pdf

-Dirk
Erin Catto
Posts: 324
Joined: Fri Jul 01, 2005 5:29 am
Location: Irvine
Contact:
Conceptually the idea of combining iteration on contacts with joint coordinates for bilateral constraints is appealling. But I wonder how well such a system would perform.

There is nothing that says the ragdoll bone has to get full control of the animation bone. You can just apply the rotation if you like.

BTW, I have worked on such a system and the external impulse just propagates through the joint coordinates. This can be done in O(n) time using the Featherstone mass matrix factorization.
Antonio Martini
Posts: 135
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London
regarding the optimized algorithm of section 4.3 of the above paper, it should be noticed that if we use the propagation method for solving the contacts the algorithm simplifies even further. We would need to know only the effect of a contact on itself informally speaking and so the last loop(but also the first) will not be from "1 to n" but rather "1 to l" if i have understood the paper's notation. We would need to traverse only the subchain starting from the body where the impulse is applied to the root node.
If PGS and the propagation methods are "equivalent", here the latter is a clear winner when it come to down to integrating the Featherstone method, given that PGS need to know the effect that the impulse applied at one contact point has on all the other contacts rather than just one.
However if by 'equivalent' we don't mean 'similar' i assume that the results for one method(preconditioning,overrelaxation,warm starting,contacts handling) transfer automatically to the other one.
Last edited by Antonio Martini on Mon Jan 09, 2006 7:42 pm, edited 6 times in total.
MrElusive
Posts: 2
Joined: Mon Jun 27, 2005 10:11 am
For articulated figures (like ragdolls) DOOM III uses a direct O(n) solver for tree-like structures with bilateral constraints. Furthermore a direct Box Constrained Mixed LCP solver based on an incremental LDLt decomposition is used for all other constraints (including contacts). The advantages of the separate solver for tree-like structures with bilateral constraints are better performance but most of all the bilateral constraints are accurately enforced without being pulled apart due to typical numerical issues of other constraints (like too many or conflicting contacts). All the matrix math is SIMD optimized and numerous tricks are used to save work in the direct solvers. The complete source code is available in the DOOM III and the Quake 4 SDK (http://www.iddevnet.com/).

For simple rigid bodies (like boxes and barrels) DOOM III uses a separate impulse based solver. This solver is rather inaccurate but also very cheap and only solves one contact per object every frame (at 60 fps). Objects sometimes feel a little sticky because so few impulses are applied over time. However, in combination with the continuous collision detection system you can still properly stack hundreds of objects without real problems.

JP,