Friction and LCP

[ngbinh]

If you want to deal with real Coloumb friction then PGS is out of question because it cannot handle non PSD matrix.

About credit of PGS, I view it like this: PGS is there for a long time, Moreau or anyone who applied PGS to handle unileteral constraints should get a credit for applying (the right) tool to this particular problem. I personally think he made a good contribution (if he's the first) to game physics.

[Dirk Gregorius]

If you want to deal with real Coloumb friction then PGS is out of question because it cannot handle non PSD matrix.

Why does the system matrix actually gets PSD in the presence of friction constraints? Is this always the case (e.g. just one contact point with friction) or does this only happen if some constraints become linear dependent? Can we maybe work around this? KenB mentioned somewhere that using something (normal impulses?) from the previous frame would help? What is this and why would this make the system matrix PD or at least better behaved? Ideas?

You mentioned to contribute / publish an implementation with a modifed Lemke algorithm? Any updates on this?

[Antonio Martini]

About credit of PGS, I view it like this: PGS is there for a long time, Moreau or anyone who applied PGS to handle unileteral constraints should get a credit for applying (the right) tool to this particular problem. I personally think he made a good contribution (if he's the first) to game physics.

i dont think anybody is trying to take credit for PGS here, the way i see it is that the Erin's work showed that SI approaches when implemented correctly and with minor but important modifications lead to a very efficient implementation of PGS. So in this respect it's a "fix" for SI approaches and takes advantage of the estabilished equivalence in order to implement PGS efficiently.

cheers,
Antonio

[ngbinh]

Why does the system matrix actually gets PSD in the presence of friction constraints? Is this always the case (e.g. just one contact point with friction) or does this only happen if some constraints become linear dependent? Can we maybe work around this? KenB mentioned somewhere that using something (normal impulses?) from the previous frame would help? What is this and why would this make the system matrix PD or at least better behaved? Ideas?

You mentioned to contribute / publish an implementation with a modifed Lemke algorithm? Any updates on this?

If you follow Coulomb friction then the big matrix is no symmetry so it's also not PSD.

You can see the formulation with contact and Coulomb friction as follow:

Stewart-Trinkle formulation with contact and Coulomb friction

mcp.PNG (27.26 KiB) Viewed 10239 times

You can see that the E and -E^T,U breaks everything here. The problem is not monotone LCP like the one you usually have but it's co-positive and (most of the time) Lemke can solve.

A better explanation of the method can be found here:

http://www.cs.rpi.edu/~nguyeb2/docs/rigid_body_lcp.pdf

About the modified Lemke, I had the theory but still can't find time to implement it. Also, John LLoyd already have a similiar optimization here:
http://www.cs.ubc.ca/spider/lloyd/fastContact.html
He even put his Java source code online.

[Erwin Coumans]

the way i see it is that the Erin's work showed that SI approaches when implemented correctly and with minor but important modifications lead to a very efficient implementation of PGS. So in this respect it's a "fix" for SI approaches and takes advantage of the estabilished equivalence in order to implement PGS efficiently.

Sequentially applying actual impulses on rigid bodies is very intuitive which makes it easy to customize it to certain needs. However it is slower compared to performing the iterations on a sparse matrix incrementally updating the scalar lambda/delta values, like quickstep does (both ODE and Bullet include PGS/SOR quickstep).

In Bullet we are moving away from the actual application of impulses towards incrementally updating scalar values that represent the applied impulse and current velocities. Several variables such as normal and relative point of application don't change during the iterations, so the cross product can be precomputed. You can interpret this as transforming the applied impulse from world space or body space into constraint axis space.

Compare those two fragments:

Code: Select all

//Optimization for the iterative solver: avoid calculating constant terms involving inertia, normal, relative position
void internalApplyImpulse(const btVector3& linearComponent, const btVector3& angularComponent,btScalar impulseMagnitude)
{
	if (m_inverseMass != btScalar(0.))
	{
		m_linearVelocity += linearComponent*impulseMagnitude;
		if (m_angularFactor)
		{
			m_angularVelocity += angularComponent*impulseMagnitude*m_angularFactor;
		}
	}
}

and the regular applyImpulse

Code: Select all

void applyCentralImpulse(const btVector3& impulse)
{
	m_linearVelocity += impulse * m_inverseMass;
}
	
void applyTorqueImpulse(const btVector3& torque)
{
		m_angularVelocity += m_invInertiaTensorWorld * torque;
}
void applyImpulse(const btVector3& impulse, const btVector3& rel_pos) 
{
	if (m_inverseMass != btScalar(0.))
	{
		applyCentralImpulse(impulse);
		if (m_angularFactor)
		{
			applyTorqueImpulse(rel_pos.cross(impulse)*m_angularFactor);
		}
	}
}

If you optimize this internalApplyImpulse a bit further, you avoid updating the linear and angular velocity vectors, and just keep track of the scalar incremental values (the angularComponent and linearComponent don't change). Then the sequential impulse based constraint solver converges towards the PGS/SOR quickstep constraint solver.
Hope this helps,
Erwin

[Dirk Gregorius]

However it is slower compared to performing the iterations on a sparse matrix incrementally updating the scalar lambda/delta values, like quickstep does (both ODE and Bullet include PGS/SOR quickstep)

I think this should read that Bullet's SI implementation is slower than ODE's Quickstep. Even the ODE Quickstep can be made much more efficient (e.g. the construction of the Jacobians and the right hand side of the linear system). Also the SI can be implemented in a more efficient way (making it more cache friendly as you already mentioned several times on the forum before). I don't necessarily see that SI is less efficient than PGS. Expecially when using NGS for position correction where the Jacobians change after *each* applied position impulse (push).

You can interpret this as transforming the applied impulse from world space or body space into constraint axis space.

I see that this works for 1D constraints, since you can use the following identities:
v' = v + P/m = v + lambda * n / m = v + n / m * lambda
w' = w + L / I = w + ( r x P ) / I = w + ( r x lambda * n ) / I = w + ( r x n ) / I * lambda

Here n / m and ( r x n ) / I can indeed be precomputed and the expensive matrix vector multiplication becomes an add-scale operation. How do you extend this to 3D constraints, e.g. back-socket joint which you usually want to solve as a block and not as three individual constraints?

[Erwin Coumans]

I think this should read that Bullet's SI implementation is slower than ODE's Quickstep.

quickstep is now native part of both Bullet and ODE. So it should read: any SI implementation that applies actual impulses modifying the linear and angular velocity of rigid bodies is slower then Russell Smith implementation of quickstep, available in Bullet under the Zlib license and in ODE under the combined BSD/GPL license. According to Vedran Klanac, Russell got some inspiration from Crotream:

Background history, how quickstep method was developed is totaly different. The exact algorithm name is Gauss-Southwell, and it is special case of SOR and Gauss-Seidel inherently.

The main idea for the algorithm evolved from improved factorization method for direct Danzig solver that I developed also earlier. At specific point in my research I figured out that iterative scheme must be more effective and less time consuming for CPU per frame. Early versions of algorithm which would solve separately all equations in stacked matrix rather then calculating sparse matrix blocks proved to be the right way to go. If you go step by step and multiply just one row of constraint equation matrices, jacobians for source and destination body with inversed mass matrix, the one can see that A matrix becomes a scalar. And all you have to do is go through the stacked matrix and perform bunch of divisions and clamp result value inside lo and hi limits. It is simple as that. I did it in first place on the piece of paper in a form of proof.
I must say that several e-mails from Mr. Richard Tonge on ODE mailing list were very resourceful to guide me towards iterative scheme.

Everything fell in its place when I read Murtay's book on iterative LCP solvers. First I implemented SOR and later on it turned out that omega factor for SOR which fit best solution convergence was 1, and that is exactley the value which makes SOR being Gauss-Seidel method. Gauss-Southwell method is different in a way that you don't calculate new solution in every frame from the beginning, but you actually maintain error of the solution and thus save some CPU time.

That's how quickstep method came to the day of light. The whole thing happened in april 2004, few months later we contributed method to ODE.

I agree with you that both PGS and SI can be optimized in identical ways. It is just that SI encourages one to apply impulses on 3D vectors. Once you avoid applying impulses I'm not sure if it should still be called Sequential Impulse.

How do you extend this to 3D constraints, e.g. back-socket joint which you usually want to solve as a block and not as three individual constraints?

Interesting, how would you perform a 'sequential impulse' on a 3D constraint? Isn't that just 3 sequential applications of a 1D constraint in SI? Once you introduce other methods, such as NGS and block-wise solving it is better to stop applying impulses within the block, and just solve a small linear system. This is probably what the Subspace minimization is about.

Thanks,
Erwin

[Dirk Gregorius]

So it should read: any SI implementation that applies actual impulses modifying the linear and angular velocity of rigid bodies is slower then Russell Smith implementation of quickstep

It depends on the implementation. If you implement it in a fashion like below I can imagine that this is quite slow:

for (all joints j )
j->Satisfy()

Still there is no need to do this. You can implement SI in the same fashion as PGS - just think of two accumulators (lambda and v). Also note that for low iterations (say 4 - 8) the constraint preparation (e.g. building Jacobians) can be slower than the actual constraint solving. Here the ODE is quite inefficient in my opinion.

What about Tokamak? If seen representative tests where it was the second fastest engine behind Havok and faster than ODE. Doesn't Tokamak use SI as well?

Interesting, how would you perform a 'sequential impulse' on a 3D constraint?

You only compute the constraint strength (lambda). As you said correctly for 1D constraints the direction is given and constant. You use this and "bake" the direction into the impulse equation which saves you the matrix multiplication for the angular impulse. More generally you solve:

J*W*JT * lambda = -J *v and P = JT * lambda

This can be of any dimension. For a contact constraint this is just one equation for one variable and P = JT * lambda is actually just the scaling of lambda with the contact normal, for a ball socket constraint this becomes a linear system of three unknowns and three equations. You could solve this as three individual constraints (the ODE does this). You can also solve the 3D linear system and apply the 3D lambda vector at the pivot point which is equivalent to building the four impulse vectors (P = JT * lambda). There is no direction anymore that could be baked into the impulse equation, and you get a matrix multiplication for the angular impulse again. The way around this is to perform some aspects of the solving process in local coordinates which makes even the Jacobians constant over the entire simulation. This is quite ugly math though and requires some assumptions about the inertia tensors. Note that the later is not equivalent to the described method above.

[Erwin Coumans]

What about Tokamak? If seen representative tests where it was the second fastest engine behind Havok and faster than ODE. Doesn't Tokamak use SI as well?

Bullet should be faster then Tokamak and ODE on Playstation 3 and other consoles, due to parallel optimizations. Expect Bullet to be faster on PC in one of the upcoming contribution, and similar performance to Havok on PlayStation 3. When switching Bullet to use quickstep and box-box, there shouldn't be a difference with ODE and concave triangle meshes in Bullet should be faster then ODE.

There you got it
Thanks,
Erwin