Number of iterations

[Dirk Gregorius]

I have heard a lot of statements in the last time about the number iterations needed to find a stable solution for stacking or articulated structures. Thanks to Erin Catto who shared a lot of information and Erwin Coumans who provides a sample implementation of this (along with his own collision library of course) for my understanding a good simulator must be able to run between 2 and 16 iterations depending on stack height or number of joints at the moment.

What is not so often mentioned is the timestep and the gravity used when measuring the iteration count. I assume that most of the time a fixed time step of 60 Hz along with a gravity of 10 is used - if not even smaller. I made the experience that the number of iterations needed to find a stable soltion drastically raises when the gravity gets higher and/or the timestep drops. So from my experience a proper simulator should run at a frequency of 30 Hz. Also higher gravity (about 30) looks much better than a gravity of 10.

I wonder what experiences others have made and how you deal with low or even worse varying timesteps? E.g. is a fixed timestep absolutely necessary or is it actually good enough to stay below some max stepsize. How do you deal with higher gravity or is a gravity of 10 simply a design principle like not stacking heavy objects on top of lighter ones. Finally I wonder how you deal with large angular velocities which can introduce a big amount of error into your system, since the linear approximation doesn't hold anymore?

Regards,

-Dirk

[kenny]

I have heard a lot of statements in the last time about the number iterations needed to find a stable solution for stacking or articulated structures. Thanks to Erin Catto who shared a lot of information and Erwin Coumans who provides a sample implementation of this (along with his own collision library of course) for my understanding a good simulator must be able to run between 2 and 16 iterations depending on stack height or number of joints at the moment.

In my opnion, if one uses a iterative projection method for solving complementarity constraints, then a competitive simulator should use around 5-10 iterations. This is easily doable for several 1000s of rigid bodies being thrown at a floor plane or similar configurations.

For jointed mechanism one would definitely need more iterations... unless one settle for tree-like structures (that is if one can break cyclic dependencies). But then again it would be much easier in my opnion to introduce specialized solvers for such articulated figures, something like the articulated body method of Featherstone would be a good candidate.

In the general case there is no telling of how many iterations that one would need. I recall having read about some people having done such experients and formulated a big O-notation for it. If my memory serve me right, the number of iterations one need with a GS-type of solver to achieve sufficient (whatever that may mean) convergence is roughly O(k d^3), where d is the diameter of corresponding graph of ones constraints ( the longest path of contacts between any two pair of rigid bodies).

The thing is that the k-coefficient depend on mass-ratios and a bunch of other stuff as well. There is no telling how a configuration would behave before one actuall try it out!

I think the point I am trying to make here is that there is no real-general purpose answer to your question:-(

What is not so often mentioned is the timestep and the gravity used when measuring the iteration count. I assume that most of the time a fixed time step of 60 Hz along with a gravity of 10 is used - if not even smaller. I made the experience that the number of iterations needed to find a stable soltion drastically raises when the gravity gets higher and/or the timestep drops. So from my experience a proper simulator should run at a frequency of 30 Hz. Also higher gravity (about 30) looks much better than a gravity of 10.

Of course the time-step size have a direct impact on the time-integration. In most cases an interleaved integration scheme is used. Creating a somewhat semi-implicit integration scheme. I believe this is a symplectic type of integrator and there is a bunch of theory out there on these things. Look up the springer book by Hairer (I forgot the book title and I am currently out of my office).

Besdies these parameters often only have a influence on the right-hand-side part of your constraint system, ie. the b-vector in

w = A lambda + b

where w and lambda is the complementary post contact velocities and the lagrange multipliers.

Exactly how the b-vector would influence the accuracy/iteration tradeoff for an iterative complementarity solver is a bit unclear to me. However, I have seen the same effects concerning gravity as you mention. Tweaking masses is although a bit more clear... This has a direct impact on the condition number of the A matrix

A = J^T W J

where

J is the constraint jacobian matrix and W the inverse mass matrix. In my opion this corresponds to increasing the diagonal parts, thus making A more diagonal dominant.

I wonder what experiences others have made and how you deal with low or even worse varying timesteps? E.g. is a fixed timestep absolutely necessary or is it actually good enough to stay below some max stepsize. How do you deal with higher gravity or is a gravity of 10 simply a design principle like not stacking heavy objects on top of lighter ones.

I try to invent better methods for time-integration and solving constraint problems:-)

