Impulse Solvers and Guaranteed Collision Free Frames

D.V.D said:
It seems odd to me that the math formulation wouldn't guarantee that we would have a collision free frame even if we solved our system of constraints perfectly, since the newly applied impulses can create new constraints that we didnt have before.

Assume we add a collision test for each solver iteration, we still have cases of potential intersections. E.g. If we put many small boxes into one large hollow box, but the summed volume of the small boxes is larger than the volume of the empty space in the hollow box. Or some ragdoll has it's lower leg in front of a thin wall, but it's upper leg on the back. The joint at the knee might pull the bodies together stronger than the contact separates them. We can construct many examples where things will go wrong, and it's not possible in practice to prevent such cases from happening.

Many of those cases can be explained by the artificial idea of ‘rigid’ bodies. They can't break, they can't compress. Such thing does not exist in reality, so we can't expect to get realistic results in any case from this model.
My goal here would be the simulation converging at a resting state, but allowing penetrations. The bodies should not jump around or jitter like hell even if constraints remain violated.

D.V.D said:
a simple 2d impulse solver

That's the problem. There are good solutions to highly constrained systems of contacts, but they're not simple. Look up “LCP solver”.

On detecting object overlap, you can cut the time step in half and retry, which can get you out of the kind of situations you describe. But not in fixed time. In some situations the physics frame rate may have to drop.

Often this isn't a problem, if you're just blowing up stuff. You can tolerate some bogus interpenetration.

@JoeJ, so I understand that we can construct such cases like a box that contains rigid bodies but doesn't have room for them, but if we ignore such cases, we still get problems in environments where we do have a infinite amount of space to work with. Like when you just have a bunch of bodies fall, then there is enough space for all of them every frame but the “perfect” solution that solves all the constraints in a given frame isn't actually guaranteed to be intersection free.

@Nagle For the LCP solver, does it rely on us specifying the constraints and then finding a solution to that? My reason for asking is that to me it feels like the problem specification is the issue here, not the PGS solver.

I guess my overall question is, why is the problem specified in a mathematical sense not actually a perfect solution? For example, in graphics we have the rendering equation which is the “perfect” solution to lighting and GI but we know we can't simulate it directly so we approximate it. For physics, I thought that the perfect solution would be the constraint vectors/matrices and then us solving the system of equations this generates so that all our constraints for a given frame are solved. This does not appear to actually be the case since solving this system perfectly doesn't result in a collision free solution. Is this a accurate description of the situation? Is there a way to specify the problem mathematically where this isn't the problem? For reference, my understanding of how the problem is specified comes from this post: https://www.toptal.com/game/video-game-physics-part-iii-constrained-rigid-body-simulation

IMO the difference is that the rendering equation doesn't have as much “feedback” in the calculations: the light might bounce around a bunch, but it's not moving anything when it does so, which makes the problem easier.

AFAIK there's no known way to solve more than a single simultaneous collision perfectly. See this article for a great overview of the problem, and some imperfect solutions: https://www.myphysicslab.com/engine2D/collision-methods-en.html

Not to mention, the whole premise of impulse-based solver is you're using linear approximations of the constraints -- that means even if you *did* perfectly solve the system of impulses, you'll inevitably have penetration/errors because of the inherent errors due to the approximation. And finally, the type of solvers typically used to find the impulses are themselves iterative and approximate by nature! So IME it's better to find method(s) you can use to correct the error/penetration, trying to completely avoid it in the context of a game sim is very hard.

