Wednesday, October 26, 2005

Constraint Force Mixing

Today I learned the magic of CFM (constraint force mixing). If CFM is not introduced, the matrix A will somehow manages to get itself into singularity at some uncertain timesteps. I think physically the system just moves to a lock-up position when A becomes singular. The simplest solution used by ODE is to add a (diminutive) positive value (CFM) to the diagonal of A such that A will keep away from singularity. A quote from the ODE manual - A nonzero value of CFM allows the original constraint equation to be violated by an amount proportional to CFM times the restoring force lambda that is needed to enforce the constraint.

However, since the global CFM is an arbitrary value added to the diagonal of A, it is inherently unfair comparison when two different As are formulated and applied with the same value of CFM. This case happens when I try to compare the results of the Smith and Baraff formulations. Without the errors introduced by the CFM, the results match pretty well while subjected to some numerical roundoff errors. But the risk is either result may blow up under different circumstances without the help of the CFM. Adding the CFM allows the errors to accumulate and my guess is the errors will become substantial after many timesteps.

The system I used for this method comparison is of course my favorite piston assembly simulation. I may just try imposing the CFM to only one of the boundary joints to reduce the errors while retaining stability. After looking around for a while, it seems that only the CFM of the limit motor can be changed ...

Linear-Time Dynamics using Lagrange Multipliers

I finally implemented Baraff's algorithm from the paper, "Linear-Time Dynamics using Lagrange Multipliers". I have compared the results of the Baraff method with those of the original method (Smith's ODE) using the piston assembly simulation. The holy grail seems to be sparse LCP, but now I am content with the current implementation. According to the paper, the idea is to first solve for the primary constraints using the sparse LDLT decomposition, and then solve for the auxiliary constraints using the LCP method. The followings are the essential steps.

Given the external force F and the velocity vt at the current time step, the goal is to find the velocity at the next time step, vt + h, where h is the time step size. The system equation can be formulated as


where M is the system inertia matrix, J1 is the primary constraint Jacobian, J2 is the auxiliary constraint Jacobian, λ is the primary constraint impulse, and μ is the auxiliary constraint impulse.

First, solve


using sparse LDLT decomposition to yield

Next, formulate A as

A = J2 M - 1 K

where K is defined as

where B is solved from the following equation

Note that the LDLT factorization has been done for the matrix

in the previous step, so solving for B is rather efficient.

Finally, formulate and solve the following equation for μ using the LCP method.

The velocity vt + h is given by