Data-Augmented Contact Model for Rigid Body Simulation

TL;DR

This paper presents a data-augmented contact model that combines analytical physics with learned data to improve the realism and expressiveness of rigid body contact simulations, handling complex behaviors efficiently.

Contribution

It introduces a novel contact modeling approach that integrates data-driven learning with analytical physics, enhancing simulation accuracy while maintaining computational efficiency.

Findings

Improved modeling of contact behaviors including bouncing, sliding, spinning.

Enhanced expressiveness over standard Coulomb models through learned transitions.

Maintains simulation speed despite increased geometric complexity.

Abstract

Accurately modeling contact behaviors for real-world, near-rigid materials remains a grand challenge for existing rigid-body physics simulators. This paper introduces a data-augmented contact model that incorporates analytical solutions with observed data to predict the 3D contact impulse which could result in rigid bodies bouncing, sliding or spinning in all directions. Our method enhances the expressiveness of the standard Coulomb contact model by learning the contact behaviors from the observed data, while preserving the fundamental contact constraints whenever possible. For example, a classifier is trained to approximate the transitions between static and dynamic frictions, while non-penetration constraint during collision is enforced analytically. Our method computes the aggregated effect of contact for the entire rigid body, instead of predicting the contact force for each contact point individually, maintaining same simulation speed as the number of contact points increases for detailed geometries. Supplemental video: https://shorturl.at/eilwX Keywords: Physics Simulation Algorithms, Dynamics Learning, Contact Learning

Equations

$${\dot{\mathbf{q}}}_{t+1}={\dot{\mathbf{q}}}_t+\mathbf{M}({\mathbf{q}}_t{)}^{-1}(h\mathbf{\tau }+\mathbf{p}),$$

$${\mathbf{q}}_{t+1}=\mathrm{Integrate}({\mathbf{q}}_t,{\dot{\mathbf{q}}}_{t+1}),$$

$${\dot{\mathbf{q}}}_{t+1}={\dot{\mathbf{q}}}_t+h{\mathbf{M}}^{-1}\mathbf{\tau }+h{\mathbf{M}}^{-1}{\mathbf{J}}^T{\mathbf{f}}_j+{\mathbf{M}}^{-1}\mathbf{p},$$

$${\mathbf{v}}_{t+1}=\mathbf{J}{\dot{\mathbf{q}}}_t+h\mathbf{J}{\mathbf{M}}^{-1}\mathbf{\tau }+h\mathbf{J}{\mathbf{M}}^{-1}{\mathbf{J}}^T{\mathbf{f}}_j+\mathbf{J}{\mathbf{M}}^{-1}\mathbf{p}.$$

$${\hat{\mathbf{v}}}_{t+1}=\hat{\mathbf{J}}{\dot{\hat{\mathbf{q}}}}_t+h\hat{\mathbf{J}}{\hat{\mathbf{M}}}^{-1}\hat{\mathbf{\tau }}-h\hat{\mathbf{J}}{\hat{\mathbf{M}}}^{-1}{\hat{\mathbf{J}}}^T{\mathbf{f}}_j,$$

$$G({\mathbf{f}}_j)={\mathbf{v}}_{t+1}-{\hat{\mathbf{v}}}_{t+1}=\mathbf{A}{\mathbf{f}}_j+\mathbf{J}{\mathbf{M}}^{-1}H({\mathbf{f}}_j)+\mathbf{c}=0,$$

$${\mathbf{f}}_j=(h\hat{\mathbf{J}}{\hat{\mathbf{M}}}^{-1}{\hat{\mathbf{J}}}^T{)}^{-1}\hat{\mathbf{J}}({\dot{\hat{\mathbf{q}}}}_t+h{\hat{\mathbf{M}}}^{-1}\hat{\mathbf{\tau }}).$$

Full text

1 Introduction

