CA and numerical integration and root finding

[Numsgil]

I'm reading the conservative advancement paper that's Mirtich's thesis. I understand that bullet uses CA, so I thought this would be a good place to talk about it.

I'm working on a 2D simulator for an alife project called Darwinbots. My idea predates reading the CA paper, but the two are similar so I'd like to run with it here.

Instead of calculating progressive conservative advancement steps and integrating the entire simulation, every body computes it's exact time of collision whenever it doesn't know its next time of collision. When the time of collision arrives, we solve the resulting island (that's the proper term, right?) and then find the next exact time of impact for each body in the island. This way a collision in one part of the simulation does not affect the bodies in another part (ie: we don't have to integrate the entire simulation when we find the TOI).

I've gotten this more or less working with a simple points vs. static AABB simulation which uses only constant acceleration terms (ie: gravity). I resolve the collision equation down to a quadratic polynomial and solve for the first positive odd root (even roots indicating glancing collisions, which aren't useful). So I think the concept is sound. Basically I avoided any numeric integration and used the actual analytical solutions for motion, and then root find using quadratic, cubic, or even quartic methods for polynomials.

My initial thinking was that I would build up an exact analytical equation for the motion of each body in my sim, and then build a quadratic (or higher) spline against it to approximate the motion, and then use the splines for root finding. However, finding an exact analytical equation when you include things like a spring mesh is just way too hard, meaning my present method doesn't scale very well. Not to mention issues with the decidedly non-polynomial rotational terms.

All of which means I probably need to delve in to numerical integrators and root finders (especially since something like a spring mesh is just too unwieldy to solve analytically). So what I'll end up doing is pretty close to doing conservative advancement on bodies in that I use a numerical root finding method. However I can use any step size and even over-shoot the root so long as I arrive at the first positive, odd root.

Now, it seems to me to be a bit backwards to take my acceleration terms for both colliding bodies, numerically integrate an approximation for position and velocity at various time steps, and use those for numerical root finding (ala most popular numerical root finding methods). Plus I know I'm only interested in the first positive, odd root (or none if there isn't a positive odd root), so maybe that will help me in some way (ie: discard any even roots somehow).

Are there any numerical root finding methods that work primarily from the second derivative? Or maybe I can build my own? I used to have a numerical analysis text book from my college course, but it was lost when I moved. I seem to remember something involving second derivatives, but I can't remember what it was called or how it worked. We also at some point went over how to build a numerical integrator, but I don't think I did very well in that section

Oh, and has this been done before? Is it called something? I've only read a handful of physics simulation papers, so for all I know I'm reinventing the wheel.

[Enki]

The rough simulation method you are describing is essentially a discrete event simulation (DES). In formal discrete event sims each event is entirely independent of another and computed from a time ordered queue which holds those events projected to happen in the future. Well independent at least until collision has been projected. Then the collision event is scheduled with the event objects starting state, and you forget about it physics-wise until it happens. It has been proven that the degenerate case of a DES simulation is a standard time-step simulation, but the paradigm has been used more for analysis modeling rather than physical modeling because the disciplines don't overlap except in a very very few schools. It also takes thinking that state is no longer where you are, but where you started and what equations project you into the future, which isn't hard but its not how physics sims have been thought about for the last few decades. I am working on something like this that takes the idea significantly farther than Mirtich's Timewarp paper and combines it with a Quantized State System DES approach. It's not trivial though, DES motions are subject to more than just collisions.

Nonlinear parameters and changing environmental variables mean you almost always have competing reasons which may require you to recompute the motion before your objects bump. Not to mention doing this one event at a time means you need to have some way to re-recompute if another collision event changes things and creates a new earlier collision that the one you just computed. That's not as bad as it sounds as long as you have a decently efficient DES system. I'm hoping this will avoid a whole bunch of numerical stability issues by doing fewer computations as well as make the collision side faster since I think the the broadphase will work bottom up rather than top down.

Specifically the issues you raise talking about the splines are the kind of thing the QSS stuff should help with. It's pretty hardcore so beware if you go looking for it.

