In the previous post, I showed how enforcing $\mathbf{Jv}=0$ imposes one or more constraints. If you have gone ahead and implemented this, you will notice that the system often collapses (especially when gravity is involved) - once the constraint is slightly violated, it never recovers, and the associated rigidbodies plummet down. However, consider what would happen if we were to add in a “bias” term (should be familiar to anyone who’s tried their hand at machine learning): $\mathbf{Jv+b}=0$. Roughly, this allows for less leeway in the positional drift, since we can feed the position error (which we had initially discounted, opting straight for velocity analysis) back into the velocity constraint. Additionally, a time-dependent bias allows us to incorporate objects with prescribed motion, like motors.
It is clear that $\boldsymbol\lambda$ will now be solved as
Let’s focus on the position error. This is represented by the violation of the original constraint $C({q_i})\ne0$, which we would like to indirectly enforce through the velocity. However, to account for the positional drift, we can set the bias to be $\mathbf b=\frac{C}{\Delta t}$, where $\Delta t$ is added into the denominator to make the expression dimensionally consistent (this is extremely important in physics engines, since computationally, missing factors of $\Delta t$ will cause $\mathcal O(10^{2})$ variations in values!). Why is this?
Look at the ordinary differential equation in terms of the constraint:
In other words, if the constraint is initially violated by $C_0$, then the bias term will assist in damping this value to zero (the desired value), exponentially, rather than letting positional errors accumulate. This is known as Baumgarte stabilisation. The bias term is typically equipped with a control parameter $\beta\in[0,1]$ to yield $\mathbf b = \frac{\beta C}{\Delta t}$, and this determined through trial and error, and may even be dependent on the joint type.
An implementation along these lines will be extremely sturdy in general. But the devil lies in contact constraints, which enforce collisoin resolution. Let’s take a look.
Contact constraints are an example of inequality constraints, since clearly we want to objects to have a penetration depth of anything less than zero. In practice, determining the penetration depth for two generic convex colliders is very complicated, but we can blackbox that for now - or alternatively work with a simplified system in which we only have circles (to avoid cheating and not implementing rotation, imagine each circle has a line drawn from the radius to a point on the circumference).
The applied velocity must be such that the bodies are just separated, with zero relative velocity along the normal, upon separation (restitution, i.e. elastic collisions, will be implemented later with minimal changes). This can be implemented through the position constraint
where $\vec p_1$ and $\vec p_2$ are the deepest penetrating points of the two rigibodies, and $\vec n$ is the normal vector separating the two bodies. In other words, a translation of either body through a displacement of $\lVert\vec p_1 -\vec p_2\rVert\vec n$ would just about separate the two - but of course we wish to implement this at the velocity level. In doing so, we can eyeball the Jacobian (this is the general procedure for any constraint).
Here the product rule was used, but the value of $\frac{\mathrm d\vec n}{\mathrm dt}$ is $0$ since the translation will occur along the normal and hence it will remain constant. Importantly, this expression does not involve the motion of the centre of mass $\vec x_1$: it is instead the velocity of another point $\vec p_1$, typically on the boundary of that collider! This means we also have to account for rotational motion, and hence $\vec v_1\vert_{\vec p_1}=\vec v_1\vert_{\vec x_1}+\vec\omega_1\times(\vec p_1-\vec x_1)$, by definition. Thus the velocity constraint becomes
Finally, using the vector identity $(A\times B)\cdot C=(B\times C)\cdot A$, the previous equation can be converted to the form
The Jacobian becomes manifest, without having to calculate a vector of partial derivatives! Finally, we can introduce our newfound Baumgarte stabilisation to regulate the position error. Adding restitution is straightforward from here: the velocity along the normal should not be zero after the update (inelastic collision, as we have done here), but rather a fraction of the incoming relative velocity during collision. In other words
The $\min(e_1, e_2)$ refers to the lower value of the restitutions of the two colliders. This is conventional, and other schemes like the average between the two can also be employed (but clearly $e\in[0,1]$).
Coupled with the solvers described in the previous post, we now have a reasonably stable, fast and accurate engine that can handle all sorts of constraints, including collisions. However, the collisions can be a bit shaky and unreliable, so I will describe some creative methods to deal with this in a future post. Stay tuned!