The use of virtual simulations has been crucial to the recent progress in building learning-based robot controllers. With recent works showing that improving the rendering realism of simulated perception [Rao et al.(2020)Rao, Harris, Irpan, Levine, Ibarz, and Khansari] or the accuracy of robot motor dynamics [Hwangbo et al.(2019)Hwangbo, Lee, Dosovitskiy, Bellicoso, Tsounis, Koltun, and Hutter] can improve the performance of learned controllers, more accurate modeling of the contact dynamics would also help the robot better understand and interact with the environment it resides. Most recent robot controllers are still developed from off-the-shelf physics simulators which use the idealized Coulomb friction model, an empirical construct to approximate the changes between two physical regimes (static friction vs dynamic friction). However, the Coulomb model assumes linear relationship between normal force and frictional force and uses a single friction coefficient to represent an isometric friction cone. The computation of contact force also involves approximation and arbitrary decisions. For example, many existing simulators formulate a Linear Complementarity Program (LCP), which solutions are not unique except for the frictionless case [Brogliato and Brogliato(1999)]. Depending on the initial guesses and numerical methods used for solving the LCP problem (e.g. Lemke method [Lemke and Howson(1964)] vs Gauss-Seidel algorithm [Jourdan et al.(1998)Jourdan, Alart, and Jean]), the resulting contact forces can be drastically different. These existing issues suggest that a computational contact model grounded by observed data can be a desired alternative.

While recent work has shown that physical phenomena can be learned from data and approximated by neural networks, precisely enforcing constraints, such as contacts, remains difficult for these function approximators learned in an end-to-end fashion. Consider a box resting on a table. If the contact force is slightly larger or smaller than the gravitational force, we will start to see the box rattling or sinking into the table without any external force. Such categorically incorrect prediction of physics outcomes are likely to have negative impacts on the development of control policies.

This paper introduces a data-augmented contact model combining analytical solutions and empirical data collected for a particular scenario (e.g. a specific robot foot colliding a specific surface), such that the simulated trajectories better match the observed data.

Our approach is built on two key insights. First, we utilize analytical solutions from first principles whenever possible and only resort to data-driven approach when the phenomenon is less well understood. For example, the contact force that prevents the objects from interpenetrating is enforced by equations instead of learned from the data. In contrast, we rely on observed data to model and validate the less-understood static friction boundary. To this end, we propose to decompose the contact problem to two steps: predicting the next contact state (i.e. static, dynamic, or detach) and determining contact forces. We solve the first step by learning a classifier from the observed data and the second step by a combination of learning a regressor and solving constrained systems. Second, we propose to compute the aggregated effect of contact at the rigid-body level, instead of predicting the contact force at each individual contact point. This decision removes the concern of simulation speed as the number of contact points increases for real-world objects with detailed geometries.

Our algorithm also applies to articulated systems, distinguishing our framework from most data-driven contact models. Once a rigid body’s contact model is trained, we can simply connect the rigid body to an articulated system without retraining the contact model. This implies that the collision data for learning the contact can be collected using replicas of the disassembled end-effector, without putting the entire robotic system at risk.

As our work focuses on a full 3D scenario in which the object can bounce, slide, or spin in all directions, existing contact data such as MIT Push [Yu et al.(2016)Yu, Bauza, Fazeli, and Rodriguez] are not diverse enough for our evaluation. Instead, we learn the data-driven models from synthesised collision data in simulations. Evaluations show that our data-augmented model matches ground-truth contact behaviors of both single and articulated rigid body systems, and better than a purely statistical model.

2 Related Work

Contact and friction is a common but extremely complex phenomenon, which continuously fascinates generations of scientists and engineers. Since Coulomb and Amontons in the 18th century made the distinction between static and dynamic frictions, inadequacies and controversies of Coulomb’s law have been extensively studied [Popova and Popov(2015)]. For example, [Oden and Martins(1985)] summarized that friction could depend on normal force or stress, normal separation distance, slip displacement, slip velocity, time of stationary contact, slip history, and vibrations. As such, many empirical models that substitute Coulomb’s law have been proposed and are summarized in [Olsson et al.(1998)Olsson, Åström, De Wit, Gäfvert, and Lischinsky]. There is no one model that is more accurate than the others in all scenarios. [Goyal et al.(1989)Goyal, Ruina, and Papadopoulos] proposed the concept of limit surface to enclose all possible friction forces on an object during planer sliding. Our 3D model also works on the rigid-body instead of the point-contact level, and loosens the assumptions that friction must oppose the direction of motion and that friction is isotropic.

