Position Correction

Post Reply
Erin Catto
Posts: 316
Joined: Fri Jul 01, 2005 5:29 am
Location: Irvine
Contact:

Position Correction

Post by Erin Catto »

I tried the several algorithms for position correction of Box2D’s revolute joint.

I looked at these systems:

- simple pendulum (1m diameter sphere on a massless 5m stick) with an initial angular velocity of 100 rad/s.
- suspension bridge with 30 1m long planks.
- multi-link chain with 30 1m long links.

Here are the algorithms:

Baumgarte - A fraction of the position error is added to the velocity error. There is no separate position solver.

Pseudo Velocities - After the velocity solver and position integration, the position error, Jacobian, and effective mass are recomputed. Then the velocity constraints are solved with pseudo velocities and a fraction of the position error is added to the pseudo velocity error. The pseudo velocities are initialized to zero and there is no warm-starting. After the position solver, the pseudo velocities are added to the positions.

This is also called the First Order World method or the Position LCP method. It is also very similar to the obsolete split impulses technique.

Modified Nonlinear Gauss-Seidel (NGS) - Like Pseudo Velocities except the position error is re-computed for each constraint and the positions are updated after the constraint is solved. The radius vectors (aka Jacobians) are re-computed too (otherwise the algorithm has horrible instability). The pseudo velocity states are not needed because they are effectively zero at the beginning of each iteration. Since we have the current position error, we allow the iterations to terminate early if the error becomes smaller than b2_linearSlop.

Full NGS or just NGS - Like Modified NGS except the effective mass is re-computed each time a constraint is solved.

I created four demos to test the different solvers. You can download them here: http://www.gphysics.com/files/PositionC ... rithms.zip

Observations:

Baumgarte - this is the cheapest algorithm but it has some stability problems, especially with the bridge. The chain links separate easily close to the root and they jitter as they struggle to pull together. This is one of the most common methods in the field. The big drawback is that the position correction artificially affects the momentum, thus leading to instabilities and false bounce. I used a bias factor of 0.2. A larger bias factor makes the bridge less stable, a smaller factor makes joints and contacts more spongy.

Pseudo Velocities - this is more stable than the Baumgarte method. The bridge is stable. However, joints still separate with large angular velocities. Drag the simple pendulum in a circle quickly and the joint will separate. The chain separates easily and does not recover. I used a bias factor of 0.2. A larger value lead to the bridge collapsing when a heavy cube drops on it.

Modified NGS - this algorithm is better in some ways than Baumgarte and Pseudo Velocities, but in other ways it is worse. The bridge and chain are much more stable, but the simple pendulum goes unstable at high angular velocities. The bias factor is 1, so no tuning is necessary.

Full NGS - stable in all tests. The joints display good stiffness. The bridge still sags, but this is better than infinite forces. The bias factor is 1, so no tuning is necessary.

Recommendations:

Pseudo Velocities are not really worthwhile because the bridge and chain cannot recover from joint separation. In other cases the benefit over Baumgarte is small. Modified NGS is not a robust method for the revolute joint due to the violent instability seen in the simple pendulum. Perhaps it is viable with other constraint types, especially scalar constraints where the effective mass is a scalar.

This leaves Baumgarte and Full NGS. Baumgarte has small, but manageable instabilities and is very fast. I don’t think we can escape Baumgarte, especially in highly demanding cases where high constraint fidelity is not needed. Full NGS is robust and easy on the eyes. I recommend this as an option for higher fidelity simulation and certainly for suspension bridges and long chains.

Full NGS might be a good choice for ragdolls, especially motorized ragdolls where joint separation can be problematic. The number of NGS iterations can be reduced for better performance without harming robustness much.

NGS is tricky to implement for contacts because re-evaluated contact points can be expensive, however it is possible to estimate the new penetration depth based on the body movement and the contact normal. Another possibility would be to cache features and quickly re-do the feature clipping to generate new points. Some points might disappear during iterations, so there are certainly some edge cases to consider.

Each joint in a can be handled differently in the position solver. So I recommend a system where the user can select the algorithm on a per joint basis. I would probably default to the slower Full NGS and let the user select the faster Baumgarte method in performance critical scenarios.
Dirk Gregorius
Posts: 861
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA

Re: Position Correction

Post by Dirk Gregorius »

I made the same observation and currently support full NGS and Baumgarte. The reasons is that Baumgarte is by far sthe fastest method (especially for simple destruction) and full NGS gives you nice stability and quality. Note that this can still become unstable, but you need sophisticated mechanisms (e.g. mobilees) and high angular velocities. The other methods (especially that crappy split impulse) gives very few (if any) imrovements at the price of some noticeable impact on the performance. Even worse the split impulse seem to effect the convergence and the system does not converge in all cases anymore. Most of the instabilities are results of the non-linearity of the position constraints (what becomes very obvious with stiff sub-solvers) and the full NGS gives you controll over this. Mostly you get away with more iterations for the position correction (what might be cheaper than smaller substeps). It was an nice observation that even with perfect satisfied velocity constraints you still need several iterations for the position constraints.