And no I don't think a soup to nuts implementation has been done before, at least I hope not! The Mirtich Timewarp paper seems to be the closest, but that's just a collision layer partially implemented discretely.

[Numsgil]

The rough simulation method you are describing is essentially a discrete event simulation (DES). In formal discrete event sims each event is entirely independent of another and computed from a time ordered queue which holds those events projected to happen in the future. Well independent at least until collision has been projected. Then the collision event is scheduled with the event objects starting state, and you forget about it physics-wise until it happens.

Yep, that does indeed sound like what I'm doing. Good to know I'm not the first person to play with it.

but the paradigm has been used more for analysis modeling rather than physical modeling because the disciplines don't overlap except in a very very few schools.

What is "analysis modeling"?

Not to mention doing this one event at a time means you need to have some way to re-recompute if another collision event changes things and creates a new earlier collision that the one you just computed. That's not as bad as it sounds as long as you have a decently efficient DES system.

Yes, this part of the algorithm is reminiscent of "Fast Local Search" (I think that's what it's called), from minmax algorithms (I think it was coupled a lot in the literature with "Guided Local Search".) I used to have a good paper link, but I can't find it now. Basically the idea was that if a node hadn't changed in the previous iteration, you removed it from the list of nodes to check. If another node interfered with that node you deactivated, you would reactivate and re-add it to the search list. The result was a rather neat cascade effect that provided remarkable speed improvements in pairwise algorithms (I used it in a Traveling Salesman problem for a CS class competition.)

And no I don't think a soup to nuts implementation has been done before, at least I hope not! The Mirtich Timewarp paper seems to be the closest, but that's just a collision layer partially implemented discretely.

Always a little terrifying to be told you're blazing new ground

...Quantized State System DES approach. It's not trivial though, DES motions are subject to more than just collisions.
...
Specifically the issues you raise talking about the splines are the kind of thing the QSS stuff should help with. It's pretty hardcore so beware if you go looking for it.

Hardcore things scare me But I'll go looking for it anyway. Is there a 10 word or less summary version?

[Enki]

The rough simulation method you are describing is essentially a discrete event simulation (DES). In formal discrete event sims each event is entirely independent of another and computed from a time ordered queue which holds those events projected to happen in the future. Well independent at least until collision has been projected. Then the collision event is scheduled with the event objects starting state, and you forget about it physics-wise until it happens.

Yep, that does indeed sound like what I'm doing. Good to know I'm not the first person to play with it.

but the paradigm has been used more for analysis modeling rather than physical modeling because the disciplines don't overlap except in a very very few schools.

What is "analysis modeling"?

Think Qualitative Analysis or Operations Research simulation & modeling, they generally don't care about the physical details, looking at events as "things" that happen in some order. They try to optimize the orderings/scheduling using the simulations as the tests for their optimization methods.

Not to mention doing this one event at a time means you need to have some way to re-recompute if another collision event changes things and creates a new earlier collision that the one you just computed. That's not as bad as it sounds as long as you have a decently efficient DES system.

Yes, this part of the algorithm is reminiscent of "Fast Local Search" (I think that's what it's called), from minmax algorithms (I think it was coupled a lot in the literature with "Guided Local Search".) I used to have a good paper link, but I can't find it now. Basically the idea was that if a node hadn't changed in the previous iteration, you removed it from the list of nodes to check. If another node interfered with that node you deactivated, you would reactivate and re-add it to the search list. The result was a rather neat cascade effect that provided remarkable speed improvements in pairwise algorithms (I used it in a Traveling Salesman problem for a CS class competition.)

Well, DES is a lot different than that, but I think you can get similar benefits coming at the problem from the different direction.

And no I don't think a soup to nuts implementation has been done before, at least I hope not! The Mirtich Timewarp paper seems to be the closest, but that's just a collision layer partially implemented discretely.

Always a little terrifying to be told you're blazing new ground

But helpful if a dissertation is involved!

...Quantized State System DES approach. It's not trivial though, DES motions are subject to more than just collisions.
...
Specifically the issues you raise talking about the splines are the kind of thing the QSS stuff should help with. It's pretty hardcore so beware if you go looking for it.