In computer animation and robotics, the simplest Coulomb’s law and the point contact representation are often used to approximate contact physics for visualizing or evaluating robotic algorithms. Existing rigid-body simulators often solve a Linear Complementarity Problem (LCP) by enforcing unilateral constraints and complementarity between relative velocities and contact forces [Baraff(1989), Stewart and Trinkle(2000), Anitescu and Potra(2002)]. In contrast, [Todorov(2014)] relaxed contact constraints to solve a convex optimization problem at each time step.

With the emergence of simulation-based controller learning, developing data-driven physics simulators capable of predicting the real world has been actively studied. [Chang et al.(2016)Chang, Ullman, Torralba, and Tenenbaum, Byravan and Fox(2017), Lerer et al.(2016)Lerer, Gross, and Fergus] trained neural networks to learn physical intuitions from vision perception, [Zhou et al.(2016)Zhou, Paolini, Bagnell, and Mason] built a data-efficient model for the limit surface of an object during planar sliding, [Bauza and Rodriguez(2017)] built a probabilistic model for planer sliding that takes into account the stochasticity of frictional forces, and [Fazeli et al.(2017)Fazeli, Zapolsky, Drumwright, and Rodriguez] proposed either to train a purely data-driven model, or to use data to learn the optimal parameters of an analytical model for planer impact. Our work is instead a data-augmented one, incorporating observed data while preserving fundamental contact constraints analytically. Contact detection (not handling) could also be replaced with learnable components via learning contact Jacobian [Pfrommer et al.(2020)Pfrommer, Halm, and Posa] or differentiable shape parameterization [Strecke and Stueckler(2021)]. System identification [Chebotar et al.(2019)Chebotar, Handa, Makoviychuk, Macklin, Issac, Ratliff, and Fox] serves as another approach to improve simulation fidelity without requiring new simulation models, with recent works demonstrating promising results by utilizing gradient-free reinforcement learning [Jeong et al.(2019)Jeong, Kay, Romano, Lampe, Rothorl, Abdolmaleki, Erez, Tassa, and Nori, Jiang et al.(2021)Jiang, Zhang, Ho, Bai, Liu, Levine, and Tan], or making simulators differentiable [Heiden et al.(2021)Heiden, Millard, Coumans, Sheng, and Sukhatme, Le Lidec et al.(2021)Le Lidec, Kalevatykh, Laptev, Schmid, and Carpentier, Degrave et al.(2019)Degrave, Hermans, Dambre, and wyffels, Jatavallabhula et al.(2021)Jatavallabhula, Macklin, Golemo, Voleti, Petrini, Weiss, Considine, Parent-Levesque, Xie, Erleben, Paull, Shkurti, Nowrouzezahrai, and Fidler].

The advantages of integrating neural networks with analytical models have been demonstrated recently. [Ajay et al.(2018)Ajay, Wu, Fazeli, Bauza, Kaelbling, Tenenbaum, and Rodriguez] introduced a method that trains a recurrent neural network to predict the deviation between real-world and simulated contact trajectories, extended in [Fazeli et al.(2020)Fazeli, Ajay, and Rodriguez] for better long-term accuracy. [Kloss et al.(2018)Kloss, Schaal, and Bohg] proposed to combine a neural network for perception with a physics model for prediction, and [Pizzuto and Mistry(2021)] added a non-penetration loss to augment the data loss. Instead of treating the physics engine as a black box and correcting its output, our work directly improves the full 3D contact handling using learned function approximators. We also show that our learned contact models can be reused in new articulated systems without retraining.

3 Method