Finally I wonder how you deal with large angular velocities which can introduce a big amount of error into your system, since the linear approximation doesn't hold anymore?

Yes, it is a big problem... I think Claude Lacoursiere have written about this problem. The thing is that in animation we really do not care about the gyroscopic forces most of the time, so who cares:-)

From a convergence theory viewpoint, you can just crank down the time-step size, but then we are no longer talking about real-time computer animation:-(

So, I think the question should be: Is accurate gyroscopic forces needed to create a plausible real-time game animation of rigid bodies? You'll have to ask some content-people or players about this:-)

[Dirk Gregorius]

Exactly how the b-vector would influence the accuracy/iteration tradeoff for an iterative complementarity solver is a bit unclear to me.

The b-vector contains the Baumgarte term and the relative velocity we like to eleminate. Take the ODE formulation for example:

J*W*JT * lambda = c/dt - J*(v/dt + W * f_ext)

We can write it slightly different:

J*W*JT * lambda * dt = c - J * ( v + W * f_ext * dt )

If you now look at the term v + W * f_ext * dt we see that this is simply the velocity at the next timestep without considering the constraints. So we get

J * W * JT * lambda = c - J * v(t+ dt)

This is the scheme Erin Catto introduced at the GDC this year and what Erwin is using in my opinion. Integrate forward in time without considering constraints. Apply sequentiel impulses to satisfy all constraints. Erin added some more tricks like e.g. accumulated impulses. But this is the basic idea. Also you see here the relation between sequentiel impulses and the LCP solved using a GS scheme. Note that J * v(t +dt) is the relative constraint velocity (e.g. the relatve velocity at the anchor of the ball joint) which we want to remove. Instead of correcting the relative velocity to zero we add some extra "speed" to catch numerical drift.

So what can happen here? I see basically two problems. First during a big timestep and assuming all constraints were not (simulaneously) satisfied we introduce quite a big geometric error - that are violations of the position constraints. When we catch this with the Baumgarte term we can add quite some overshoot. Projection (Cline) helps here, but only to some extent I guess. The system oscilates and needs some time to settle.

The second problem I see are big angular velocities. If you take for example the ball joint. The position constraint is:

C(x1, q1, x2, q2) = x1 + r1 - x2 - r2

We normally find the velocity constraint building the derivative w.r.t time and get:

dC/dt = v1 + w1 x r1 - v2 - w2 x r2

The suspicious part is the derivative of the lever arms to the anchor - here w1 x r1 and w2 x r2. IIRC from university the is only an approximation and not valid in general. If I find time I take out my old records and try to find it.

Tweaking masses is although a bit more clear... This has a direct impact on the condition number of the A matrix

If this is true we could easily factor the big (inverse) mass out. Given the above formulation:

J*W*JT * lambda = c/dt - J*(v/dt + W * f_ext) and lo < lambda < hi

We can factor out the biggest mass and get:

J* s * W' * JT * lambda = c/dt - J*(v/dt + W * f_ext) and lo < lambda < hi

J * W * JT * (s* lambda) = c/dt - J*(v/dt + W * f_ext) and lo < lambda < hi

Substitute s*lambda = lambda' we get

J * W * JT * lambda' = c/dt - J*(v/dt + W * f_ext) and lo < lambda' / s < hi

And finally arrive at:

J * W * JT * lambda' = c/dt - J*(v/dt + W * f_ext) and s * lo < lambda' < s * hi

The trick here is to scale the limits. I tried this and it works. Actually I needed this to deal with big mass (not big mass ratios). I didn't recognized that this improved convergence. I will look at it...

Yes, it is a big problem... I think Claude Lacoursiere have written about this problem. The thing is that in animation we really do not care about the gyroscopic forces most of the time, so who cares:-)