Hardcore things scare me But I'll go looking for it anyway. Is there a 10 word or less summary version?

"time discretization is replaced by state variables quantization arriving(sic) to discrete event simulation models"

"they can reduce the number of calculations, their parallel implementation is straightforward and they can exploit structural properties like sparsity in a very efficient fashion."

[the paper is: Kofman, Discrete Event Simulation of Hybrid Systems]

[Numsgil]

Okay, I've read through a couple paper on DES and QSS. Certainly a great deal of formalism involved that I don't understand, but from what I can tell it looks like the general idea with QSS is to piecewise create the motion for a body assuming constant velocity/slope in each piece? Specifically I'm referring to this paper. Wouldn't that mean that you couldn't properly model ballistic motion?

Also, I didn't see any mention of rotational issues. For example, something like a a block falling/rolling down stairs, which is in my mind the primary issue preventing a robust implementation. Finding the first point of collision between two bodies exactly when either body is rotating is extremely difficult, because the screw motion of a point on a (quickly) rotating body can't easily be modeled with polynomials. Or at least I don't think it can, since it involves a rotation matrix (and hence trig functions).

Plus, once you add things like spring meshes to the mix, things get ugly real fast.

Perhaps a hybrid approach is in order. Use the traditional time step approach, but use that timestep window as a bounds on the piecemeal approximation of the bodies. You lose the benefit of not checking for collisions against every body every frame, but you gain the ability to properly use numerical integration techniques. Compute the current position and the speculative next position, and the current and speculative velocities (or maybe some arbitrary point in the window [0, timestep]), and fit a spline to those control points. You then have a tightly bounded approximation for the motions of bodies, and you can build up speculative collision lists and resolve impacts at the exact time of impact. Using the splines allows easier root finding since you don't need to use a potentially error-prone numerical method.

[Enki]

Not a problem for ballistic motion at all. Like I said earlier traditional time step simulation is just the degenerate case of discreet event simulation -- they behave equivalently if the discrete event steps are all identical. So if you have a problem in one, you generally have the same problem in the other. No free lunch there. For example, in a time step sim that ballistic travel is a whole bunch of short linear approximations of equivalent time duration. The QSS stuff just says sometimes the slope changes so slowly we can take aggressive advantage of adaptive time steps, and it gives some formal proofs of how that is correct and accomplished inside a discrete event sim. The payoff is potentially doing far fewer computations, and by doing fewer computations having fewer numerical issues for the same tolerances. The price is how you organize the equations and the data, but in the end they are the same equations and the same data. That is the hybrid approach. It is just a formal adaptation of adaptive time step done one item at a time, vice all the items at the same time.

The DES does avoid forced synchronous evaluations which can create a ton of extra work, that was the point of Mirtich's Timewarp paper & implementation. But he just did the discrete parts in the collision subsystem, not in general. Kofman is doing high quality engineering production type physics on individual isolated components, not a big chaotic world. My work is trying to mesh that and more typical game/virtual environment physics together. And you are right it can be a mess, but parallel computation is getting pretty close to making that a worthwhile thing to do now.

[ewjordan]

Instead of calculating progressive conservative advancement steps and integrating the entire simulation, every body computes it's exact time of collision whenever it doesn't know its next time of collision. When the time of collision arrives, we solve the resulting island (that's the proper term, right?) and then find the next exact time of impact for each body in the island. This way a collision in one part of the simulation does not affect the bodies in another part (ie: we don't have to integrate the entire simulation when we find the TOI).

If I understand your explanation, you've outlined roughly how CA is implemented in Box2d, except that the TOI calculation is conservative instead of exact, which avoids the difficulties in getting exact TOIs. It doesn't really matter, though, because as long as you get the core shapes within the true collision envelope, the contact solver can take over, and in practice this doesn't take many iterations as long as you're not overly conservative. But the island handling is performed basically as you suggest (to my understanding, at least - Erin, please correct me if this is wrong!): TOI islands are built and solved individually rather than stepping the whole set of bodies for each time step. The island is seeded each substep by finding the minimum TOI over all contacts, and islands with TOIs greater than the full step are not touched except to read off the cached TOI and see that it is irrelevant. In fact, Box2d limits the amount of computation even further by capping the size of the TOI islands themselves, which leaves open the potential for mistakes but performs well in practice.