We propose a method to predict the contact impulse between a specific pair of near-rigid objects. An ideal contact model in a physics simulator should at least guarantee the following properties:

Non-penetration: The geometries of the objects in contact should not overlap. 2. 2.

Repulsive force: The contact force should only push the objects away instead of pulling them together. 3. 3.

Workless condition: The contact force becomes zero at the instance when the bodies begin to separate. 4. 4.

Two friction regimes: There is an unsmooth switch between static and dynamic friction forces, depending on the materials of the objects and other factors. 5. 5.

Dynamic friction model: The dynamic friction force depends on the normal contact force, the relative velocity of the objects, and other factors.

Among these five properties, (4) and (5) depend on empirical models because their mechanics are not well understood. As such, our method will use a data-driven approach to achieve (4) and (5), while maintaining the analytical solution that satisfies (1)-(3).

3.1 Assumptions

We address a perfectly inelastic (i.e. the coefficient of restitution is zero), non-adhesive contact phenomenon between two near-rigid objects, which is a common scenario in robotic applications. We assume both objects have convex shapes and that one of them is stationary. The deformation during the collision is negligible comparing to the overall rigid body motion. We expect the range of collision impulses encountered at test time to lie within the range used for training. For collision between articulated rigid body systems, we assume that the distal link is the only part of the system that is in contact.

We compute the aggregated effect of contact for the entire rigid body, instead of using point-contact representation, which increases computational complexity and leads to issues with over-parameterization. Our method represents the contact geometry as a patch with non-zero area on the surface of the object in contact, denoted as $\mathcal{P}$. We assume that the normal direction of the contact patch is well-defined.

3.2 Single rigid body

The contact computation in physics engines typically consists of two separate processes, contact detection that identifies the contact location on the surface of the object, and contact handling that calculates the contact force such that the contact constraints and the equations of motion are satisfied. Using a standard contact detector, $D(\mathbf{q}_{t})$, we can compute a contact patch $\mathcal{P}$, represented by the convex hull of contact points, given shape and current position of the rigid body $\mathbf{q}_{t}\in\mathbb{R}^{6}$ expressed in the generalized coordinates.

The main challenge of contact computation lies in contact handling. In its most general form, a contact handling routine can be expressed as a function $\mathbf{p}=H(\mathbf{q}_{t},\dot{\mathbf{q}}_{t},\boldsymbol{\tau},\mathcal{P})$, which maps the pre-contact state (i.e. $\mathbf{q}_{t}$ and $\dot{\mathbf{q}}_{t}$), applied forces $\boldsymbol{\tau}$, and contact patch $\mathcal{P}$ to the 6D contact impulse, $\mathbf{p}$. Precisely, $\mathbf{p}$ is the integrated pressure on the contacting surface over the entire contact patch $\mathcal{P}$ and over the time step interval $h$. During the collision process, the contact pressure might not be constant, but its aggregated effect in one time step is equivalent to the impulse $\mathbf{p}$, which can be used to integrate the state forward to the next time step by

[TABLE]

where $\mathbf{M}(\mathbf{q}_{t})$ is the generalized mass matrix of the rigid body. At every time step, if $D$ detects a non-empty $\mathcal{P}$, we invoke the contact handler $H$, described in Algorithm 1.

We first update the current velocity $\dot{\mathbf{q}}_{t}$ to an intermediate velocity by explicitly integrating the applied force $\boldsymbol{\tau}$: $\dot{\mathbf{q}}^{(1)}=\dot{\mathbf{q}}_{t}+h\mathbf{M}^{-1}\boldsymbol{\tau}$. Directly training a regressor to predict the contact impulse $\mathbf{p}$ is likely to violate Properties (1)-(3), as they require precise satisfaction of constraints. Instead, we use observed data to train a classifier $C(\mathbf{q}_{t},\dot{\mathbf{q}}^{(1)})$ that predicts one of the following outcomes for $\mathcal{P}$: static, dynamic, or detach. Based on the predicted outcome, we will calculate the contact impulse differently in order to satisfy (1)-(3).