(I think Doom3 used a "never let things overlap" approach, however the fact that this wasn't widely adopted by Bullet/Havok/etc. suggests that there are some bad tradeoffs to such an approach.)

See this thread for some discussion on handling-vs-avoiding penetration: https://pybullet.org/Bullet/phpBB3/viewtopic.php?p=460

And this thread about some ways to deal with penetration: https://pybullet.org/Bullet/phpBB3/viewtopic.php?f=4&t=11675&p=39301&hilit=penetration#p39301

D.V.D said:
why is the problem specified in a mathematical sense not actually a perfect solution?

Because impulse/constraint physics has discontinuities. Either you have contact, or you don't. An infinitesimal movement can change the contact situation. So trying to solve the system by numeric hill climbing is not sufficient.

If you do spring/damper physics, with actual contact forces, this can be avoided. You have to have multi-point contact, so you don't get that teetering-between-two-contact points discontinuity when a cube lands on a flat surface. Now you have a system where you can hill-climb along the gradient until the error becomes infinitesimal and the objects are in a valid state. The compute load is higher, you have to cut the step size to tiny values at the beginning of a collision, and so it's hard to do in real time.

Thus, you can have accurate physics or constant-time physics, but not both.

Here's my original spring/damper ragdoll, from 1997. This is spring/damper for the contacts, has multiple contact points between pairs of bodies when necessary, and, uses Featherstone's algorithm for the character. The time step is cut until an error measure on the integrator indicates that the time step is short enough that a linear approximation is valid.

Nobody does it this way any more. But it's sound and stable.

This may be the first successful ragdoll demo.

D.V.D said:
For physics, I thought that the perfect solution would be the constraint vectors/matrices and then us solving the system of equations this generates so that all our constraints for a given frame are solved.

That's much harder than your example of the rendering equation. For rendering we can assume our geometry to be static and constant for the duration of our frame.
No longer the case for physics: Here everything changes, including the ‘relationships’ between bodies, e.g. a contact. And we don't know about future relationships, so we can't have a single and compact system of equations to represent all relations. I use the word ‘relationship’ because i forgot the proper word - maybe it was ‘bilateral’ joint, meaning a contact joint which has no effect once bodies separate from each other.

A similar situation is joint limits. A hinge joint may violate it's limit or not, and again we don't know all states our joint might have during the frame.

So we have two options:
1. Order all events (like a collision) sequentially in time, so we can update our system of equations every time it changes. We will get a correct solution. But in case of paradox situations like my boxes example, the solver would get stuck in an infinite loop. And worse: Worst case will have much higher computational cost than average case.
So that's interesting mainly if you can guarantee your number of bodies remains small.
2. Treat all events as simultaneous and lasting over the whole (sub)frame. Gives only an approximation, but we can bound the error by our timestep.

(I just see there's already a similar answer)

@raigan Thanks for the link on resolving collisions, I haven't realized that this was such a hard problem for rigid bodies to be honest and it helped put that in perspective with the examples (seems to be similar to what @nagle was mentioning). I agree that we are better off handling penetrations gracefully rather than enforcing they cannot happen due to computational constraints, I guess a lot of this is more theory for me. But you did mention that even if we solved the system of constraints, you would have penetrations due to the constraints being approximations. My argument is you would have penetrations even with perfect constraints since we are only adding constraints to the matrix for objects that collide, but the resolved positions might still collide with objects they didn't previously collide when we setup the matrix.

@joej has mentioned this as well, but the difference between the rendering equation and collision responses is that our geometry changes as we resolve more collisions (or at least positions and orientations do, so it feedbacks onto itself). I agree that this is the case, I think from @joej ‘s post, my working idea was the second option of treating all events as simultaneous and lasting for the whole frame. If we treat them this way, can we build a system of equations that guarantees its solutions won’t penetrate after it is solved. My working idea on how to do this would be the following;

• Represent all geometry of each object as some sort of equation (maybe a SDF?)
• Add n^2 penetration constraints for each object vs all others such that each object doesn't overlap (we need the equation/SDF here to specify the constraint)
• Solve this matrix

The above can be done for small demos with simple geometry I think pretty easily if our geometry is a circle since then we can specify penetration constraints without a contact point. The part that was confusing me is that the source code on github I was seeing for physics engines don't appear to try to approximate this. Instead, they do the following:

• Find all collisions in the scene
• Specify a constraint for each collision/contact point in our matrix
• Solve the matrix

This was confusing me because solving the above does not guarantee you don't have collisions in a particular frame, even if the matrix is solved perfectly. You would need to represent all your geometry as some sort of math equation so that you can use it to write out a constraint that doesn't rely on having a contact point. So if object A collides with object B, resolving this contact point might make object B now collide with object C which it wasn't colliding with before so its constraint wasn't in our matrix when we solved it. So the physics engines I see on github are sorta “distributing” this computation across frames, and hoping that overtime, we won't have new collisions created as we resolve contact points since things will stabilize. It makes sense to do it this way, but I guess I was surprised this was the case.

Thanks everyone for the info! I hope I correctly reiterated how this all works in the above.

D.V.D said:
So the physics engines I see on github are sorta “distributing” this computation across frames, and hoping that overtime, we won't have new collisions created as we resolve contact points since things will stabilize.

Yeah, and there are some arguments why this should do fine:
Timestep is small, so error is small too.
Chaotic scene (explosion, many collisions, high velocities) might not need an accurate solution - it's just chaos.
Resting scene (stack of boxes) needs an accurate solution and is the main problem for us to solve. For that we get a complete graph of all contacts anyway.

D.V.D said:
Add n^2 penetration constraints for each object vs all others such that each object doesn't overlap (we need the equation/SDF here to specify the constraint)

Assuming this algorithm aims for a more complete collision graph than the above, we could add pairs of potentially colliding bodies so we can keep checking them within each solver iteration.
We could implement it with a speculative broad phase with extended bounds of the bodies. This way we don't have to run the complete collision detection for each solver iteration.
Result should be more accurate, but also more work because there are more potentially colliding pairs than actually colliding pairs, and we need to check for penetrations constantly.
The penetration check remains expensive, even if you use SDF. Because you have to intersect two SDFs for each pair, it's not just a single lookup but a search for the deepest penetration out of many, or some weighted average of all penetrations.
And finally, there is still no guarantee we detect all collisions as the system changes. Things could go out of our extended bounds, which was just a guess we did for the broad phase.

I assume it's no win in most cases. Performance suffers and errors still remain. Due to errors, you need good enough results from an incomplete collision graph anyway, so i would implement the faster algorithm first. Maybe it's good enough and you'll see it's not necessary to handle each and every contact which might happen within one frame. (But that's just my personal belly feeling assuming generic scenes.)

@D.V.D I would recommend looking up Erin Catto's talks about how the Box2D solver works (see https://box2d.org/publications/​ ), and/or you could read the paper “Position-Based Dynamics” (or these slides http://mmacklin.com/EG2015PBD.pdf​ ), to get some good references for why and how most people solve game physics using iterative Gauss-Seidel methods (which amount to solving pairwise constraints and hoping the error propagates/diffuses and things eventually converge on a solution).

There are other more direct solvers, eg IIRC Doom 3 solved the big-matrix problem directly, and Cel Damage used a Conjugate Gradient-based solver.

About “the resolved positions might still collide with objects they didn't previously collide when we setup the matrix.”, in an impulse-based solver this is a problem for next frame (ie those collisions will be found next frame), because during the current frame nothing changes position, only velocities are changed. (This is what lets impulse-based solvers be so fast, the constraint gradients (which are usually a function of position) are fixed so they can be computed once and re-used over multiple solver iterations during the same frame.) Position errors are typically corrected by feeding some of the error into the velocities, ie generating extra forces which push the positions out of error, and/or you can directly solve the position constraints (see below).

In a position-based simulator, it's indeed a problem that the solver can move things into collision when there weren't any collisions initially. The typical solution for this is “speculative contacts”: you generate a collision constraint between anything that might collide this frame (ie generate a constraint if the shapes are within some distance apart, not only if they currently overlap), and then when solving these constraints you skip them if they're not currently overlapping. This moves collision detection into the solver loop a bit, which is more expensive. Also you can get errors depending on how much you approximate the collision constraint (ie typically they're point-to-plane constraints, while the actual shapes that generated them are finite (not infinite planes), which means you can get ghost collisions). This paper has lots of useful references for PBD simulation: http://mmacklin.com/smallsteps.pdf

Anyway

@joej @raigan Thanks for the helpful posts, I think it really cleared up a lot of stuff for me. I understand why we don't do these kind of perfect solvers with each potential collision due to computation cost but it does clear up what the “perfect” solution would look like, and where we choose approximations.

I'm looking to implement a PGS solver right now, but I'm running into the following issue. I'm looking at the code here as a reference: https://github.com/granj2020/Cirobb-Engine/blob/master/cirobb/Manifold.cpp​ ← from my understanding, this is a simplification of Box2DLite. Sequential Impulses and PGS should be the same thing, and so I went through the math of deriving what how we solve the below matrix:

J M^-1 J^t lambda = -(Jq2 + b)

where q2 is the integrated velocities for position/angle that are not yet corrected, and b is a bias. My understanding is that we want to throw all our constraints into J since we can have multiple constraints operating on a single body so we want the matrix to capture that interaction. Going through the algebra, I got a pretty complicated equation that solves for lambda using Gauss-Seidel (decompose the matrix on the left into a upper and lower triangular, then propagate the results backwards to solve for lambda). The equation is pretty complicated, but the main thing that comes out of it is that we have gradients of different constraints dotted with each other. In practice, the gradients of constraints are sparse since they affect only one or two bodies, but this is still a interaction that needs to be captured.

When I look at the source for https://github.com/granj2020/Cirobb-Engine/blob/master/cirobb/Manifold.cpp​ however, it seems to build a separate J matrix for each constraint and deal with all of them separately (non penetration constraints that deal with the same body should be dotted with each other). Now I understand that PGS is all about handling one constraint at a time irrespective of the others, but doesn't that make it not equivalent to Gauss-Seidel method? Or is there a reason why we build J separately for each constraint instead of throwing it all into one big constraint matrix? I ask because the Box2D GDC presentations appear to handle constraints this way, but going through the article here, they mention that you want to throw the constraints all together: https://www.toptal.com/game/video-game-physics-part-iii-constrained-rigid-body-simulation​ and I'm not really sure which is correct now.

The equation for each x term or in our case lambda term in Gauss-Seidel is given here: https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method​ and the A matrix that we use is J M^-1 J^t which does dot products between the gradients of constraints against each other, so if we have more than one constraint, we will have these gradient dot products in our Gauss-Seidel solver but those don't seem to be present in the source code. Not sure if I'm just messing up my math somewhere or not getting something thats simple here.

I honestly don't know the math well enough to answer you, aside from you're correct that Box2D (and, every other physics engine AFAIK) only solves constraints pair-wise. I think “relaxation” is the general term. (Recently Roblox added a direct solver to their simulation, and there are some other direct solvers (eg IIRC PhysX has a Featherstone solver for articulated trees of joints), but in general things are solved row-by-row. The “small steps position-based dynamics” paper I linked above should explain this and have some references to papers where they expand on the details.

If I'm reading it correctly, Wikipedia suggests that solving each row of the constraint matrix is how these methods are meant to work. See the section containing the text “can be computed sequentially using forward substitution” on this page: https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method​

Or, if instead of changing the state when solving each row, you instead accumulate and average the changes, you get this method: https://en.wikipedia.org/wiki/Jacobi_method

(Jacobi converges more slowly, but behaves better in the presence of contradictory/impossible pairs of constraints since it produces a midpoint instead of ping-ponging between the the solutions)

I would recommend reading through the archives of the Bullet forum, there are a bunch of threads which cover physics simulation (in around 2009 when I was learning, I posted a lot of questions and some kind people helped walk me through the answers): https://pybullet.org/Bullet/phpBB3/viewforum.php?f=4