You mention that you want to mix both methods. I like the idea, especially since I don't want to loose the stiff spring (soft distance constraints) that you get with Baumgarte. On the other side I doubt that this really will work (e.g. mixing both methods in a ragdoll). Has anybody tried this?
User avatar
Erwin Coumans
Site Admin
Posts: 4221
Joined: Sun Jun 26, 2005 6:43 pm
Location: California, USA
Contact:

Re: Position Correction

Post by Erwin Coumans »

Thanks for sharing the results of your experiments.

I got 3 questions, applying this discussion to contact constraints:

1) Did you notice problems for the 'split impulse' approach for contacts constraints (apart from performance)? It seems to solve the problem that Baumgarte suffers from, adding unwanted velocity due to position correction.

2) How does NGS compare with subdividing the timestep, and running the entire simulation loop at higher frequency at one solver iteration with warmstarting.

3) Can you compare your NGS approach to Gino van den Bergen's dynamics loop, for contacts? See description below.

It would be interesting to know the similarities and fundamental differences for both (2) and (3).
Thanks,
Erwin




Here is a brief description of Gino's method, for your reference, below are some snippets from in solid-3.5.6\examples\physics.cpp, see moveAndDisplay():

------------------------
stepSimulation:

[...]

for number of velocity iterations
solve pairwise velocity error (impulse based)

for number of positional iterations
solve pairwise positional (penetration) error (without taking mass/inertia into account, called relax)
update position, and perform collision detection
------------------------

Code: Select all

void moveAndDisplay() 
{
   int i;
   for (i = 0; i < NumObjects; ++i) 
   {
	   object[i].applyForces(TimeStep);
	   object[i].proceed(TimeStep);
   }
   
   int k = MaxImpulseIterations;
   while (k-- && DT_Test(g_scene, g_respTable))
   {
	   for (i = 0; i < NumObjects; ++i) 
	   {
		   object[i].backup();
		   object[i].proceed(TimeStep);
	   }
   }
   
   DT_Test(g_scene, g_respTable);
   
   k = MaxRelaxIterations;
   while (k-- && DT_Test(g_scene, g_fixRespTable))
   {
	   for (i = 0; i < NumObjects; ++i) 
	   {
		   object[i].relax();
	   }
   }
   
   display();
}
Notice that the backup/proceed combination inside the velocity correction iteration loop is a way to 'sense' all contact points. It is updating the position using the latest velocity starting from the initial backed-up transform.

Within this velocity correction iteration loop, the combinatio of 'proceed' and 'DT_Test' is performing velocity correction, by using collision detection which calls the response callback, which is a pair-wise 'body2body' resolveCollision:

Code: Select all

void resolveCollision(RigidBody& body1, const MT_Vector3& pos1,
                      RigidBody& body2, const MT_Vector3& pos2,
                      MT_Scalar depth, const MT_Vector3& normal, MT_Scalar restitution)
{
	MT_Vector3 rel_pos1 = pos1 - body1.getPosition(); 
	MT_Vector3 rel_pos2 = pos2 - body2.getPosition();
	
	MT_Scalar rel_vel = normal.dot(body1.getVelocity(rel_pos1) - body2.getVelocity(rel_pos2));
	if (rel_vel < -MT_EPSILON) 
	{
		MT_Vector3 temp1 = body1.getInvInertiaTensor() * rel_pos1.cross(normal); 
		MT_Vector3 temp2 = body2.getInvInertiaTensor() * rel_pos2.cross(normal); 
                          MT_Scalar rest = restitutionCurve(rel_vel, restitution);
		MT_Scalar impulse = -(1.0f + rest) * rel_vel / 
			(body1.getInvMass() + body2.getInvMass() + normal.dot(temp1.cross(rel_pos1) + temp2.cross(rel_pos2)));
		
		MT_Vector3 impulse_vector = normal * impulse;
		body1.applyImpulse(impulse_vector, rel_pos1);
		body2.applyImpulse(-impulse_vector, rel_pos2);
	}
}
'relax' is performing pairwise positional correction, without taking mass/inertia into account.

Thanks,
Erwin
melax
Posts: 42
Joined: Mon Apr 24, 2006 4:55 pm

Re: Position Correction

Post by melax »