Static case

The static case indicates that $\mathcal{P}$ will remain in the same position and orientation at the beginning of next time step, but the rigid body is not necessarily stationary. Fig. 2 illustrates an example in which the rigid body is moving while $\mathcal{P}$ is static. The contact impulse in this case must ensure that $\mathcal{P}$ has zero velocity at the end of this time step. In physics engines, this is easily achieved by setting positional constraints at $\mathcal{P}$ and solve for the constraint impulse $\mathbf{p}$ that satisfies $\mathbf{v}_{p}=\mathbf{0}$, where $\mathbf{v}_{p}$ is the generalized velocity of $\mathcal{P}$ at the end of time step. Since $\mathbf{v}_{p}$ and $\mathbf{p}$ have the same effective degrees of freedom, the linear system $\mathbf{v}_{p}(\mathbf{p})=\mathbf{0}$ has a unique solution. Therefore, if a state is classified as ”static” by a highly accurate classifier $C$, the solution of $\mathbf{v}_{p}(\mathbf{p})=\mathbf{0}$ will satisfy Properties (1)-(3).

Dynamic case

In the dynamic case, the contact patch $\mathcal{P}$ will change its position or orientation, or lose some area at the next time step (Fig. 2). To compute the motion of $\mathcal{P}$, we need to predict the contact impulse $\mathbf{p}$. We propose to train a regressor from the observed data because the idealized Coulomb friction model is limited when approximating the complex dynamic friction phenomena. However, directly using a regressor to predict $\mathbf{p}$ will still suffer from the same problem of failing to satisfy Properties (1)-(3) precisely.

Instead, our algorithm first reparameterizes and decouples $\mathbf{p}=(p_{x},p_{y},p_{z},m_{x},m_{y},m_{z})^{T}$, where $p$ and $m$ indicate linear impulse and impulsive torque respectively, into a frictional impulse $\mathbf{p}_{f}=(p_{x},p_{z},m_{y})^{T}$ and a normal impulse $\mathbf{p}_{n}=(p_{y},\tilde{m}_{x},\tilde{m}_{z})^{T}$. We define $\tilde{m}_{x}$ and $\tilde{m}_{z}$ to be the impulsive torques induced by the normal linear impulse $p_{y}$. With this decoupling, in a perfectly inelastic case (i.e. restitution is zero), given any $\mathbf{p}_{f}$, there exists one unique $\mathbf{p}_{n}$ such that Properties (1)-(3) are satisfied 111The actual distribution of normal force over the contact patch is still undetermined. [Baraff(1992)]. Therefore, we train a regressor to only predict $\mathbf{p}_{f}$ and analytically calculate the unique solution for $\mathbf{p}_{n}$ based on the predicted $\mathbf{p}_{f}$. Specifically, Algorithm 1 in the dynamic case first predicts $\mathbf{p}_{f}$ using the trained regressor $R(\mathbf{q}_{t},\dot{\mathbf{q}}^{(1)})$ (Line 8) and then integrates $\mathbf{p}_{f}$ to obtain a second intermediate velocity: $\dot{\mathbf{q}}^{(2)}=\dot{\mathbf{q}}^{(1)}+\mathbf{M}^{-1}\mathbf{T}_{f}\mathbf{p}_{f}$ (Line 9). Here $\mathbf{T}_{f}\in\mathbb{R}^{6\times 3}$ transforms $\mathbf{p}_{f}$ to the generalized coordinates. The unique solution for $\mathbf{p}_{n}$ can be solved by any routine that respects normal complementaries. Our algorithm uses a Danzig-like positive definite LCP solver for frictionless contacts [Cottle and Dantzig(1968”)] (Line 10). Finally, we combine the analytical $\mathbf{p}_{n}$ and predicted $\mathbf{p}_{f}$ to obtain $\mathbf{p}$ in generalized coordinates (Line 11).