What do you expect to gain over CA? If it's a speed thing, I'm not sure that a spline-based approximation will necessarily save that much time, since you're resorting to an iterative search anyhow. I could be wrong, but I would think carefully about what the more standard methods lack that you really need before diving into what could be a difficult research project, especially if you're working on a real world project as opposed to an academic inquiry.

[Numsgil]

Box 2D cheats, IIRC, how it resolves collision response after CCD returns a TOI. It locks both bodies to their positions at the time of impact, and resolves the resulting collision after the entire timestep has been integrated. My main issue with this method is that it isn't "correct". Even if it's just for a video game, there's a frame where the two bodies are exactly touching, which can give a weird look to the collision. What should happen is that once the collision is detected inside the time step, it's resolved, and then both bodies are continued to be simulated for the rest of the time step. In practice this issue probably doesn't matter very often. Still, it's just not something I want to accept.

As far as CA goes, it should arrive at the same answer as other TOI algorithms, but it's iterative. If I can approximate the motion of the body with something that's easier to find roots for (like a polynomial), I can entirely avoid any iterative methods (or at least severely reduce the number of iterations needed). Hopefully meaning more speed at the cost of a little accuracy. Which is important when you need to use CCD on all bodies in your world.

[Numsgil]

@Enki:

Correct me if I'm wrong, but let's consider ballistic motion with no resistance (constant acceleration over time). If I use r(t) := r(0) + t * v(0) + a / 2 * t ^2 from introductory physics, I can compute the motion exactly. However if I use a linear piecewise approximation (which is what QSS does?), I can no longer compute the motion exactly. That is, it would seem that QSS assumes that, for each discrete event, no accelerations are taking place. Which is a rather large limitation. Unless I've misunderstood something.

The method I've been using is to construct a vector quadratic spline for the motion of each body over some time interval. You can then do something like:

(A2 - A1) t^2 + (B2 - B1) t + (C2 - C1) = at^2 + bt + c

where An, Bn, and Cn are vector coefficients for bodies 1 and 2, and a, b and c are scalar coefficients (for things like a radius term for two colliding spheres that varies over time). I then square both sides to get all the coefficients in to scalar form (using the dot product for the vector coefficients), and combine to form a quartic polynomial, which can be solved analytically for its roots (so I avoid any iteration), and gives me the exact time of collision(s), barring machine accuracy issues with the root finding method. (BTW, is there a better way to solve an equation of the form I presented?) It's then easy to disregard even roots and negative roots, and arrive at the first positive odd root. This would seem to give better results than what QSS is offering, right? Unless I'm terribly misunderstanding something.

The DES does avoid forced synchronous evaluations which can create a ton of extra work

