For reasons of stability and accuracy, we typically work with velocity-based constraints to introduce interesting features into our physics engine, like joints and collisions. This approach outstrips its peers, since:
- Acceleration-based updates lead to extreme instabilities as the accumulated acceleration can go to infinity
- Position-based updates either look jittery or unrealistic (“hugging”), without a good balance in between. You can’t even implement friction at a velocity level.
These constraints are well-defined, and are essentially 100% physically accurate. That’s why the same story carries forth to rendered physics, realistic simulations and robotics. The key difference is that game physics engines are time-critical, with ever-increasing demands for “higher FPS” by gamers leading to ever-decreasing amounts of self-contentment in low-level game programmers. In other words, it’s not the constraints themselves that are the problem, but rather solving them - the time complexity of obtaining the exact solution is cubic in the number of constraints, which is a nightmare for any physics-based game after the 1990s. As the programmer’s saying goes: “It’s $\mathcal O(n)$ or nothing”.
So what is a constraint, mathematically? We have at our disposal the kinematic variables ${q_i}$ of all the rigid bodies in the scene. In practice, $q_i \equiv{\vec x_i, \vec v_i, \vec\theta_i, \vec\omega_i}$ (okay, the vectors on the angular variables are gratuitous for 2D. Whatever, it looks cooler). Then a constraint is simply $C({q_i})\overset!=0$ (equality constraint) or $C({q_i})\overset!\ge0$ (inequality constraint). For instance, a ball-and-socket or fixed-distance joint between $\mathbf 1$ and $\mathbf 2$ is given by $C(\vec x_1, \vec x_2)=\lVert\vec x_1-\vec x_2\rVert=0$. Nothing to it! Unfortunately, apropos the disclaimer I issued above, rigidly enforcing $C=0$ is a mistake, so we work with $\partial_t C=0$ instead. In words, we want the constraint to remaining unchanging over time, and correct the velocity of the rigidbodies in question in order to enforce this. Through the kindergarten multivariable chain rule, we have
Specialising to constraints between two rigidbodies (i.e. practically all of them), we can adapt this to a more linear algebraic form by defining the the stacked velocity vector $\mathbf v = [\vec v_1\ \vec\omega_1\ \vec v_2\ \vec\omega_2 ]$ - this is a column vector with length 12 (or 6 in 2D). Similarly, $\mathbf J\equiv\frac{\partial C({q_i})}{\partial q_i}$, which is a row vector or the same length, known as the Jacobian (a familiar term from multivariable calc). If you are aware of Jacobians acting roughly as coordinate transformations, then you can view it here as a projection map down from ordinary Euclidean space to constraint space. Then equation $(1)$ becomes $\mathbf{Jv}=0$. We can interpret this as meaning that the Jacobian vector is perpendicular to the velocity. But this is precisely the same property obeyed by the constraint force!
Why? Recall the principle of virtual work. Since constraint forces are nudging the system back into static equilibrium, they crucially cannot do any work. Consequently for the constraint forces, $\mathbf F\cdot\mathbf v=0$, which implies that $\mathbf F=\mathbf J^\top\lambda$, where the transpose is merely due to the distinction between row and column vectors, and $\lambda$ acts as a Lagrange multiplier - it is this value we need to solve to obtain the requisite update, since the value of the Jacobian can be easily computed using the current velocities and the type of constraint.
Now at any given time step, the constraint will typically be violated due to the presence of external forces acting on the body, say due to user input, powered motion, etc. A correctional velocity must be applied to force $\partial_tC$ back to $0$. The velocity update (this is shorthand for the updates to the velocity and angular velocity of both bodies) is by definition $\delta\mathbf v=\mathbf M^{-1}\mathbf F=\mathbf M^{-1}\mathbf J^\top\lambda$. Here $\mathbf M=\mathrm{diag}[m_1\mathbb{I}^3, \mathbf{I}_1, m_2\mathbb{I}^3, \mathbf{I}_2]$, where $m$ and $\mathbf I$ are the mass and inertia tensor respectively - this definition should be intuitively obvious from the definition of the generalised velocity. Thus $\lambda$ is determined by
In other words, at each time step, for each constraint, we update $\mathbf v$ of the involved bodies as
The problem here of course is that updating one constraint at a time generically breaks all of the other ones, so naïvely it appears we cannot enforce all the constraints: are we doomed?
Of course not. The matrix notation above is very suggestive. Rather than having a length-12 row vector as a Jacobian, we can simultaneously incorporate multiple constraints by horizontally stacking the two individual Jacobians to form a $2\times12$ matrix. $\lambda$ is also promoted to a matrix, but $\delta\mathbf v$ remains a column vector, as you can check. Thus, $\mathbf{Jv}=-\mathbf{JM}^{-1}\mathbf J^\top\boldsymbol\lambda\sim\mathbf{b}=\mathbf A\boldsymbol\lambda$ is now a matrix equation in $\boldsymbol\lambda$. As with any other matrix equation, it is best to solve it using something like LU decomposition rather than inverting $\mathbf A$, which would be extremely inefficient.
This will essentially yield perfect results, provided a global solution does indeed exist for all the constraints. However, this method is downright terrible for game physics engines! This is because the size of the stacked Jacobian scales as $\mathcal O(N^2)$, where $N$ is the number of bodies, and so you will have to solve extremely large matrix equations to determine $\boldsymbol\lambda$, with a time complexity of $\mathcal O(N^{\sim2.8})$. Also, from experience, it can lead to single-frame collapse, where there is no global solution, so it just ends up destroying all the constraints in the absence of time-consuming mitigation tactics. That said, this method is perfect for pre-built rendered physics, where speed is not a major concern.
However this is just one method. The second, espoused by Box2D creator Erin Catto, is sequential impulse, where we solve the constraints sequentially and iteratively, in the hope that after several runs, the accumulated velocity correction converges to the global solution. The theoretical time complexity is $\mathcal O(N)$ space and time, though as with most Big-O notations, this hides a constant prefactor. The linear time is because there are no matrix equations to be solved, just fast, constant-size matrix multiplications. This works fantastically … provided you iterate enough times - in a naïve implementation, this may still be too slow. Nevertheless, plenty of work has been channeled into optimisation and stability schemes in this regime, which makes it the most popular framework today.
However, there are certain cases where you actually do want to solve multiple constraints simultaneously for stability reasons - the best example I can think of is a resting inequality constraint: there are two contact normals generated with the following constraint: $C=(\vec x_2-\vec x_1)\cdot\vec n=0$, but solving them sequentially will lead to the top rigidbody oscillating back and forth rather than remaining still. To counter this, the block-solving method identifies “blocks” of objects that you want to solve together as a matrix equation. Then, each block is solved iteratively as in sequential impulse above, bringing together the advantages off both previous methods.
All of these methods, in their raw forms, still raise questions of stability. In the next post, I will satiate the enthusiastic reader with concrete examples of the concepts and methods discussed above, and also describe inequality constraints and various methods to improve stability.