minor point:
It appears the Solid library's physics example just translates rigid bodies to deal with penetration. i.e. no rotation. There could be situations where some rotation is required to remove penetration. Imagine a dynamic box between two walls where the box somehow got rotated and is penetrating the walls on both sides. See attached picture. If the walls and box aren't spinning/moving, there is nothing to encourage the velocity part of the solver to reorient the box. In contrast any of the solutions that work at the contact points (fire an impulse along the contact normal or whatever) will correct this case.
brick wedged between two walls requires rotation to avoid penetration.
brick wedged between two walls requires rotation to avoid penetration.
wedged_brick2.jpg (13.32 KiB) Viewed 33852 times
Erin Catto
Posts: 316
Joined: Fri Jul 01, 2005 5:29 am
Location: Irvine
Contact:

Re: Position Correction

Post by Erin Catto »

Erwin,

You raise a good point. My discussion has mainly centered on joints, not on contacts.

1) The split impulse method (aka Pseudo Velocities, aka First Order World) seems to work well for contacts. But the method is worse than Baumgarte for joints. NGS is my favorite method for joints (so far), but I don't have a favorite method for contacts. That is an issue I'll leave for another day.

2) With the demos you can do this. I set the period to 120Hz. The Baumgarte and Pseudo Velocity solvers were no better or worse on the Bridge. You should try this yourself using the GUI.

3) Gino's position correction algorithm is a stripped down NGS. The algorithm doesn't consider mass, rotation, or Lagrange multipliers. The algorithm only projects translations.

Gino's velocity solver is a stripped down sequential impulse algorithm. It doesn't use warm starting or accumulate impulses for clamping. The basic structure seems okay (not sure about the backup part).

What is your motivation for this comparison Erwin? Did you notice some interesting qualities in Gino's solver? (I haven't run the demos).
Antonio Martini
Posts: 126
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London

Re: Position Correction

Post by Antonio Martini »

Erin Catto wrote:I tried the several algorithms for position correction of Box2D’s revolute joint.
...
interesting comparisons, if i read it correctly they are also consistent with what intuition suggests, the more up to date the used information is(effective mass, positions, jacobians etc...) the better the solution behaves.

cheers,
Antonio
Antonio Martini
Posts: 126
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London

Re: Position Correction

Post by Antonio Martini »

AntonioMartini wrote:
Erin Catto wrote:I tried the several algorithms for position correction of Box2D’s revolute joint.
...
interesting comparisons, if i read it correctly they are also consistent with what intuition suggests, the more up to date the used information is(effective mass, positions, jacobians etc...) the better the solution behaves.

cheers,
Antonio
the idea of using the most up to date information is as we know also the difference between a Jacobi iteration and a GS iteration. A jacobi solver uses only information from the previous iteration when solving for a new variable.
Now let's say we start with a square cloth having many constraints and at each iteration we pick and solve constraints at random(like it has been suggested in the past) accordingly to a uniform distribuition. Given that the distribuition is uniform it's very unlikely that initially we will pick nearby constraints and so when we solve for the first n constraints we will be very likely use only values from the previous iteration. So as a matter of fact it would look more like a Jacobi style solution in the initial stage. This fact alone would seem to suggest that picking constraints at random isn't a good idea after all.

a constraint joining bodies B0 and B1 for can be written C(0,1)

for simplicity let's build a chain made of 5 bodies:

B0,B1,B2,B3,B4,B5

and constraints:

C(0,1), C(1,2), C(2,3), C(3,4), C(4,5)

if we solved the above constraints in the given order it would be a GS style solution, we also know that a Jacobi solver uses only the information from the last iteration, so if we change the order of solving as follow for example:

C(0,1), C(2,3), C(4,5), C(2,3), C(1,2), C(3,4)

we see that the first 3 constraints use no "fresh" information, if we just rename the indices of the constraints we can see that 3 constraints are Jacobi style steps with 3 GS style steps. Luckily chains are usually built serially so it's very unlilkely that we have that order. However for Ragdolls, clothes, tetrahedral meshes it maybe very different. so what's the optimal sequence? in general a body can be connected by n contraints to other bodies(or particles), let's say we add a counter n_i for each body i, everytime we solve a contraint connected body i we increase n_i , so a mesaure of how up to date the information for that body is can be defined by:

f(i) = n_i/max(n_i) ( max(n_i) = maximum number of connected constraints to body i)

so when the body is fully updated f(i) = 1 when it's using only "old" information(and so Jacobi style) is equal to zero.

so when choosing the next optimal constraint to solve for we should aim at picking the constraint C(a,b) with maximum g(a,b) = f(a)*f(b). Note that f(i) can be interpreted as a linear fuzzy membership function and "*" can be seen as a fuzzy AND operator.see also:

http://www.atp.ruhr-uni-bochum.de/rt1/s ... de118.html
http://www.atp.ruhr-uni-bochum.de/rt1/s ... de117.html

so we pick and solve the constraint that connects the 2 most up to date bodies(or particles)
However this choice could be only locally optimal, a global solution could be the one of finding a constraints ordering that minimise the sum of all g(i,j) for a single iteration (there is one g(i,j) for each contraint). for a fixed connectivity(ragdoll, clothes, volumes) this could be done offline.

