Speculative contacts#

Because standard collision detection only inspects geom poses at discrete timesteps, geoms can tunnel through each other whenever their displacement over a single step exceeds their thickness, passing from one side to the other without overlapping at either step. This is especially problematic for simulations with fast-moving bodies, thin geoms, or large timesteps.

Speculative contacts work by speculating contacts that would have occurred between the current and the next timestep. This feature can be enabled using the speculative flag. Speculative contacts work like normal contacts and are treated the same by the solver. They are produced by the following pipeline:

  1. Broadphase and sphere filtering expands the sweep-and-prune AABBs and bounding-sphere tests to account for each fast-moving geom’s linear and rotational motion over the step.

  2. Conservative advancement sweeps candidate geom pairs to determine the time of impact (TOI) \(0 < t < 1\) as a fraction of the timestep.

  3. Line-of-sight raycasting casts a ray between the two geoms to exclude speculative contacts obscured by another geom in the middle.

  4. Effective contact distance inverts the solver’s reference acceleration equation (\(\ari = -b_i (J v)_i - k_i r_i\)) to compute a synthetic contact distance dist so that the constraint solver arrests \((1 - t)\) of the relative normal velocity over the step and lands the geom on the contact surface at \(t = 1\).

Example models demonstrating speculative contacts can be found in model/speculative/.

Enabling speculative contacts adds scene-dependent collision overhead per step, scaling with the number of geoms whose single-step travel exceeds their thickness and the candidate pairs they sweep across. At the same time, preventing tunneling often allows a larger timestep than discrete collision detection would tolerate, reducing overall wall-clock time.

Note

Speculative sweeps extrapolate geom motion using only pre-step velocities (mjData.qvel) and ignore acceleration. This is both because collision detection runs in the position stage before forces and accelerations are evaluated, and because acceleration contributes an \(O(\Delta t^2)\) displacement over a single step compared to the \(O(\Delta t)\) velocity term. Consequently, the timestep \(\Delta t\) cannot be made arbitrarily large: if a slow or stationary geom is accelerated across an obstacle within a single step by a large applied force, gravity, or another body pushing it during the constraint solve, its pre-step velocity will not predict the motion and it can still tunnel.

Broadphase and sphere filtering#

To avoid overhead on slow-moving objects, a geom is only treated as fast-moving if the maximum distance traveled by any point on its surface over the timestep \(\Delta t\) exceeds the minimum dimension of its local AABB half-extents. Given pre-step linear displacement \(\Delta x = \Delta t \, v\) and angular displacement \(\Delta \omega = \Delta t \, \omega\) about the geom center \(c\), any point \(p\) on the geom surface moves over the step by

\[\Delta p = \Delta x + \Delta \omega \times (p - c) \]

Since \(\|p - c\| \le r_{\text{bound}}\), where \(r_{\text{bound}}\) is the geom’s bounding-sphere radius, this surface travel is bounded above via the triangle inequality by the sweep distance

\[s = \|\Delta x\| + \|\Delta \omega\| \, r_{\text{bound}} \]

Geoms whose sweep distance \(s\) is smaller than the minimum dimension of their local AABB half-extents cannot tunnel in a single step and use standard discrete bounds.

Swept bounding-sphere filter#

For each candidate geom pair, the bounding-sphere filter tests whether the two bounding spheres can overlap at any fraction \(t \in [0, 1]\) of the timestep. Writing \(p = x_2 - x_1\) for the initial center separation and \(v = \Delta x_2 - \Delta x_1\) for the relative linear displacement over the step, the relative center trajectory is \(r(t) = p + t v\). We find the global minimizer \(\tau\) of \(\|r(t)\|^2\) on \([0, 1]\) by finding the root \(t^*\) of the derivative \(\frac{d}{dt}\|r(t)\|^2 = 2(p + t v) \cdot v = 0\):

\[\tau = \begin{cases} 0, & \phantom{0 \lt {}} t^* \le 0 \\ t^*, & 0 \lt t^* \lt 1 \\ 1, & 1 \le t^* \end{cases} \qquad \text{where } t^* = -\frac{p \cdot v}{\|v\|^2} \]

The pair is pruned if \(\|r(\tau)\| > r_{\text{bound},1} + r_{\text{bound},2} + \text{margin} + \text{gap}\). When one of the geoms is a plane with unit normal \(n\), the filter checks the minimum signed distance at the endpoints \(t \in \{0, 1\}\), pruning the pair if \(\min(d_0,\, d_0 + n \cdot \Delta x) > r_{\text{bound}} + \text{margin} + \text{gap}\).

Swept broadphase AABBs#

