This post is part of my Game Physics Series.

### Constraints

The resolution phase of a constraints-based physics engine uses the concept of __constraints__. A free rigid body in 3D has 6 __degrees of freedom__: 3 positional and 3 rotational; a rigid body in 2D has 3 degrees of freedom: 2 positional and 1 rotational. A constraint decreases the degrees of freedom of a rigid body. For instance, a constraint that pins an object in space at its center of mass decreases the object’s degrees of freedom by 3: all the positional degrees of freedom are removed and the object can thus now only rotate with 3 degrees of freedom.

In a constraint-based physics engine, we model everything as a constraint: including collision contacts, frictions, springs, pulleys, you name it. A __joint__ is a constraint that involves the interaction of 2 rigid bodies; the examples enumerated here are all technically joints.

Mathematically, a constraint is of the form , where is the expression representing the constraint. We would specifically express as a linear combination of positional properties (linear position and angular position, a.k.a. orientation); for joints, would contain positional properties from both rigid bodies.

### Sequential Impulse

The part of code during the resolution phase that makes sure all constraints are satisfied is referred to as the __constraint solver__. There are numerous ways to perform this task. One of the most popular methods is Sequential Impulse, proposed and popularized by Erin Catto, author of the famous Box2D. Unless specifically stated otherwise, from here on I will be talking about Sequential Impulse in terms of resolution. Later on, I will draw an analogy between the process of solving constraints with Sequential Impulse and point projections.

Directly manipulating position properties of rigid bodies to satisfy constraints is a bad idea, because it would result in excessive jitters and your game would look really bad. Instead, Erin Catto proposed that we derive the expression for , differentiate it to get the corresponding velocity constraint , and then solve for .

A velocity constraint between two rigid bodies is modeled as a linear combination of linear and angular velocities of both rigid bodies (body A and body B) plus a bias term.

where is a 12-by-1 column vector called the __velocity vector__ that contains the linear and angular velocities of both rigid bodies:

and are 3-by-1 column vectors representing linear velocities of rigid body A and B, respectively; and and are 3-by-1 column vectors representing angular velocities of rigid body A and B, respectively.

is a 12-by-1 row vector called the __Jacobian__ that contains the coefficients of linear combination of velocity components:

where , , and are all 3-by-1 column vectors that represent coefficients of linear combination of the components from the velocity vector.

The bias term is only present in certain constraint scenarios. It is a means to “introduce energy” into the system. More on this later.

### Solving Constraints

So, for a single constraint that involves two rigid bodies (a joint), how will you go about and manipulate the linear and angular velocities of both rigid bodies to satisfy it?

First, please take a leap of faith here and believe that the change we need to apply to the velocity vector in order to satisfy the constraint, denoted , is proportional to , where is the 12-by-12 __mass matrix__:

where each is a 3-by-3 zero matrix; is the mass matrix of rigid body A (with mass ), and is the inertia tensor of rigid body A; is the mass matrix of rigid body B (with mass ), and is the inertia tensor of rigid body B.

The inverse of the mass matrix is thus:

Now that we know (or assume at this point) that the required change to the velocity vector is proportional to , in order to satisfy the velocity constraint , lets denote the ratio factor (a scalar) by , called the __Lagrangian Multiplier__. So the required change to the velocity vector is . We know that after applying the change, , to the velocity vector, the constraint is satisfied, so we have the following equation:

and is equal to :

So, with a few algebraic manipulations, we get the following formula:

We call the __effective mass__.

Great. We now have a way to solve for . All is left to do is apply the change, , to the velocity vector, , and then we will have satisfied the constraint .

If we have multiple constraints in the physics system, we then just iteratively repeat this process on all constraints, and the system would converge to a global solution, hence the name Sequential Impulse.

### The Point-Projection Analogy

Remember the leap of faith of believing that is proportional to ? Now I am going to explain the reason behind this with a geometric analogy that involves point projection.

Recall your high school math and look at the equation below:

What does it look like? Right, a plane equation in 3D.

If a point is not on the plane, meaning , the fastest way to bring the point onto the plane is by projecting the point onto the plane, along the plane’s normal .

If we group the coefficients and variables into separate matrices, the plane equation becomes:

where is a 1-by-3 matrix representing the plane normal :

and is a 3-by-1 matrix representing the point :

is just a “bias term”.

The required change to in order to project it onto the plane, denoted , is proportional to the plane’s normal,

Do you see it now? and has the same exact form, except that the former formula deals with a 3-dimension plane, and the latter deals with a 12-dimensional plane. The “shortest” change to the velocity vector, , is along the “plane normal” . However, different rigid bodies have different masses and inertia tensors, so we cannot just apply the change uniformly to the entire velocity vector, hence is modified by an inverse mass matrix, and finally we get the proportion .

Using Sequential Impulse, we iteratively repeat the process of “point projection”. If there are not conflicting constraints, we would eventually converge to the “intersection of all planes”, which is the global solution.

On a side note, if our physics system contains conflicting constraints (say, an impossible set-up of distant joints), then Sequential Impulse cannot find a proper global solution, which is equivalent to trying to find an intersection of planes that don’t actually have a common intersection.

### Solving Multiple Constraints in One Shot

Sometimes, we would like to solve multiple constraints at once for stability considerations. For instance, it is better to solve a prismatic joint as solving a distance joint and angular joint simultaneously rather than solving them separately.

In the aforementioned example, the Jacobian matrix would be a 2-by-12 matrix, instead of a 1-by-12 matrix. And the formula for the Lagrangian Multiplier will change from this:

to this:

where the effective mass is no longer a scalar but a 2-by-2 matrix. And the Lagrangian Multiplier is now a 2-by-1 matrix instead of a scalar.

However, is still a 12-by-1 matrix, because is still a 12-by-1 matrix.

### End of Constraints & Sequential Impulse

Now you understand the theory behind velocity constraints and how to solve them with Sequential Impulse. In the following post, I will go on and derive one of the most important constraints for a physics engine: the contact constraint.

I think what you’re calling the “effective mass” is actually the “inverse effective mass”.

I’ve seen both definitions. I’m using the one I saw on the Bullet forums. I can see how “inverse effective mass” would make more sense, since there’s an inverse mass matrix in the definition.

Hello, i have implemented the principles in your articles and it works great for a pairwise constrain, but how can i join more particles togheter without exploding the simulation ?

I tried to constrain 3 particles togheter and they races agaiinst each other.

Iterative impulses didn’t work either.

Are you trying to “constrain” 3 particles together so that they don’t have relative motion? If so, you should try a technique called “composition”. Essentially, a single rigid body owns the three particles, and the rigid body is integrated as a single object. Each particle stores relative position from the rigid body’s center of mass (computed from the particles) and update its position after the rigid body is integrated, no constraints involved.