i haven't got a chance to try this method yet so take it more as an idea rather than a solution, however i would be curious to see what kind of difference it would make, especially for highly connected systems with a fixed topology like clothes or tetrahedral meshes...

cheers,
Antonio
Erin Catto
Posts: 316
Joined: Fri Jul 01, 2005 5:29 am
Location: Irvine
Contact:

Re: Position Correction

Post by Erin Catto »

When I was studying constraint ordering, I would simulate a stack of boxes or a chain and export the MLCP to a text file and load it into Matlab. It was quite easy to try different orderings and measure convergence. You may get better results, but the best I could do was shave off one iteration.

In my experience the performance killer for constraints is cache misses. Random access to constraints can be eliminated, giving a huge performance gain. As far as I can tell, random access to body state data (position, velocity) can not be eliminated.
Antonio Martini
Posts: 126
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London

Re: Position Correction

Post by Antonio Martini »

Erin Catto wrote:When I was studying constraint ordering, I would simulate a stack of boxes or a chain and export the MLCP to a text file and load it into Matlab. It was quite easy to try different orderings and measure convergence. You may get better results, but the best I could do was shave off one iteration.
i was thinking more about clothes and tetrahedral meshes, for rigid bodies i dont see many large systems with a fixed topology, usually chains and bridges are already arranged in a optimal way by construction i believe and as they are serially connected there are not many combinations available to be explored, just the bad ones:) then also for clothes and solids it could well be that there is not that much to be gained, however saving a couple of iterations by just rearraging the data is already a good result.

cheers,
Antonio
Erin Catto
Posts: 316
Joined: Fri Jul 01, 2005 5:29 am
Location: Irvine
Contact:

Re: Position Correction

Post by Erin Catto »

Ahh. Yeah, for cloth constraint ordering can make a big difference, especially if you use a stability factor as in Andrew Megg's cloth presentation.

https://www.cmpevents.com/GD05/a.asp?op ... essID=3931
mewert
Posts: 52
Joined: Sat Oct 08, 2005 1:16 am
Location: Itinerant

Re: Position Correction

Post by mewert »

I found that using a stability factor caused some bad artifacts when simulating non-trivial cloth. So I left it out. For console performance we were getting killed by load-hit-stores. This is caused by reading from a vertex that has just been changed ( still on it's way to being stored in the cache ). Which is what accessing your mesh serially will do. So we ended up trying to randomize our vertex access as much as possible. How quickly optimization strategies change with new hardware!
User avatar
Dennis
Posts: 12
Joined: Sat Jun 23, 2007 9:13 pm
Location: San Francisco
Contact:

Re: Position Correction

Post by Dennis »

Hi all, I also experimented with constraint ordering for rigid body problems some time ago, but I never found a method that was noticably better than random ordering. Multiple iterations tend to smoothen out the difference a lot. Consider four iterations. The difference (in terms of dereferenced body velocities) is not huge:

C(0,1) C(2,3) C(1,2) C(0,1) C(2,3) C(1,2) C(0,1) C(2,3) C(1,2) C(0,1) C(2,3) C(1,2)
C(0,1) C(1,2) C(2,3) C(0,1) C(1,2) C(2,3) C(0,1) C(1,2) C(2,3) C(0,1) C(1,2) C(2,3)

Cheers,
Dennis
Antonio Martini
Posts: 126
Joined: Wed Jul 27, 2005 10:28 am
Location: SCEE London

Re: Position Correction

Post by Antonio Martini »

Dennis wrote:Hi all, I also experimented with constraint ordering for rigid body problems some time ago, but I never found a method that was noticably better than random ordering. Multiple iterations tend to smoothen out the difference a lot. Consider four iterations. The difference (in terms of dereferenced body velocities) is not huge:

C(0,1) C(2,3) C(1,2) C(0,1) C(2,3) C(1,2) C(0,1) C(2,3) C(1,2) C(0,1) C(2,3) C(1,2)
C(0,1) C(1,2) C(2,3) C(0,1) C(1,2) C(2,3) C(0,1) C(1,2) C(2,3) C(0,1) C(1,2) C(2,3)

Cheers,
Dennis
i believe that if there is an optimal ordering(why not?) surely real-time rigid bodies would not benefit from it for the the following reason:

the only large systems are typically stacks of objects with a highly varying connectivity. So the optimal sequence would not be computable in real-time.

then you have also multiple iterations for Jacobi but we know that GS is better for most cases. What you gain/loose in one iteration is propagated to multiple iterations i suppose.

so clothes and tetrahedral meshes with a fixed topology would be the field of application. Here we have many constraints in the order of thousands and so the difference maybe noticeable especially if the optimal ordering is compatible with the architecture we are implementing it on.

cheers,
Antonio
Post Reply