During sweep-and-prune (Selection), each geom’s oriented bounding box is projected onto the three principal axes \(\hat{u}_j\) of the scene covariance frame. For a fast-moving geom, the static half-projection along \(\hat{u}_j\) is first inflated by its rotational sweep radius \(\|\Delta \omega\| \, r_{\text{bound}}\) and intersected with its bounding sphere radius \(r_{\text{bound}}\) (since rotation about the geom center cannot move any point outside \(r_{\text{bound}}\)). The resulting interval \([\text{lo}_j, \text{hi}_j]\) is then extended directionally by the linear displacement projection \(\Delta x \cdot \hat{u}_j\):

\[\big[\text{lo}_j + \min(0, \Delta x \cdot \hat{u}_j) - \text{margin},\; \text{hi}_j + \max(0, \Delta x \cdot \hat{u}_j) + \text{margin}\big] \]

Because per-body bounding volume hierarchies (mjModel.bvh_aabb) are precomputed statically in the body frame at compile time, mid-phase BVH tree traversal is bypassed for body pairs containing a fast-moving geom so that swept geom pairs are not culled by uninflated internal tree nodes.

Conservative advancement#

There are many approaches to computing the time of impact (see Erin Catto’s GDC 2013 presentation for an overview). In MuJoCo, when discrete narrowphase collision detection finds no active surface contact (dist < margin) and at least one geom in the pair exceeds the fast-sweep threshold, we use a unilateral, directed conservative advancement algorithm to find the time of impact \(t \in (0, 1)\).

At each iteration, GJK evaluates the minimal distance \(\text{dist}(t)\) between the two geoms, closest witness points \((x_1(t), x_2(t))\), and separating unit normal \(n = (x_2(t) - x_1(t)) / \text{dist}(t)\) (pointing from geom 1 to geom 2) at the current advanced poses. By the Separating Axis Theorem (SAT), the supporting planes perpendicular to \(n\) at \(x_1(t)\) and \(x_2(t)\) bound every pair of points \(p_1(t) \in \text{geom } 1\) and \(p_2(t) \in \text{geom } 2\) along \(n\):

\[(p_2(t) - p_1(t)) \cdot n \;\ge\; (x_2(t) - x_1(t)) \cdot n \;=\; \text{dist}(t) \]

As each geom translates and rotates about its center \(c_i\) with step displacements \((\Delta x_i, \Delta \omega_i)\), any point \(p_i\) on geom \(i\) moves with velocity \(\dot{p}_i = \Delta x_i + \Delta \omega_i \times (p_i - c_i)\) with respect to \(t\). Because every point on geom \(i\) lies within its bounding sphere of radius \(r_{\text{bound},i}\) centered at \(c_i\) (\(\|p_i - c_i\| \le r_{\text{bound},i}\)) and \(\|n\| = 1\), each rotational term satisfies \(\big|n \cdot (\Delta \omega_i \times (p_i - c_i))\big| \le \|\Delta \omega_i\| \, r_{\text{bound},i}\). Projecting the relative closing velocity \(\dot{p}_1 - \dot{p}_2\) onto the fixed axis \(n\) therefore gives an upper bound \(v_{\text{bound}}\) on how fast any two points on the geoms can close across the separating plane:

\[\begin{aligned} n \cdot (\dot{p}_1 - \dot{p}_2) &= n \cdot (\Delta x_1 - \Delta x_2) + n \cdot \big(\Delta \omega_1 \times (p_1 - c_1)\big) - n \cdot \big(\Delta \omega_2 \times (p_2 - c_2)\big) \\ &\le n \cdot (\Delta x_1 - \Delta x_2) + \|\Delta \omega_1\| \, r_{\text{bound},1} + \|\Delta \omega_2\| \, r_{\text{bound},2} \;\equiv\; v_{\text{bound}} \end{aligned} \]

Unlike the worst-case isotropic sweep \(s_1 + s_2\)—which bounds linear motion in all directions and stalls when a geom moves rapidly tangentially while grazing a surface—projecting \(\Delta x_1 - \Delta x_2\) onto \(n\) discards tangential linear motion parallel to the separating plane.

If \(v_{\text{bound}} \le 0\), the geoms cannot close along \(n\) and the search terminates without contact. Otherwise, advancing by \(\Delta t \ge 0\) can reduce the projection along \(n\) by at most \(\Delta t \, v_{\text{bound}}\), so the true separation at \(t + \Delta t\) is bounded below by

\[\text{dist}(t + \Delta t) \;\ge\; \min_{p_1, p_2}\,(p_2(t + \Delta t) - p_1(t + \Delta t)) \cdot n \;\ge\; \text{dist}(t) - \Delta t \, v_{\text{bound}} \]

Setting this lower bound equal to a target clearance \(\text{dist}_{\text{target}} = \text{margin} + 0.02\,(\text{dist}(0) - \text{margin})\) just above margin (\(\text{dist}(t) - \Delta t \, v_{\text{bound}} = \text{dist}_{\text{target}}\)) yields the largest guaranteed-safe step:

