This part is hairy, and the reason why I wrote this tutorial/summary in the first place. Much of the scientific literature explains the mathematics in a rather abstract way. After finally understanding both the theoretical and the practical side, I believe it's better to gather an intuitive, less founded understanding first, then proceed to the theory behind it such as LCPs.
Constraints are needed for bodies to interact - to avoid penetration when bodies rest on each other, or collide; or when bodies are connected at specific points (joints as for doors, ragdolls, ropes, ...). I found it very instructive to understand different classes of constraint solvers - penalty methods in force, velocity or position space; or direct/simultaneous solvers. Garstenauer gives a really good overview over these. I will present one specific method - the most widespread one - here, then later on describe what class it belongs to and what other classes there exist. For now, it's enough to know it's a sequential impulse-based constraint solver.
Constraints are formulated as constraint equations (or inequations). For instance, to model a joint between two bodies connected at points of body A and of body B, one starts with the constraint equation and tries to maintain a value of , i.e. to keep the two points at the same location. A distance constraint would be modeled with .
Impulse and Force Based Constraints
One way to ensure the constraint is fulfilled is to start at a valid configuration of positions (i.e. ) and ensure it stays like this by aiming for (i.e. aiming for a valid configuration in the velocity domain); this is called a velocity or impulse based constraint solver. An alternative is an acceleration or force based constraint solver - here, the constraint equation is differentiated twice, and one aims for . I'll focus on the more common velocity based solvers. The different types of constraint solvers used in contemporary physics engines - impulses vs forces, sequential vs simultaneous, were discussed in this insightful thread on the Bullet forum.
In velocity based constraint solvers, one searches for a valid configuration by directly changing the (linear and angular) velocities of the involved bodies. The velocities are changed by applying constraint impulses. Impulses - denoted by - lead to a discontinuous change in velocities, weighed by the mass or moment of inertia of the body. It is easy to convert between forces and impulses .
Impulses change velocities; they are applied in a similar manner to forces. There, we had
Velocities are changed by impulses via , with
Like forces, an impulse really is a 2D vector that is expanded to a 3D vector to simplify the formulas.
This might look a bit complex at first, but the actual implication is rather straightforward: Physics engines start at an invalid configuration (i.e., where ), then perform an iteration over all constraints, for each computing the right impulse , changing the velocities via . This process is repeated about 4 times, and hopefully, the errors have become small, and the corrected velocities are used to compute the new positions and orientations.
Note for native german (also: danish, dutch & russian) speakers: impulse is a false friend; impulse (eng.) != Impuls (dt.); impulse (eng.) == Kraftstoss (dt.), wheres momentum (eng.) == Impuls (dt.)
Computing Constraint Impulses
This section shows how, given a constraint equation , we can derive the required impulse to apply to the velocities; it is heavily based on Erin Catto's GDC2009 talk. Concerning notation - thus far, the impulse and velocity vectors had 3 elements, and the mass matrix was 3x3. Depending on whether the constraint involves one or two bodies, the vectors in the following formulation may have 3 or 6 elements, by combining the impulses or velocities for both bodies in one vector.
The temporal derivative of a constraint equation , according to the chain rule, is ; is the so-called Jacobian, a row vector. is a state vector representing the position and orientation of the involved bodies (usually 2). In practice, one adds a bias term, to require ; this term counteracts numerical errors / drifting, more about it at the end of this section. Part of the bias term may also arise from the constraint equation directly (for motors and collision constraints).
To compute the constraint impulse , we need three ingredients:
1. To ensure the constraint,
2. Constraint forces or impulses may not add energy to the system, i.e. do work - force times distance: . Since , the impulse we search must be a multiple of the transposed Jacobian:
NB, some tutorials use the bead-on-a-wire example and state that the constraint force or impulse stands orthogonally on the velocity. This may lead to a misconception; I wouldn't attach much geometrical meaning to , since both are accumulated vectors. For a joint constraint, the individual impulses roughly point from one point to the other.
3. As discussed further above, the impulse is applied to the previously predicted, tentative velocity via
Replacing the impulse in equation 3 with , then replacing the velocity in equation 1 gives .
Solving for lambda gives the magnitude of the constraint impulse:
Then, a new velocity is predicted by
The above computations are performed in two nested loops; the outer one performing a small, usually fixed set of e.g. 4 iterations; the inner loop iterating over all constraints, for each computing the new and updating the velocities of the involved bodies. Since only velocities change in this loop, the Jacobian can usually be precomputed before starting the iterations.
Some notes on the bias term :
using the constraint equations as described above can lead to position drift, this can be avoided by a method called Baumgarte stabilization, where the original constraint equation is fed back into the velocity constraint , by using a modified , with (needs to be determined manually).
I don't have an intuitive explanation for this method, but mathematically, the term is a first order linear differential equation, with a solution at - a function that decays exponentially with time to 0, i.e. a valid configuration
In practice, one computes the original constraint , sets and proceeds with the equations above. I use and it works well. Why the constant is divided by I don't know.
Example: Complete Joint Constraint Derivation
This section infers the precise equations for joint constraints used for the rope simulation at the end of this chapter. I think it makes a lot of sense to walk through the complete derivation and discuss some details on the way; that being said, the mathematics here are lengthy, but straightforward linear algebra.
The final formulation above does not look complex, but computing the formula of the Jacobian directly (by differentiating with respect to positions and orientations of both bodies) may actually be rather gory. I find it easier, and Erin Catto recommends this, too, to compute the temporal derivative directly, then rearrange the formula to isolate and .
The derivation below is not difficult really but lots of writing; I would recommend to go through it or compute it yourself to get the necessary intuition.
Let's take the joint constraint equation: :
: the point on body A we want to apply the constraint to, and its linear velocity
: the center of mass / rotation of body A, and its linear velocity
: the angular velocity of body A (in radians/s)
: a pseudo 2D cross product used in computing a point's velocity: . simply computes the normal vector of in counter-clockwise direction, and scales it with
and the same notation for the point on body B of course
this ugly formula needs to be split up into the x and y components to make it even uglier:
then, the dot product is computed and all the velocity components are put into one vector, giving rise to the following form:
That was a bit lengthy but important. One can now simply plug in the above formula for into the generic constraint solver discussed above, in addition to the bias parameter . This is the general way to go when adding new constraints: one only needs to calculate the Jacobian and the bias factor, then it is straightforward to plug into a constraint solver! Note that in spite of the complex formulas, the computation should be pretty efficient by reusing computed values and using vector instructions.
I noted that the constraint impulse has an intrinsic dimensionality of 2 (or four for two bodies) only - - but is expanded to a vector with 3 elements to incorporate the angular momentum, giving . As you can see, this intrinsic dimensionality is fulfilled by the Jacobian we computed above: Since , we have , . If you arrive at a Jacobian that does not have such a form, you likely made a mistake in your calculation.
As mentioned during the overvie, modern physic engines seem to use mostly iterative solvers that work as follows (pseudocode):
for i = 1 to nIterations
for c in constraints
P = computeImpulses(c);
The set of all constraints are usually coupled: a body that shares two constraints with two other bodies, such as a link in the middle of a chain, is affected (and affects) both connected links. However, this is neglected when solving for pairwise constraints individually, and satisfying one constraint will usually violate the other constraint(s).
This type of solving the constraints - similar to the iterative Gauss-Seidel method for linear systems of equations - is just one way, but the most widespread and efficient one apparently.
Alternatively, one could attempt to solve the collisions directly, a process that can be sped up by isolating uncoupled submatrices etc., but it seems like the iterative methods are good enough and most efficient.
Note, I did not add parameters to this simulation;
frankly, because using links with less mass leads to jitter and I didn't bother taking care of it.
I simply went to 1/4th of the previous time step size - 4ms.
I suppose that rope simulations with a larger number of links always lead to problems; if one uses
not sequential impulses but springs to simulate it, one gets a stiff system of equations.
There's probably a way to take better care of it though; but keep in mind a rope is a rather extreme
application of joints because of the many couplings. Other, similarly complex applications of
joint chains like ragdolls are treated specifically as far as I know.