If the regressor $R$ were perfectly accurate, the uniqueness of $\mathbf{p}_{n}$ ensures that the decoupling treatment in Algorithm 1 does not affect the true solution of the contact impulse $\mathbf{p}$. When $R$ is not perfectly accurate, the frictionless LCP (Line 10) serves as a corrective step on $\mathbf{p}$ that prioritizes the satisfaction of Properties (1)-(3) over Property (5).

Detach case

In the detach case, $\mathcal{P}$ is predicted to have positive normal velocity and leave the surface at the next time step. We thus set $\mathbf{p}=\mathbf{0}$, ensuring that Property (3) is satisfied.

Remarks: Our method addresses Property (4) by learning a classifier from the observed data. Similarly, the regressor incorporates the observed data to address Property (5). Assuming that the classifier is highly accurate, the resultant contact impulse for the static and detach cases will closely match reality and satisfy Properties (1)-(3) exactly. The classifier will be less accurate near the decision boundary, which coincides with the poorly-understood region where transitions between different physical regimes occur.

3.3 Articulated rigid bodies

Our method can be extended to articulated rigid body systems. We decompose the state of the system into the distal body $(\mathbf{q}_{t},\dot{\mathbf{q}}_{t})$ and all other bodies in the upstream system $(\hat{\mathbf{q}}_{t},\dot{\hat{\mathbf{q}}}_{t})$ (Fig. 3). The joint force $\mathbf{f}_{j}$ transmitted between the distal body and the upstream system is unknown and must be solved simultaneously with the contact impulse $\mathbf{p}$.

Algorithm 4 starts with expressing the velocity of the distal body at the next time step:

[TABLE]

where $\mathbf{M}(\mathbf{q}_{t})$ is the generalized mass matrix for the distal body and $\mathbf{J}(\mathbf{q}_{t})$ is the Jacobian transforming from the generalized coordinates of the rigid body to the Cartesian space at the joint. The Cartesian velocity at the joint at the next time step is then given by

[TABLE]

Similarly, the velocity of the upstream system evaluated at the joint can be expressed as

[TABLE]

where $\hat{\mathbf{M}}$ and $\hat{\mathbf{J}}$ are the mass matrix and Jacobian for the upstream system.

Since the joint constraint is satisfied at the beginning of the time step $t_{0}$, we only need to ensure that the velocity of the constraint is satisfied so that at $t_{1}$ the distal body and the upstream system still coincide at the joint. Therefore, we need to solve for a $\mathbf{f}_{j}$ such that $\mathbf{v}_{t+1}-\hat{\mathbf{v}}_{t+1}=\mathbf{0}$:

[TABLE]

where $\mathbf{A}=h\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^{T}+h\hat{\mathbf{J}}\hat{\mathbf{M}}^{-1}\hat{\mathbf{J}}^{T}$ and $\mathbf{c}=\mathbf{J}\dot{\mathbf{q}}_{t}+h\mathbf{J}\mathbf{M}^{-1}\boldsymbol{\tau}-(\hat{\mathbf{J}}\dot{\hat{\mathbf{q}}}_{t}+h\hat{\mathbf{J}}\hat{\mathbf{M}}^{-1}\hat{\boldsymbol{\tau}})$ are constants in the equation given $\mathbf{q}_{t}$ and $\dot{\mathbf{q}}_{t}$. $H(\mathbf{f}_{j})$ is a shorthand for $H(\mathbf{q}_{t},\dot{\mathbf{q}}_{t},\boldsymbol{\tau}+\mathbf{J}^{T}\mathbf{f}_{j},\mathcal{P})$, which outputs $\mathbf{p}$ depending on $\mathbf{f}_{j}$.

We solve Eq. 6 using Powell hybrid method [Powell(1970)], which uses finite difference to approximate the Jacobian matrix and is less sensitive to the initial guess to the problem. Powell’s method only requires a routine to evaluate $G(\mathbf{f}_{j})$ and an initial guess. Using the heuristic that assumes $\hat{\mathbf{v}}_{t+1}=\mathbf{0}$, we compute the initial $\mathbf{f}_{j}$ by