What do you mean by gyroscopic forces?

[Erwin Coumans]

So what can happen here? I see basically two problems. First during a big timestep and assuming all constraints were not (simulaneously) satisfied we introduce quite a big geometric error - that are violations of the position constraints. When we catch this with the Baumgarte term we can add quite some overshoot.

See Erin's comment here: http://www.continuousphysics.com/Bullet ... .php?t=273

I haven't tried this yet, but you should be able to propagate two velocity states. One velocity state includes the bias and is used to update position. The second velocity state does not include the biase and is carried over to the next step.

I'll add this to Bullet soon, should improve quality a lot.
See this comparison between PhysX and Open Dynamics Engine in XSI:
http://softimage.wiki.avid.com/index.ph ... ics_Engine

However, in just about every other way physX is superior. The "bullet points" are:

physX uses a faster physics integration algorithm, which allows the stable simulation of many, many more colliding objects than ODE. We're talking tens of thousands rigid bodies for physX, versus maybe 50 for ODE.
physX is more stable than ODE. This is related to the superior integration technique that physX uses, but the upshot on the user side is that deeply interpenetrating objects do not fly apart at a zillion miles an hour; difficult collisons (for example, an object pinched between two other objects) can be resolved cleanly, and large piles of highly constrained rigid bodies (a pile of hundreds of ragdolls, for example) will be stable and not twitch around.
physX has convex hull collision types, as well as a very stable and accurate actual-shape collision type. ODE's actual shape collision handling is "heuristic" (to put it politely) and fails when more than a few actual-shape objects are colliding with each other simultaneously.

Notice that Bullet already has proper convex support, and handles 500 object (I think ODE can handle more then 50?!?).
So fixing this 'do not fly apart at a zillion miles an hour' would bring Bullet closer to PhysX

What do you mean by gyroscopic forces?

It is the velocity dependent force, also known as Coriolis I think. Bullet doesn't add it, and in Open Dynamics Engine (version 0.6) you can switch it off with a command line (--disable-gyroscopic).

[kenny]

Exactly how the b-vector would influence the accuracy/iteration tradeoff for an iterative complementarity solver is a bit unclear to me.

The b-vector contains the Baumgarte term and the relative velocity we like to eleminate. Take the ODE formulation for example:

J*W*JT * lambda = c/dt - J*(v/dt + W * f_ext)

We can write it slightly different:

J*W*JT * lambda * dt = c - J * ( v + W * f_ext * dt )

If you now look at the term v + W * f_ext * dt we see that this is simply the velocity at the next timestep without considering the constraints. So we get

J * W * JT * lambda = c - J * v(t+ dt)

This is the scheme Erin Catto introduced at the GDC this year and what Erwin is using in my opinion. Integrate forward in time without considering constraints. Apply sequentiel impulses to satisfy all constraints. Erin added some more tricks like e.g. accumulated impulses. But this is the basic idea. Also you see here the relation between sequentiel impulses and the LCP solved using a GS scheme.

Yes, there is definitely a relation between solving a constraint-based method with a GS method and using sequential impulses. However, there is one major difference.

In the GS approach the right-hand-side b-vector remains constant, while
using sequential impulses would correspond to always updating the J^T u term in the b-vector with the best-known velocities.

One could see it as a way to improve convergence of GS even more... but I think this adds an element of randomness to the simulation....

Note that J * v(t +dt) is the relative constraint velocity (e.g. the relatve velocity at the anchor of the ball joint) which we want to remove. Instead of correcting the relative velocity to zero we add some extra "speed" to catch numerical drift.

So what can happen here? I see basically two problems. First during a big timestep and assuming all constraints were not (simulaneously) satisfied we introduce quite a big geometric error - that are violations of the position constraints. When we catch this with the Baumgarte term we can add quite some overshoot. Projection (Cline) helps here, but only to some extent I guess. The system oscilates and needs some time to settle.

The second problem I see are big angular velocities. If you take for example the ball joint. The position constraint is:

C(x1, q1, x2, q2) = x1 + r1 - x2 - r2

We normally find the velocity constraint building the derivative w.r.t time and get:

dC/dt = v1 + w1 x r1 - v2 - w2 x r2

The suspicious part is the derivative of the lever arms to the anchor - here w1 x r1 and w2 x r2. IIRC from university the is only an approximation and not valid in general. If I find time I take out my old records and try to find it.

I think you are mixing together two different topics

1) The time-stepping mehtod

2) Solving for lagrange multipliers in a constraint based method

In my opinion both of the problems you mention are related to the actual time-stepping scheme and has nothing to do with how you solve for the lagrange multipliers.

I was talking about iteration/accuracy tradeoffs in an iterative GS type solver used to solve for lagrange multipliers.

Secondly I believe both the problems you mention is rooted at the same ``core'' problem. The time-stepping scheme is derived from ``linear-approximations''

The semi-implicit time-stepping is a linear approximation to the Newton-Euler equations, kinematic constraints is basically a linear differential constraint, if you need ``higher order'' you should consider second-derivatives etc.. Also it is quite true that one can not just blindly integrate positional-constraints to derive kinematic constraints. One must ensure that the kinematic constraints are ``consistent''.

I find the discusion you have started very interesting and important, but I see it as a modelling issue. Are the assumptions we use to derive the usual time-stepping scheme sensible or correct for our purpose (real-time game physics)?

[Erin Catto]

In the GS approach the right-hand-side b-vector remains constant, while
using sequential impulses would correspond to always updating the J^T u term in the b-vector with the best-known velocities.

True, but sequential impulses also adjusts the left-hand side. Instead of

J*W*JT * lambda,

sequential impulses computes

J*W*JT * delta_lambda.

Viewed another way, at each iteration, sequential impulses computes an adjustment to lambda (delta_lambda), then applies that incremental impulse immediately to update the velocities. A necessary ingredient is tracking the accumulated impulses (lambda) and clamping delta_lambda correctly.

An aside:
In PGS, lambda is typically computed as a force. Due to the discretization of time, it is really an average force. In many situations the exact force is infinite. For a given time step, impulses and average forces are proportional. Thus, they are interchangeable.

[Dirk Gregorius]

I find the discusion you have started very interesting and important, but I see it as a modelling issue. Are the assumptions we use to derive the usual time-stepping scheme sensible or correct for our purpose (real-time game physics)?

Thanks. You brought it to the point. I want do understand where the limitations of the model are so I have solutions if I run into problems. I agree with your numbers. You should be able to run at 4 itertions for a competetive solver. If you follow Erins recommendations you can get there quickly without anoying pain. They work very well for jointed mechanism as well.

Well I agree that there is a need for specialized solvers, but I think of them if I want to combine ragdolls with the animation system (e.g. ragdolls that can get up). Here Featherstone comes into mind, but also Baraff's linear time solver. Then there are finally ropes, cloth and fluids. I haven't looked into fluids, but Erin's method gave the same results for cloth like the Jacobsen Verlet/Position method and I think it works well for ropes as well. So the question is can we solve anything with a projective solver - well I think you can do quite a lot...

One other thing that gets important when you actually got down to 4 iterations is two reduce the amount of work you do during each iteration. Again Erin's or Erwin's suggestion to apply the friction in the manifold center is a good example. It gives better results and you get down from 12 to 7 constraints for a manifold patch of 4 points, which nearly halfes the workload for each manifold in the inner relaxation loop. There are other ideas as well, but I haven't tried them out yet...

-Dirk

[Dirk Gregorius]

In many situations the exact force is infinite.

I don't understand this. Could you give an example please. Why does this happen and how does it affect the simulation?

-Dirk

[Erin Catto]

I don't understand this. Could you give an example please. Why does this happen and how does it affect the simulation?

Collisions for example. Also, David Baraff has some example of a rod sliding across a horizontal plane. The friction causes the rod to dig in, which increases the normal force, which increases the friction. To get a rigid body solution, a friction impulse is necessary.