\[\Delta t = \frac{\text{dist}(t) - \text{dist}_{\text{target}}}{v_{\text{bound}}}, \qquad t \leftarrow t + \Delta t \]

Unlike bilateral advancement—which can overshoot into overlap and step both forward and backward using signed penetration depth queries—this unilateral advancement approaches \(\text{dist}_{\text{target}}\) strictly from the separated side (\(\text{dist}(t) > \text{margin}\)), so GJK always returns well-defined witness points and a separating normal without invoking EPA. When multiccd is enabled (or one geom is a plane), polygon clipping extracts up to four contact points at the advanced impact pose.

Line-of-sight raycasting#

The witness points \(x_1(t)\) and \(x_2(t)\) found at impact time \(t\) are mapped from each geom’s advanced local frame back to its \(t = 0\) world pose \(x_1(0)\) and \(x_2(0)\), and interpolated by each geom’s sweep fraction \(w_i = s_i / (s_1 + s_2)\) to form the contact position:

\[x_{\text{contact}} = w_1 \, x_1(0) + w_2 \, x_2(0) \]

Because narrowphase evaluates each geom pair independently, a fast-moving geom could otherwise generate a speculative contact with a second object sitting behind an intervening obstacle (for example, a falling box registering a speculative contact with the floor underneath a thin tabletop it is about to strike). To enforce line-of-sight, a ray segment is cast between the \(t = 0\) witness points \(x_1(0)\) and \(x_2(0)\) and tested against all other collidable geoms compatible with both bodies; if any intervening geom intersects the segment, the occluded speculative contact is discarded.

Effective contact distance#

Because a speculative contact is created at \(t = 0\) before geometric overlap occurs, its true separation exceeds margin. Passing the true distance to the constraint solver would either leave the contact inactive or, if activated with zero penetration, cause the velocity damper in the reference acceleration

\[\ari = -b_i (J v)_i - k_i r_i \]

to arrest \(100\%\) of the closing speed immediately—freezing the body in mid-air at \(t = 0\).

Instead, collision detection inverts this reference acceleration equation to solve for an effective contact distance \(\text{dist} = \text{margin} - \delta\) such that the constraint solver arrests \((1 - t)\) of the incoming normal speed \(|v_n|\) over the timestep, leaving the remaining fraction \(t\) of the velocity to carry the geom onto the contact surface at the end of the step (\(z_1 = z_0 + h v_1 = 0\)).

Under the standard Euler integrator, a contact with dimensionless single-step stiffness \(k_0 = h^2 k\) and damping \(b_0 = h b\) (derived from solref and solimp) produces a single-step normal velocity change

\[h \ari = b_0 |v_n| + k_0 \frac{\delta}{h} \]

Setting \(h \ari = (1 - t)|v_n|\) and solving for the effective penetration depth \(\delta\) yields

\[\delta = \frac{1 - t - b_0}{k_0} \, h |v_n|, \qquad \text{dist} = \text{margin} - \delta \]

Under the implicit discrete integrator, where the spring-damper is evaluated at the end-of-step position and divided by the implicit factor \(1 + b_0 + k_0\), the numerator \(1 - t - b_0\) is replaced by \(1 - t - t(b_0 + k_0)\):

\[\delta_{\text{discrete}} = \frac{1 - t - t(b_0 + k_0)}{k_0} \, h |v_n| \]

Notice that for default critically damped contacts (\(b_0 = 1\)), \(1 - t - b_0 = -t < 0\), so \(\delta < 0\) and \(\text{dist} > \text{margin}\): the spring term pulls back against the damper by \(t |v_n|\) so that only \((1 - t)|v_n|\) is removed during the approach step. Speculative contacts remain active even when \(\text{dist} \ge \text{margin}\).

Modeling notes#

  • Elastic restitution: Standard speculative contacts split impact across two steps—arresting \((1 - t)|v_n|\) to land on the surface on the first step, and resolving surface contact on the next. When modeling nearly elastic collisions via direct negative solref (e.g. \(b_0 \approx 2\) so the damper reverses \(v_n \to -v_n\) in a single step, as in newton.xml), set the stiffness to zero (e.g. solref="0 -499") with a stiff impedance (solimp="0.995 0.995 0.001") so the speculative position spring term does not cancel out the elastic damper during the approach step (noting that zero stiffness provides no static position spring for resting contacts).

  • Single-pass pre-step velocities: Because speculative sweeps use pre-step velocities (mjData.qvel), a stationary intermediate body struck by a fast neighbor acquires its post-impact velocity only after the solver runs. In multi-body impact chains (such as Newton’s cradle), keep a small clearance between adjacent bodies or choose a timestep where a single step’s travel does not jump across two gaps at once.