Position Based Dynamics
Published:
Position Based Dynamics
A brief introduction to position based dynamics (PBD).
Classical Simulation and Examples
Physically-based simulation of mechanical effects is an important topic in computer graphics. We want to simulate the contact, collision and other dynamics of shapes. Classical methods simulate the dynamic of shapes with explicit Newton’s second law.
\[F = m \dot{v}\]For element in the system, the position, momenta and forces are all explicitly maintained and updated during the time integration.
If we just want to use the classical method to simulate the motion of a charged particle in a given electromagnetic field, everything works well. But once collision, or in another view, constraints are involved, the numerical simulation becomes difficult.
Example 1: Harmonic Oscillators
Let us consider a harmonic oscillator. The force is explicitly modeled by
\[F = -kx.\]As the simplest physical model, this problem can be solved explicitly. By Newton’s second law, this is a second order differential equation with general solution
\[x(t) = A \cos(\omega t + \phi),\]where $\omega = \sqrt{k / m}$ is the frequency and $\phi$ is the initial phase.
Consider how to simulate this process numerically. At each time step, we maintain the current configuration $(x, v)$. The acceleration is
\[\dot{v} = \frac{F}{m} = \frac{-k x}{m}\]Then the configuration is updated by the semi-implicit Euler method:
\[v^{n + 1} = v^n + \dot{v} \Delta t = v^n - \frac{k x^n}{m} \Delta t,\] \[x^{n + 1} = x^n + v^{n + 1} \Delta t.\]Example 2: Two Particles Connected by a Rod
Consider a more complex scenario. We have two charged particles connected by a rigid rod in a 3-dimensional electric field. Suppose the length of the rod is $L$.
The system is defined over the configuration space $(x_1, x_2, v_1, v_2)$. Different from the harmonic oscillator, there is a constraint imposed by the rigid rod, applying internal forces to the particles.
To simulate the evolution of this system, we model it by a “soft string” with stiffness $k$ (in some sense, modeled as a harmonic oscillator). The internal force is
\[F_1 = k (\Vert x_1 - x_2 \Vert - L) \frac{x_2 - x_1}{\Vert x_2 - x_1 \Vert}.\]We can use the same semi-implicit Euler method to update the configurations.
A new challenge appears. For example, if the stiffness value is too large, the rod may be flipped due to improperly large forces.
Example 3: Rigid Body Collision
Consider the simulation of a set of rigid bodies. During the simulation, the solids may collide and we can detect it. Once solids collide, we apply an impulsive force to them. The force will update the velocity and hence the position.
Like Example 2, numerical overshooting issues can arise.
Position Based Dynamics (PBD)
Different from force-based simulation, position-based dynamics does not explicitly calculate internal forces. PBD maintains only the system configuration $(x_i, p_i)$ and resolves constraints via projection.
In geometric terms, the constraints define a submanifold of the configuration space. The projection step projects the state from the total configuration space to the constrained submanifold.
In the original paper by Müller et al. (2006), the PBD algorithm proceeds as follows:
- Initialization: Set initial positions $x_i^0$, velocities $v_i^0 = 0$, and masses $m_i$.
- Predict: For each timestep, predict new positions: $p_i = x_i + \Delta t (v_i + \Delta t g)$, where $g$ is external acceleration (e.g., gravity).
- Constraint Projection: Iteratively project positions to satisfy all constraints $C_j(x) = 0$.
- Velocity Update: $v_i^{n+1} = \frac{x_i^{n+1} - x_i^n}{\Delta t}$.
- Position Update: $x_i^{n+1} = p_i + \Delta x_i$.
Constraints Solver
In the iteration, we update the positions from $x_i$ to $p_i$, and then project $p_i$ to satisfy
\[C(p_1, \cdots, p_n) = 0\]Now, given $p$, we want to find $\Delta p$, such that $C(p + \Delta p) = 0$. We use a first order approximation
\[C(p + \Delta p) = C(p) + \nabla_p C(p) \cdot \Delta p\]We specifically choose the direction as $\Delta p = \lambda \nabla_p C(p)$, since it is perpendicular to the level set $C(p) = \text{const}$.
With this strategy, we have \(\Delta p = - \frac{C(p)}{\Vert \nabla_p C(p) \Vert} \nabla_p C(p).\)
Stiffness
We can incorporate stiffness $k$ for constraints. The simplest approach is to multiply the corrections $\Delta p$ by $k \in [0, 1]$. There is an exponential corrector presented in the PBD paper; we omit those details here.
Inequality Constraints: Beyond equality constraints $C(x) = 0$, PBD can handle inequality constraints $C(x) \geq 0$ by only applying corrections when the constraint is violated (i.e., $C(x) < 0$). This enables contact handling and collision response in a unified framework.
Extended Position Based Dynamics (XPBD)
XPBD addressed the problem of time-step and iteration count dependent stiffness in PBD.
Define the so-called compliant $\alpha$ as the inverse of stiffness: \(\alpha = \frac{1}{k}.\) We define a potential energy
\[U(x) = \frac{1}{2} C(x)^T \alpha^{-1} C(x),\]which defines an elastic force
\[f = - \nabla_x U = \nabla C(x)^T \alpha^{-1} C(x)\]Consider the equation of motion \(M \ddot{x} = \nabla C(x)^T \alpha^{-1} C(x)\)
By the implicit Euler method, we have
\[v^{n + 1} = v^n + \Delta t M^{-1} \nabla_x C(x^{n + 1})^T \alpha^{-1} C(x^{n + 1})\] \[x^{n + 1} = x^n + \Delta t v^{n+1} = x^n + \Delta t v^n + \Delta t^2 M^{-1} \nabla_x C(x^{n+1})^T \alpha^{-1} C(x^{n+1})\]Define $x^* = x^n + \Delta t v^n$ as the predicted position.
\[M (x^{n + 1} - x^*) - \Delta t^2 ~ \nabla_x C(x^{n + 1})^T \alpha^{-1} C(x^{n + 1}) = 0\]Define
\[\lambda = \Delta t^2~ \alpha^{-1} C(x) = \tilde{\alpha}^{-1} C(x),\]the discretized equation of motion is
\[\begin{aligned} M(x^{n + 1} - x^*) - \nabla C(x^{n + 1})^T \lambda^{n + 1} =: h(x, \lambda) = 0 \\ C(x^{n + 1}) + \tilde{\alpha} \lambda^{n + 1} =: g(x, \lambda) = 0 \end{aligned}\]Linearizing these equations, we have
\[\begin{pmatrix} K & -\nabla C^T(x_i) \\ \nabla C(x_i) & \tilde{\alpha} \end{pmatrix} \begin{pmatrix} \Delta x \\ \Delta \lambda \end{pmatrix} = - \begin{pmatrix} h(x_i, \lambda_i) \\ g(x_i, \lambda_i) \end{pmatrix}\]where $K = \partial h / \partial x$.
Approximation: XPBD makes the following approximations:
- $K = M$, with $O(\Delta t^2)$ error
- $g = 0$
Finally, we arrive at the position update rule
\[\Delta x = M^{-1} \nabla C(x)^T \Delta \lambda\]And the multiplier update rule: \(\Delta \lambda = \frac{-C(x) - \tilde{\alpha} \lambda}{\nabla C M^{-1} \nabla C^T + \tilde{\alpha}}\)