[TABLE]

In our experiment a solution can always be found at each time step with 0.5% convergence tolerance. The number of evaluations of $G$ is often fewer than 10.

3.4 Implementation

A contact patch $\mathcal{P}$ in the real world will always be a 2D surface. However, when $\mathcal{P}$ degenerates to nearly an edge or a point, in practice, the dimension of the controllable space of the friction impulse will reduce. Since the dimension of $\mathcal{P}$ is available from the collision detector $D$, we utilize this information to improve learning accuracy by treating three types of $\mathcal{P}$ separately: a surface (2D), a line (1D), or a point (0D). For each type, we train a specific classifier and a regressor. Using separate neural networks allows the regressors to have different output dimensions according to the dimension of controllable space of the friction impulse.

The same set of training data can be used to train the classifiers and the regressors. The data collection involves throwing objects to each other with different initial velocities. Since $\boldsymbol{\tau}$ is not part of the input of the learned models, we do not need to apply various $\boldsymbol{\tau}$ during training sample generation, greatly simplifying the data collection process. We record the entire trajectory for each throw and extract the state of every contact instance: $\mathbf{q}_{t},\dot{\mathbf{q}}_{t}$. By evaluating $\dot{\mathbf{q}}_{t}$ at the patch, we can identify and label “static” and “dynamic” cases. To determine “detach” cases and to calculate the training output for the regressors, we need to recover the contact impulse $\mathbf{p}$ for each contact instance: $\mathbf{p}=\mathbf{M}(\dot{\mathbf{q}}_{t+1}-\dot{\mathbf{q}}_{t})-h\mathbf{g}$, where $\mathbf{g}$ is the gravitational force. If $\mathbf{p}$ is near zero, we label this contact instance ”detach”.

The range of initial velocities is chosen to cover the range of the anticipated collision impulses during testing. We found that the choice of rotation representation significantly affects the accuracy of the learned models. Our experiments show that representing the 3D orientation of the rigid body as a rotation matrix outperforms other representations of SO(3). We also found that including two redundant features—the position of the center of $\mathcal{P}$ in the body frame and the velocity at the center of $\mathcal{P}$—reduces the errors of the regressors.

4 Evaluation

We evaluated our data-augmented contact model on rigid bodies with different 3D shapes, a 3-linked articulated rigid body chain, and an object with anisotropic friction coefficient. To demonstrate the complexity of 3D collision, we also included one 2D example for comparison. The ground-truth data were collected in a simulated environment using a generic physics engine DART [Lee et al.(2018)Lee, Grey, Ha, Kunz, Jain, Ye, Srinivasa, Stilman, and Liu].

Although our method is agnostic to the representation of classifiers or regressors, we used feed-froward neural networks for their expressiveness as function approximators. A standard cross-entropy error or MSE was used as the loss function. We used the same collision data to train classifiers and regressors. The hyper-parameters are shown in Fig. 5.

We evaluated our results against two baselines. The first baseline is the ground truth (GT), i.e. contact handling of DART. The second baseline is a purely data-driven (PDD) approach which learns regressors to directly predict impulse from the rigid body state, without using classifiers or solving analytical constraints. The input representation, network architecture and learning algorithm are the same between PDD and our regressors. However, we gave PDD $50\%$ more training data to reach comparable one-step accuracy (Fig. 5).

4.1 Single step prediction

We generated $40,000$ individual collisions to test each learned classifier and regressor. The accuracies are shown in Fig. 5. Note that PDD can reach $90\%+$ accuracy when predicting individual collisions.

4.2 Predicting trajectories of single rigid body

We threw a box and a pentagon prism to the ground under gravity from various initial positions and velocities. The range of initial states and the geometries of the rigid bodies are detailed in Fig. 6, where $h_{0}$, $\mathbf{R}$, $\omega$, $v$ indicate the ranges of the initial height, 3D orientation, angular velocity, and linear velocity respectively. Each simulated trajectory contains $800$ time steps, equivalent to $1.6$ second of motion. For each rigid body, we simulated $100$ trajectories with random initial states and reported three metrics: the average errors of the final horizontal distance, of the final orientation, and of the first collision impulse 222The later collisions cannot be compared to the ground truth because the motions start to deviate after the first collision. of each trajectory.

