Hi,
I currently have a physics engine that deals with contacts using regular sequential impulses prioritized by the penetration depth and closing velocity. It also has a friction model.
Such a system is very fast of course, but suffers from vibrating (jittering) bodies while resting, so I found Box2d_Lite of Erin Catto, which uses a modified version of the sequential impulses method, but with impulse accumulation, coherence of contacts using the features and some position correction method.
My problem is, I can't get to the exact idea of Erin Catto. I looked at his GDC2006 presentation and Box2d_Lite, but I still can't figure out what to modify/add to my system.
If you have more resources regarding this method it would be very helpful
Thanks,
Atef

**0**

# Erin Catto's Sequential Impulses with accumulation and coherence

Started by AHashem, Jul 27 2009 06:40 AM

12 replies to this topic

Sponsor:

###
#2
Members - Reputation: **501**

Posted 27 July 2009 - 06:54 AM

One thing he mentions on his site which I found to work quite well is to form the MLCP of the full contact manifold per collision pair and solve it explicitly.

Basically, suppose you are in 2D and you have two boxes stacking. You form the full Jacobian of this system as:

[ -n

[ -n

(That is, it's a single large rectangular matrix).

Then you multiply out J * M

If they're both > 0, you're done. If one is < 0, clamp it to 0 and remove it from the system and solve the resulting 1x1 matrix for the other contact in isolation. If

You can likewise expand this to include friction, but it's more complicated (still working through how it would work myself). Either way, it gives good stacking behavior and far less jitter.

Basically, suppose you are in 2D and you have two boxes stacking. You form the full Jacobian of this system as:

[ -n

_{1}, -n_{1}x r^{AP}_{1}, n_{1}, n_{1}x r^{BP}_{1}][ -n

_{2}, -n_{2}x r^{AP}_{2}, n_{2}, n_{2}x r^{BP}_{2}](That is, it's a single large rectangular matrix).

Then you multiply out J * M

^{-1}* J^{T}to form a 2x2 matrix. Then you explicitly solve this 2x2 matrix to find your two contact impulses.If they're both > 0, you're done. If one is < 0, clamp it to 0 and remove it from the system and solve the resulting 1x1 matrix for the other contact in isolation. If

*that*one is less than 0, then clamp that one as well and you know that there are no impulses involved.You can likewise expand this to include friction, but it's more complicated (still working through how it would work myself). Either way, it gives good stacking behavior and far less jitter.

Darwinbots - Artificial life simulation

###
#3
Members - Reputation: **134**

Posted 27 July 2009 - 07:29 AM

Thanks for your reply!

So, do you suggest any resources on how to solve that MLCP? I believe it should be somehow easier because it's a 2x2 matrix, but I never implemented one myself?

I just started considering handling constraints and contacts using a Jacobian and LCP (or MLCP), but I haven't implemented one yet, so any tutorials, papers or source code would be very helpful, but please for the beginner's level in the subject...

Thanks again,

Atef

So, do you suggest any resources on how to solve that MLCP? I believe it should be somehow easier because it's a 2x2 matrix, but I never implemented one myself?

I just started considering handling constraints and contacts using a Jacobian and LCP (or MLCP), but I haven't implemented one yet, so any tutorials, papers or source code would be very helpful, but please for the beginner's level in the subject...

Thanks again,

Atef

###
#4
Members - Reputation: **501**

Posted 27 July 2009 - 07:46 AM

Ignore for the moment that it's an LCP. You have a 2x2 system of equations. You can solve it using cramer's rule. You'll get two solution impulses: one for each contact point.

Now remember that it's an LCP, so the answers must all be >0. Check your two impulses. If they're both > 0, you're done. If one of them is <0, you need to feed your problem into an LCP solver. Since the system is so small, you can basically brute force it (direct enumeration method).

So you take the impulse that was <0, and you set it to 0 (using the complementarity condition, you're basically setting a slack variable to some non zero, positive value instead). Then you solve for the remaining impulse term as if it were the only contact. If it's <0, you clamp it to 0 as well.

...

Here's a web book on LCP. If you understand section 1.1 you'll be in good shape. It's in postscript format, so you'll need to download and install ghostscript and ghostview. If the download links are broken (they are sometimes), look through the download mirrors they have listed.

But here's the crib notes for it: In LCP format, basically you have Iw - Mz = q. M is your 2x2 matrix that you found from the Jacobian. q is the vector of initial conditions (ie: each element is (1+e)*(delta V * n)). I is the identity matrix, and w is a "slack" variable. If you remove Iw from the equation, you basically have your system of equations Mz = -q. Which you can solve normally using Cramer's rule. (The M matrix here is your J * M

The LCP has a non negativity condition (which is what you want in this case for collision impulses). So elements in z must be >0. If any are <0, they get clamped to 0.

The LCP also has what's called a complementarity condition. Which means that if an element in z is 0, the corresponding element in w is >= 0. So when you have to clamp an element in z to 0, you increase the corresponding element in w so the system will still solve correctly.

Solving the system is basically a combinatorial problem. You want to find what elements in z can be >0, and which need to be clamped to 0 and have the slack vector w handle it.

At the end of the process, you just discard w and the z vector are your impulses. For larger problems, you have to use more sophisticated solvers than direct enumeration, since it has an exponential order of complexity. But since this is a simple 2x2 system, you can do the algorithm I outlined at the beginning of the post: solve the system normally, and if any of the impulses are <0, clamp it to 0 and solve the remaining contact point in isolation.

Using this method in my little homebrew engine I was able to stack over 200 boxes in real time (~55-60 FPS) and have them be stable. You have to add friction to keep them from slipping out laterally, though (imagine stacking wet soap). And you need Baumgarte stabilization, too, or the stack will start to sag (imagine stacking blocks of depleted uranium on top of blocks of styrofoam). But that was it in my engine. No sequential impulses, etc.

I'm not sure if you need to add friction to the LCP, though, or if you can just sort of fake it. I tried adding it to my LCP, but I had to hack it to get it to work right, and even then it only works when the two bodies have a two point contact manifold. If it's vertex vs. edge it won't apply any friction for some reason.

Now remember that it's an LCP, so the answers must all be >0. Check your two impulses. If they're both > 0, you're done. If one of them is <0, you need to feed your problem into an LCP solver. Since the system is so small, you can basically brute force it (direct enumeration method).

So you take the impulse that was <0, and you set it to 0 (using the complementarity condition, you're basically setting a slack variable to some non zero, positive value instead). Then you solve for the remaining impulse term as if it were the only contact. If it's <0, you clamp it to 0 as well.

...

Here's a web book on LCP. If you understand section 1.1 you'll be in good shape. It's in postscript format, so you'll need to download and install ghostscript and ghostview. If the download links are broken (they are sometimes), look through the download mirrors they have listed.

But here's the crib notes for it: In LCP format, basically you have Iw - Mz = q. M is your 2x2 matrix that you found from the Jacobian. q is the vector of initial conditions (ie: each element is (1+e)*(delta V * n)). I is the identity matrix, and w is a "slack" variable. If you remove Iw from the equation, you basically have your system of equations Mz = -q. Which you can solve normally using Cramer's rule. (The M matrix here is your J * M

^{-1}* J^{T}matrix).The LCP has a non negativity condition (which is what you want in this case for collision impulses). So elements in z must be >0. If any are <0, they get clamped to 0.

The LCP also has what's called a complementarity condition. Which means that if an element in z is 0, the corresponding element in w is >= 0. So when you have to clamp an element in z to 0, you increase the corresponding element in w so the system will still solve correctly.

Solving the system is basically a combinatorial problem. You want to find what elements in z can be >0, and which need to be clamped to 0 and have the slack vector w handle it.

At the end of the process, you just discard w and the z vector are your impulses. For larger problems, you have to use more sophisticated solvers than direct enumeration, since it has an exponential order of complexity. But since this is a simple 2x2 system, you can do the algorithm I outlined at the beginning of the post: solve the system normally, and if any of the impulses are <0, clamp it to 0 and solve the remaining contact point in isolation.

Using this method in my little homebrew engine I was able to stack over 200 boxes in real time (~55-60 FPS) and have them be stable. You have to add friction to keep them from slipping out laterally, though (imagine stacking wet soap). And you need Baumgarte stabilization, too, or the stack will start to sag (imagine stacking blocks of depleted uranium on top of blocks of styrofoam). But that was it in my engine. No sequential impulses, etc.

I'm not sure if you need to add friction to the LCP, though, or if you can just sort of fake it. I tried adding it to my LCP, but I had to hack it to get it to work right, and even then it only works when the two bodies have a two point contact manifold. If it's vertex vs. edge it won't apply any friction for some reason.

Darwinbots - Artificial life simulation

###
#5
Members - Reputation: **134**

Posted 27 July 2009 - 07:55 AM

Thanks ALOT for your helpful reply! :)

I'll go through with what you said and hopefully get the engine working properly..

About *faking* friction..could you just elaborate more on that?

Thanks again..I was stuck for three days trying to get my head round this...:)

Atef

I'll go through with what you said and hopefully get the engine working properly..

About *faking* friction..could you just elaborate more on that?

Thanks again..I was stuck for three days trying to get my head round this...:)

Atef

###
#7
Members - Reputation: **501**

Posted 27 July 2009 - 08:15 AM

When I first thought about friction, I was going to form a friction contact at each real contact (basically form the perpendicular to each contact normal for each contact point, and add that constraint). So for a two point contact manifold that would give 4 constraints. I would then use my matrix library based on SIMD and invert the matrix to solve the system.

However, if you try doing that you quickly run into a problem I didn't understand at the time: the matrix will be singular. The reason is actually sort of interesting.

So imagine you're in 2D. There are three degrees of freedom per body. So imagine now you have two bodies. There are six degrees of freedom for the system. Now add a constraint between the bodies. The constraint, no matter what it is, removes a single degree of freedom from both bodies. So each (binary) constraint removes two degrees of freedom from the system.

So in the 2D friction case I was working on, I was trying to add 4 constraints, which would remove a total of 8 degrees of freedom. But there's only 6 degrees of freedom total, so the system ends up being over constrained and no solution exists.

In my engine I just dropped one of the friction constraints arbitrarily. And if any of the contact impulses were <0, I clamped all of them to 0. A really terrible hack, but it worked well enough.

If you're building on my work, there are a couple of ways I can think to approach the problem:

1. First, solve the 2x2 LCP without friction, just from the contact impulses. Then solve the 2x2 friction linear system (without any constraints. Pretend you have infinite friction coefficients) for your friction impulses. Then clamp the friction impulses based on the values you got from the contact impulses. So if(contact impulse * friction coefficient < abs(friction impulse)) -> clamp friction impulse so its absolute value is equal to contact impulse * friction coefficient. It's not mathematically correct to handle it this way, but it might work just fine in practice.

2. Take the two constraints from the two point contact manifold, and determine a new constraint which removes the last two degrees of freedom. This is your implicit friction term. You can now solve the critically constrained (by which I mean you've removed exactly all degrees of freedom from the system) system. You then enforce the constraints that the two contact impulses must be >0, and the friction impulse must be bounded by +/- (coefficient friction * sum of contact impulses).

Unfortunately I have no idea how to solve this problem. It's not a LCP. I don't even think it's a MLCP (mixed linear complimentarity problem). It might be something individual solvers can handle (Lemke, Simplex, etc.), though.

If you graph the problem, the "feasible region" (the area in which a solution might exist given the constraints) is the intersection of cones and half planes, so it's very well defined. I'm thinking you might be able to do something like "find the closest point in the feasible region to what the solution would be if it were unconstrained", which is a least squares quadratic programming problem that's fairly well defined. You might even be able to use a modified GJK algorithm for it. I'm not really sure to be honest.

3. The ultimate in faking it: just take your sequential impulse solver engine exactly as is, and just replace solving the contact impulses in a manifold sequentially with the 2x2 LCP solver. Your friction terms are found sequentially just like normal. You'd probably be better off solving the friction term as if there's only one per contact manifold though (centroid (center) of the contact manifold).

...

Beyond that I have no idea.

[Edited by - Numsgil on July 27, 2009 3:15:46 PM]

However, if you try doing that you quickly run into a problem I didn't understand at the time: the matrix will be singular. The reason is actually sort of interesting.

So imagine you're in 2D. There are three degrees of freedom per body. So imagine now you have two bodies. There are six degrees of freedom for the system. Now add a constraint between the bodies. The constraint, no matter what it is, removes a single degree of freedom from both bodies. So each (binary) constraint removes two degrees of freedom from the system.

So in the 2D friction case I was working on, I was trying to add 4 constraints, which would remove a total of 8 degrees of freedom. But there's only 6 degrees of freedom total, so the system ends up being over constrained and no solution exists.

In my engine I just dropped one of the friction constraints arbitrarily. And if any of the contact impulses were <0, I clamped all of them to 0. A really terrible hack, but it worked well enough.

If you're building on my work, there are a couple of ways I can think to approach the problem:

1. First, solve the 2x2 LCP without friction, just from the contact impulses. Then solve the 2x2 friction linear system (without any constraints. Pretend you have infinite friction coefficients) for your friction impulses. Then clamp the friction impulses based on the values you got from the contact impulses. So if(contact impulse * friction coefficient < abs(friction impulse)) -> clamp friction impulse so its absolute value is equal to contact impulse * friction coefficient. It's not mathematically correct to handle it this way, but it might work just fine in practice.

2. Take the two constraints from the two point contact manifold, and determine a new constraint which removes the last two degrees of freedom. This is your implicit friction term. You can now solve the critically constrained (by which I mean you've removed exactly all degrees of freedom from the system) system. You then enforce the constraints that the two contact impulses must be >0, and the friction impulse must be bounded by +/- (coefficient friction * sum of contact impulses).

Unfortunately I have no idea how to solve this problem. It's not a LCP. I don't even think it's a MLCP (mixed linear complimentarity problem). It might be something individual solvers can handle (Lemke, Simplex, etc.), though.

If you graph the problem, the "feasible region" (the area in which a solution might exist given the constraints) is the intersection of cones and half planes, so it's very well defined. I'm thinking you might be able to do something like "find the closest point in the feasible region to what the solution would be if it were unconstrained", which is a least squares quadratic programming problem that's fairly well defined. You might even be able to use a modified GJK algorithm for it. I'm not really sure to be honest.

3. The ultimate in faking it: just take your sequential impulse solver engine exactly as is, and just replace solving the contact impulses in a manifold sequentially with the 2x2 LCP solver. Your friction terms are found sequentially just like normal. You'd probably be better off solving the friction term as if there's only one per contact manifold though (centroid (center) of the contact manifold).

...

Beyond that I have no idea.

[Edited by - Numsgil on July 27, 2009 3:15:46 PM]

Darwinbots - Artificial life simulation

###
#8
Members - Reputation: **501**

Posted 27 July 2009 - 08:18 AM

Quote:

Original post by AHashem

I'd also like to know your opinion in accumlating the impulses and using the features of collision to store collision points among frames...do you do this in your engine??

No, but my engine was pretty ghetto :) Basically I was applying impulses as soon as I found them instead of storing them across iterations. I wasn't doing any warm starting, either. The only caching I was doing was only performing collision detection once per frame.

You'll get better results (probably) with both methods.

Darwinbots - Artificial life simulation

###
#11
Members - Reputation: **501**

Posted 27 July 2009 - 02:07 PM

You

So in practical use cases, there's probably a sweet spot. For a hinge constraint, you probably just end up replacing one of the points in the contact manifold with a hinge, and so the system is still 2x2.

Doing what I've suggested, solving the contact manifold between bodies iteratively, is sort of like an "iterative kernel solver". It might be called something else, I'm sort of grasping at terminology here. But basically you can increase the size of the "kernel" (increase the size of the system you solve at each step), and it'll decrease the number of iterations you need to arrive at an answer, but increase the number of operations you need to perform per kernel solve. So in a theoretic sense there's some sweet spot between kernel size and iterations, but I don't have any idea how to approach it.

For certain problems, where constraints can be arranged into a tree structure (ragdolls, for instance), there's a way to solve the system in linear time. Clicky. And incidentally, I'm currently working my way through the math to find a global solver for constraints in linear time. You can follow my progress in this thread. But neither Baraff's paper or my work in the latter part of that thread deal with how to handle box constraints on the values (ie: constrain the impulses to the range [-10, 20] or something like that). Which is another open question I can't figure out and doesn't seem to have been adequately addressed anywhere that I can tell.

*could*. In theory you can form the Jacobian for your entire physics world and solve it all directly that way. The problem is that in practice the computational complexity increases to be more than you can do in the 16ms you have to do a frame in. That's not even counting the combinatorial part of the LCP, which can potentially be exponential if you use a stupid solver (direct enumeration), or the problem is ill formed in ways I don't entirely understand.So in practical use cases, there's probably a sweet spot. For a hinge constraint, you probably just end up replacing one of the points in the contact manifold with a hinge, and so the system is still 2x2.

Doing what I've suggested, solving the contact manifold between bodies iteratively, is sort of like an "iterative kernel solver". It might be called something else, I'm sort of grasping at terminology here. But basically you can increase the size of the "kernel" (increase the size of the system you solve at each step), and it'll decrease the number of iterations you need to arrive at an answer, but increase the number of operations you need to perform per kernel solve. So in a theoretic sense there's some sweet spot between kernel size and iterations, but I don't have any idea how to approach it.

For certain problems, where constraints can be arranged into a tree structure (ragdolls, for instance), there's a way to solve the system in linear time. Clicky. And incidentally, I'm currently working my way through the math to find a global solver for constraints in linear time. You can follow my progress in this thread. But neither Baraff's paper or my work in the latter part of that thread deal with how to handle box constraints on the values (ie: constrain the impulses to the range [-10, 20] or something like that). Which is another open question I can't figure out and doesn't seem to have been adequately addressed anywhere that I can tell.

Darwinbots - Artificial life simulation

###
#12
Members - Reputation: **665**

Posted 29 July 2009 - 04:23 AM

Quote:

Original post by Numsgil

So imagine you're in 2D. There are three degrees of freedom per body. So imagine now you have two bodies. There are six degrees of freedom for the system. Now add a constraint between the bodies. The constraint, no matter what it is, removes a single degree of freedom from both bodies. So each (binary) constraint removes two degrees of freedom from the system.

I think that this may be incorrect -- a single normal/axial constraint (constrains the distance between a point on each body projected onto some axis; used in collision constraints) should only remove a single DOF from the system, i.e it removes one relative DOF, not one from each body.

Example: describe body B relative to body A; B has 3 DOF relative to A.

Constraining the position of B relative to A along some axis only removes one DOF from the whole system -- B is still free to move perpendicular to n, and to rotate about pB, and A is free to move and rotate in any way.

A "pin"/joint constraint does remove 2 DOF, but that's because it fixes the position of one body relative to the other -- analogous to two axial constraints.

###
#13
Members - Reputation: **501**

Posted 29 July 2009 - 04:40 AM

Hmm, looking it over again I think you're right. I musta got confused somewhere. I guess that doesn't explain why my 4x4 LCP with friction was singular.

Darwinbots - Artificial life simulation