Yes, but you only get that advantage in a specific case (the ODEs aren't stiff). Of course, there are plenty of situations where that applies, and it's a problem for traditional methods anyway...

Specifically for my uses, articulated bodies and spring meshes would be common. Once you introduce them, I don't see a way to get around doing iterative updates of their state information. Which I think, as you point out, is the same thing as a traditional time step simulation.

However, you can now compute exact times for collisions. If you update the collision information as you go, instead of waiting until the end of the frame for things to get messy, you can reduce the size of your NxN matrix for solving contact forces. You can also avoid issues with trying to guess what the contact normal should be...

[Enki]

@Enki:

Correct me if I'm wrong, but let's consider ballistic motion with no resistance (constant acceleration over time). If I use r(t) := r(0) + t * v(0) + a / 2 * t ^2 from introductory physics, I can compute the motion exactly. However if I use a linear piecewise approximation (which is what QSS does?), I can no longer compute the motion exactly. That is, it would seem that QSS assumes that, for each discrete event, no accelerations are taking place. Which is a rather large limitation. Unless I've misunderstood something.

You are simply, oversimplifying. QSS doesn't make you leave out anything, although QSS is only first order. There are QSS2 and QSS3 variations which yield even more favorable update rates, although I don't see an efficient way to optimize a collision broadphase based on a parabolic sections. It seems trivial to treat the accelerations as first order derivations of velocity and go from there that way. I think that will allow me to generate collision capsules for each object as it's event is computed. Then just compare it from the bottom up inside an octtree following insertion. Since I am only moving one thing at a time I should have a minimal set to test against as driven by the octtree and if I compute a collision I just shorten the event lifetime so that is will recompute at the time of impact and never "overshoot". So the CCD is done to schedule the collision event on the DES queue.

The method I've been using is to construct a vector quadratic spline for the motion of each body over some time interval. You can then do something like:

(A2 - A1) t^2 + (B2 - B1) t + (C2 - C1) = at^2 + bt + c

where An, Bn, and Cn are vector coefficients for bodies 1 and 2, and a, b and c are scalar coefficients (for things like a radius term for two colliding spheres that varies over time). I then square both sides to get all the coefficients in to scalar form (using the dot product for the vector coefficients), and combine to form a quartic polynomial, which can be solved analytically for its roots (so I avoid any iteration), and gives me the exact time of collision(s), barring machine accuracy issues with the root finding method. (BTW, is there a better way to solve an equation of the form I presented?) It's then easy to disregard even roots and negative roots, and arrive at the first positive odd root. This would seem to give better results than what QSS is offering, right? Unless I'm terribly misunderstanding something.

I doubt you are doing significantly better than RK4 and that can be implemented in QSS as a set of state equations for each variable and using an adaptive timestep to boot. And I really don't want to think about the comparative complexity of doing collision testing against a twisty spline as compared to simple bounding boxes or capsules.

Yes, but you only get that advantage in a specific case (the ODEs aren't stiff). Of course, there are plenty of situations where that applies, and it's a problem for traditional methods anyway...

No free lunch!

Specifically for my uses, articulated bodies and spring meshes would be common. Once you introduce them, I don't see a way to get around doing iterative updates of their state information. Which I think, as you point out, is the same thing as a traditional time step simulation.

True, but why use all those unneeded extra computations BEFORE you touch that spring mesh?

However, you can now compute exact times for collisions. If you update the collision information as you go, instead of waiting until the end of the frame for things to get messy, you can reduce the size of your NxN matrix for solving contact forces. You can also avoid issues with trying to guess what the contact normal should be...

This is a big advantage of DES like I am trying to implement it. You can NEVER have more than one collision being evaluated potentially changing another computation that was just completed since EVERYTHING is strictly done in time order. One big thing I will need to find out is whether this method ends up killing other optimizations making it less attractive, or not. That's a long ways off yet.

[Erwin Coumans]

It is better to keep TOI calculation separate from continuous collision response.

Bullet uses CA to calculate the exact time of impact. This is just a query, returning the fraction and contact points at the time of impact. To calculate the time of impact, several internal CA iterations are performed, without changing the transform of the actual bodies.
Please read the FAST paper and CATCH paper on more details on CA, time of impact calculation.

Continuous physics response is a separate topic. Time warp rigid body simulation is one way of dealing with CP, taking the discontinuities at the time-of-impact into account.

A hybrid approach for continuous physics is better in my opinion: allowing for a small amount of penetration, and entirely preventing deep penetrations using the time of impact. Such hybrid approach is implemented in the work-in-progress btContinuousCollisionWorld.

Hope this helps,
Erwin

[Erin Catto]

Box 2D cheats, IIRC, how it resolves collision response after CCD returns a TOI. It locks both bodies to their positions at the time of impact, and resolves the resulting collision after the entire timestep has been integrated. My main issue with this method is that it isn't "correct". Even if it's just for a video game, there's a frame where the two bodies are exactly touching, which can give a weird look to the collision. What should happen is that once the collision is detected inside the time step, it's resolved, and then both bodies are continued to be simulated for the rest of the time step. In practice this issue probably doesn't matter very often. Still, it's just not something I want to accept.

That was my initial implementation. The version released with 2.0 finishes the time step, so there is no visual hiccup.

The first TOI is found, then an island is built for the bodies involved. That island is solved and integrated forward to the end of the time step. Then the next TOI is found. Repeat.

[Enki]

It is better to keep TOI calculation separate from continuous collision response.

Of course. This discussion of discrete event simulation doesn't change that, but it does change how some of the TOI is done because there is a restricted subset of objects to worry about whenever the computation is done. This also gives good potential for parallelization. The response results are then used to compute and schedule the next events for the involved objects. The DES method is only generally trying to limit the number of computations needed by bypassing "unnecesary" ones, not changing to wholly different algorithms.

Bullet uses CA to calculate the exact time of impact. This is just a query, returning the fraction and contact points at the time of impact. To calculate the time of impact, several internal CA iterations are performed, without changing the transform of the actual bodies.
Please read the FAST paper and CATCH paper on more details on CA, time of impact calculation.

Continuous physics response is a separate topic. Time warp rigid body simulation is one way of dealing with CP, taking the discontinuities at the time-of-impact into account.

We are quite familiar with the Time warp method, it has been cited several times in this thread already. Not everyone needs the same boilerplate responses and provided links. Time warp actually is a DES layer within that version of Mirtich's collision system. We have been discussing the extension of that down into the rest of the engine. Matter-of-fact Mirtich got the Time Warp name and method from Jefferson's Discrete Event Simulation papers published in the early to mid '80s. Only problem for DES was few simulations partitioned nicely so Time Warp never really caught on in a big way within the DES community and they still gnash teeth and wring hands over mass rollback or eggregious synch overhead. Mirtich had a good idea and poked Time Warp a bit, but I haven't read much more since. I plan on poking it a little more strongly although without holding to Jefferson's original division of Time Warp scheduling.

A hybrid approach for continuous physics is better in my opinion: allowing for a small amount of penetration, and entirely preventing deep penetrations using the time of impact. Such hybrid approach is implemented in the work-in-progress btContinuousCollisionWorld.

Hope this helps,
Erwin

Every little bit helps. I don't want to reinvent the world, I want to make physics far more accessible within the DES paradigm and to the families of projects I am targeting. And adapting existing algorithms is far more preferable than outright invention.

[Erwin Coumans]

Indeed, those ideas are not new, I discussed them with Brian Mirtich back in 1998.

mass rollback or eggregious synch overhead

Allowing a small amount of penetration, reduces the amount of rollback and synchronization to a minimum. Contact situations usually don't cause deep penetration, so in that case the time step doesn't need to be subdivided.

Thanks,
Erwin

[Numsgil]

Box 2D cheats, IIRC, how it resolves collision response after CCD returns a TOI. It locks both bodies to their positions at the time of impact, and resolves the resulting collision after the entire timestep has been integrated. My main issue with this method is that it isn't "correct". Even if it's just for a video game, there's a frame where the two bodies are exactly touching, which can give a weird look to the collision. What should happen is that once the collision is detected inside the time step, it's resolved, and then both bodies are continued to be simulated for the rest of the time step. In practice this issue probably doesn't matter very often. Still, it's just not something I want to accept.

That was my initial implementation. The version released with 2.0 finishes the time step, so there is no visual hiccup.

The first TOI is found, then an island is built for the bodies involved. That island is solved and integrated forward to the end of the time step. Then the next TOI is found. Repeat.

Ah, that's good to know. I might poke around your code some more (I don't have a strong understanding of joints, so I was going to look at your implementation again anyhow). I do have two questions relating to Box 2D's specific implementation:

1. When you detect a mid-frame collision, do you advance/regress the entire simulation back to the time of impact? That's my biggest gripe with most continuous collision response methods (or my understanding of most continuous collision response methods, at least), since the problem has such a highly localized nature.

2. Are your results deterministic independent of the order that the bodies are resolved? That is, are you actually resolving the collisions in the order that they happen to arrive at the "correct" motion?