Fig. 7 shows the results in comparison with the two baselines, where for each metric, we show the mean of GT, the average error of PDD compared to GT, and the average error of our method compared to GT. In most cases, our method matches GT closely and is significantly better than PDD, demonstrating the advantages of using a classifier and analytical solutions. It is worth noting that both our method and PDD achieve similar accuracy in learning the regressors, but our method has much lower error in predicting the impulses. This is because when a collision instance is correctly classified as ”static” or ”detach”, our method solves for an analytical solution which adds no error to the simulation. We also notice that PDD produces large errors in distance and orientation. This is because the small but persistent errors in impulse often result in perpetual movements instead of letting the rigid body come to rest. The erroneous behavior further highlights the advantage of identifying static cases and enforcing analytical constraints for those cases. For comparison, we also tested PDD and our method on a 2D rectangle. PDD performs much better on the planar throwing problem and is comparable to our method. This seems to suggest that PDD can only simulate consecutive bounces well when its regressor has a very high accuracy ($\approx 98\%$), which is much harder to achieve in 3D problems.

Since 3D motions involve much more complex contact behaviors than 2D, we also compared the sequence of contact events (rolling, spinning etc.) in addition to the final state of the trajectory. Fig. 8 shows the contact event sequences of the 3D box from five random throws. To represent a sequence, we used integers (1-8) to label the contacting vertices, lower case letters ($a$ to $l$) to label the contacting edges, and upper case letters ($A$ to $F$) to label contacting faces. The results show that our method produces similar contact events to GT while PDD produces wildly different contact events.

4.3 Predicting trajectories of articulated rigid bodies

We demonstrated our method on an articulated three link system connected by two revolute joints. The top of the first link is pinned to a fixed world space location. The chain started at a horizontal position and swang passively to the ground under gravity. We compared our method to ground truth and showed the motion sequences in the supplementary video. Though our method is only trained on the contact instances between an isolated distal body (i.e. the third link) and the ground, we show that both contact impulses and joint constraint forces can be predicted or solved accurately.

4.4 Anisotropic friction cone

We created a fictitious material which has an anisotropic friction cone. The friction coefficient of the ground is $1.5$ along z-axis and $0.75$ along x-axis. We collected training data from this simulated scenario and learned the classifiers and regressors using the same algorithms. During testing, we threw a box to the ground in eight directions. For each direction, we oriented the initial orientation and velocity to align with the throwing direction. Fig. 9 shows the distances traveled for each direction using our contact model (top) and using the GT simulator (bottom). The results show that our method is able to predict the outcome of collision for anisotropic materials.

5 Conclusion and Limitations

We introduced a data-augmented contact model that predicts the contact behaviors for a particular pair of near-rigid bodies or articulated rigid body systems. We evaluated our method on a set of 3D examples using simulated training data. The promising results indicate that our work could be a first step toward a contact model capable of predicting the 3D contact behaviors of real-world, near-rigid materials.

When collecting the training data from the real world, an isolated distal part of the robot will be used to create the collisions with the surface. The range of collision impulses should cover that of the anticipated collisions during the operation of the full robot. Although the data collection can be conducted in isolation without involving the entire robot, the data-efficiency remains a major concern. Since the dimension of the input and output space is relatively low, it is possible to use other function approximators, such as support vector machines or Gaussian processes, which might be more data-efficient than neural networks.

Our current algorithm assumes that one of the objects is stationary, and does not handle simultaneous contacts of multiple bodies, which limits its usage in many manipulation tasks. As immediate future directions, we plan to extend Algorithm 1 to two moving bodies by including states of both objects as input. Finally, combining our method with learnable contact detectors could further increase expressiveness of the model.