Game Physics: Resolution – Constraints & Sequential Impulse

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 C = 0, where C is the expression representing the constraint. We would specifically express C as a linear combination of positional properties (linear position and angular position, a.k.a. orientation); for joints, C 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 C, differentiate it to get the corresponding velocity constraint \dot{C}, and then solve for \dot{C}.

A velocity constraint \dot{C} 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.

    \[ \dot{C} : JV + b = 0,  \]

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

    \[ V =  {\left[ {\begin{array}{ccc}    \overrightarrow{V_A} \\    \overrightarrow{\omega_A} \\    \overrightarrow{V_B} \\    \overrightarrow{\omega_B} \\   \end{array} } \right]} \]

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

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

    \[ J =  {\left[ {\begin{array}{cccc}   \overrightarrow{J_{V_A}}^T & \overrightarrow{J_{\omega_A}}^T & \overrightarrow{J_{V_B}}^T & \overrightarrow{J_{\omega_B}}^T \\   \end{array} } \right]},  \]

where \overrightarrow{J_{V_A}}, \overrightarrow{J_{\omega_A}}, \overrightarrow{J_{V_B}} and \overrightarrow{J_{\omega_B}} are all 3-by-1 column vectors that represent coefficients of linear combination of the components from the velocity vector.

The bias term b 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 \dot{C} : JV + b = 0 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 V in order to satisfy the constraint, denoted \Delta V, is proportional to M^{-1} J^T, where M is the 12-by-12 mass matrix:

    \[ M =  {\left[ {\begin{array}{cccc}   M_A & 0 & 0 & 0 \\   0 & I_A & 0 & 0 \\   0 & 0 & M_B & 0 \\   0 & 0 & 0 & I_B \\   \end{array} } \right]},  \]

where each 0 is a 3-by-3 zero matrix; M_A is the mass matrix of rigid body A (with mass m_A), and I_A is the inertia tensor of rigid body A; M_B is the mass matrix of rigid body B (with mass m_B), and I_B is the inertia tensor of rigid body B.

    \[ M_A =  {\left[ {\begin{array}{ccc}   m_A & 0 & 0 \\   0 & m_A & 0 \\   0 & 0 & m_A \\   \end{array} } \right]}, \, M_B =  {\left[ {\begin{array}{ccc}   m_B & 0 & 0 \\   0 & m_B & 0 \\   0 & 0 & m_B \\   \end{array} } \right]},  \]

The inverse of the mass matrix is thus:

    \[ M^{-1} =  {\left[ {\begin{array}{cccc}   M_A^{-1} & 0 & 0 & 0 \\   0 & I_A^{-1} & 0 & 0 \\   0 & 0 & M_B^{-1} & 0 \\   0 & 0 & 0 & I_B^{-1} \\   \end{array} } \right]},  \]

Now that we know (or assume at this point) that the required change to the velocity vector is proportional to M^{-1} J^T, in order to satisfy the velocity constraint \dot{C} : JV + b = 0, lets denote the ratio factor (a scalar) by \lambda, called the Lagrangian Multiplier. So the required change to the velocity vector is M^{-1} J^T \lambda. We know that after applying the change, \Delta V, to the velocity vector, the constraint is satisfied, so we have the following equation:

    \[ J(V + \Delta V) + b = 0,  \]

and \Delta V is equal to M^{-1} J^T \lambda:

    \[ J(V + M^{-1} J^T \lambda) + b = 0. \]

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

    \[ \lambda = \frac {-(JV + b)} {J M^{-1} J^T} \]

We call J M^{-1} J^T the effective mass.

Great. We now have a way to solve for \lambda. All is left to do is apply the change, \Delta V = M^{-1} J^T \lambda, to the velocity vector, V, and then we will have satisfied the constraint \dot{C} = JV + b.

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 \Delta V is proportional to M^{-1} J^T? 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:

    \[ ax + by + cz + d = 0 \]

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

If a point (x, y, z) is not on the plane, meaning ax + by + cz + d \ne 0, the fastest way to bring the point onto the plane is by projecting the point onto the plane, along the plane’s normal (a, b, c).

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

    \[ NP + d = 0,  \]

where N is a 1-by-3 matrix representing the plane normal (a, b, c):

    \[ N =  {\left[ {\begin{array}{ccc}   a & b & c   \end{array} } \right]},  \]

and P is a 3-by-1 matrix representing the point (x, y, z):

    \[ P =  {\left[ {\begin{array}{c}   x \\   y \\   z \\   \end{array} } \right]} \]

d is just a “bias term”.

The required change to P in order to project it onto the plane, denoted \Delta P, is proportional to the plane’s normal, N^T

Do you see it now? NP + d = 0 and JV + b = 0 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, \Delta V, is along the “plane normal” J^T. However, different rigid bodies have different masses and inertia tensors, so we cannot just apply the change uniformly to the entire velocity vector, hence J^T is modified by an inverse mass matrix, and finally we get the proportion M^{-1} J^T.

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:

    \[ \lambda = \frac {-(JV + b)} {J M^{-1} J^T} \]

to this:

    \[ \lambda = (J M^{-1} J^T)^{-1} \, (-JV - b)},  \]

where the effective mass J M^{-1} J^T is no longer a scalar but a 2-by-2 matrix. And the Lagrangian Multiplier \lambda is now a 2-by-1 matrix instead of a scalar.

However, \Delta V is still a 12-by-1 matrix, because M^{-1} J^T \lambda 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.

About Allen Chou

Physics / Graphics / Procedural Animation / Visuals
This entry was posted in Gamedev, Physics. Bookmark the permalink.

4 Responses to Game Physics: Resolution – Constraints & Sequential Impulse

  1. Andrew Meadows says:

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

    • Allen Chou says:

      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.

  2. josh says:

    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.

    • Allen Chou says:

      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.

Leave a Reply