Learning to Slide Unknown Objects with Differentiable Physics Simulations

I Introduction

Nonprehensile manipulation of objects is a practical skill used frequently by humans to displace objects from an initial configuration to a desired final one with a minimum effort. In robotics, this type of manipulation can be more advantageous than the traditional pick-and-place when an object cannot be easily grasped by the robot, due to the design of the end-effector and the size of the object, or the obstacles surrounding the manipulated object. For example, combined pushing and grasping actions have been shown to succeed where traditional grasp planners fail, and to work well under uncertainty .

The mechanics of planar pushing was extensively explored in the past . Large datasets of images of planar objects pushed by a robot on a flat surface were also recently presented . Recent techniques for planar sliding mechanics focus on learning data-driven models for predicting the motions of the pushed objects in simulation. While this problem can be solved to a certain extent by using generic end-to-end machine learning tools such as neural networks , model identification methods that are explicitly derived from the equations of motion are generally more efficient . A promising new direction is to directly differentiate the prediction error with respect to the model of the object’s mechanical properties, such as its mass and friction distributions, and to use standard gradient descent algorithms to search for values of the properties that reduce the gap between simulated motions and observed ones. Unfortunately, most popular physics engines do not natively provide the derivatives of the predicted poses , and the only way to differentiate them is through numerical finite differences, which are expensive computationally.

In the present workThis work was supported by NSF awards 1734492, 1723869 and 1846043. , we propose a method for safely sliding an unknown object from an initial configuration to a target one. The proposed approach integrates a model identification algorithm with a planner. To account for non-uniform surface properties and mass distributions, the object is modeled as a large set of small cuboids that may have different material properties and that are attached to each other with fixed rigid joints. A simulation error function is given as the distance between the centers of the objects in simulated trajectories and the true observed ones. The gradient of the simulation error is used to search for the object’s mass and friction distributions.

The main contribution is the derivation of the analytical gradient of the simulation error with respect to the mass and friction distributions using the proposed cuboid representation of objects. The second contribution is the use of the derived analytical gradients for identifying models of unknown objects by using a robotic manipulator, and demonstrating the computational and data efficiency of the proposed approach. The third contribution is the use of the proposed integrated method for pre-grasp sliding manipulation of thin unknown objects.

A video of the data collection process, the pre-grasp sliding and grasping experiments using the robotic setup shown in Figure 1 can be viewed at https://bit.ly/371Y6Y1 .

II Related Work

Algorithms for Model-based reinforcement learning explicitly learn the unknown dynamics, often from scratch, and search for an optimal policy accordingly . The unknown dynamics are often modeled using an off-the-shelf statistical learning algorithm, such as a Gaussian Process (GP) , or a neural network . This approach was recently used to collect images of pushed objects and build models of their motions . While the proposed method belongs to the category of model-based RL, it differs from most related methods by the explicit use of the dynamics equations, which drastically improves its data-efficiency. The identified mass distribution can also be used to predict the balance and stability of the object in new configurations that are not covered in the training data. For instance, a GP or a neural net cannot predict if an object remains stable when pushed to the edge of a support table, unless such an example is included in the training data, with the risk dropping the object and losing it.

The mechanics of pushing was explored in several past works , from both a theoretical and algorithmic point of view. Notably, Mason derived the voting theorem to predict the rotation and translation of an object pushed by a point contact. A strategy for stable pushing when objects remain in contact with an end-effector was also proposed in . Yoshikawa and Kurisu proposed a regression method for identifying the support points of a pushed object by dividing the support surface into a grid and approximating the measured frictional force and torque as the sum of unknown frictional forces applied on the grid’s cells. A similar setup was considered in with a constraint to ensure positive friction coefficients. The limit surface is a convex set of all friction forces and torques that can be applied on an object in quasi-static pushing. The limit surface is often approximated as an ellipsoid , or a higher-order convex polynomial . An ellipsoid approximation was also used to simulate the motion of a pushed object to perform a push-grasp . In contrast with our method, these works identify only the friction parameters, and assume that the mass distribution is known or irrelevant in a quasi-static regime.

There has been a recent surge of interest in developing natively differentiable physics engines . A combination of a learned and a differentiable simulator was used to predict effects of actions on planar objects , and to learn fluid parameters . Differentiable physics simulations were also used for manipulation planning and tool use . Recently, it has been observed that a standard physical simulation, formulated as a Linear Complementary Problem (LCP), is also differentiable and can be implemented in PyTorch . In , a differentiable contact model was used to allow for optimization of several locomotion and manipulation tasks.

III Problem Setup and Notation

We consider the problem of displacing a rigid object on a flat homogeneous surface from an initial pose $x_{0}$ to a desired final pose $x^{d}_{T}$. The object has an unknown shape and material properties. We assume that a depth-sensing camera provides a partial 3D view of the object. The partial view contains only the object’s upper surface. A 3D shape is then automatically constructed by assuming that the occluded bottom side is flat. Most of the objects used in our experiments are not flat. However, the autonomously learned friction model simply assigns near-zero friction coefficients to the regions of the bottom surface that do not actually touch the tabletop. Thus, learned near-zero friction forces compensate for wrongly presumed flat regions in the occluded bottom part of the object.

We approximate the object as a finite set of small cuboids. The object is divided into large number of connected cells $1,2,\dots,n$, using a regular grid structure. Each cell $i$ has its own local mass and coefficient of friction that can be different from the other cells. The object’s pose $x_{t}$ at time $t\in[0,T]$ is a vector in $[SE(2)]^{n}$ corresponding to the translation and rotation in the plane for each of the $n$ cells. In other terms, $x_{t}=[p^{1}_{x,t},p^{1}_{y,t},\theta^{1}_{t},\dots,p^{n}_{x,t},p^{n}_{y,t},\theta^{n}_{t}]^{T}$, where $(p^{i}_{x,t},p^{i}_{y,t})$ is the $i^{th}$ cell’s 2D position on the surface, and $\theta^{i}_{t}$ is its angle of rotation. Similarly, we denote the object’s generalized velocity (a twist) at time $t$ by $\dot{x}_{t}=[\dot{p}^{i}_{x,t},\dot{p}^{i}_{y,t},\dot{\theta}^{i}_{t}]_{i=0}^{n}$, where $(\dot{p}^{i}_{x,t},\dot{p}^{i}_{y,t})$ is the $i^{th}$ cell’s linear velocity on the surface, and $\dot{\theta}^{i}_{t}$ is its angular velocity. The object’s mass matrix $\mathcal{M}$ is a diagonal $3n\times 3n$ matrix, where the diagonal is $[\mathcal{I}_{1},\mathcal{M}_{1},\mathcal{M}_{1},\mathcal{I}_{2},\mathcal{M}_{2},\mathcal{M}_{2},\dots,\mathcal{I}_{n},\mathcal{M}_{n},\mathcal{M}_{n}]$, $\mathcal{I}_{i}$ is the moment of inertia of the $i^{th}$ cell of the object, and $\mathcal{M}_{i}$ is its mass. $\mathcal{I}_{i}=\frac{1}{6}\mathcal{M}_{i}w^{2}$ where $w$ is the width of a cuboid. $\mu$ is a $3n\times 3n$ diagonal matrix, where the diagonal is $[\mu_{1},\mu_{1},\mu_{1},\mu_{2},\mu_{2},\mu_{2},\dots,\mu_{n},\mu_{n},\mu_{n}]$. $\mu_{i}$ is the coefficient of friction between the $i^{th}$ cell of the object and the support surface. We assume that: $\forall i\in\{1,\dots,n\}:\mu_{i}\in[0,\mu_{max}]$, where $\mu_{max}$ is a given upper bound. An external generalized force (a wrench) denoted by $F$ is an $1\times 3n$ vector $[f^{1}_{x},f^{1}_{y},\tau^{1},\dots,f^{n}_{x},f^{n}_{y},\tau^{n}]^{T}$, where $[f^{i}_{x},f^{i}_{y}]$ and $\tau^{n}$ are respectively the force and torque applied on cell $i$. External forces are generated from the contact between the object and a fingertip of the robotic hand used to push the object. We assume that at any given time $t$, at most one cell of the object is in contact with the fingertip. Therefore, $F=[0,0,0,\dots,f^{c(t)}_{x},f^{c(t)}_{y},\tau^{c(t)},\dots,0,0,0]^{T}$ where $c(t)\in\{0,\dots,n\}$ is the index of the contacted cell at time $t$.

A ground-truth trajectory $\mathcal{T}^{g}$ is a state-action sequence $(x^{g}_{0},\dot{x}^{g}_{0},F_{0},\dots,x^{g}_{T-1},\dot{x}^{g}_{T-1},F_{T-1},x^{g}_{T},\dot{x}^{g}_{T})$, wherein $(x^{g}_{t},\dot{x}^{g}_{t})$ is the observed pose and velocity of the pushed object, and $F_{t}$ is the external force applied at time $t$, as defined above. A corresponding simulated trajectory $\mathcal{T}$ is obtained by starting at the same initial state $\hat{x}_{0}$ in the corresponding real trajectory, i.e., $\hat{x}_{0}=x_{0}$, and applying the same control sequence $(F_{0},F_{1},\dots,F_{T-1})$. Thus, the simulated trajectory $\mathcal{T}$ results in a state-action sequence $(x^{g}_{0},\dot{x}^{g}_{0},F_{0},x_{1},\dot{x}_{1},F_{1},\dots,x_{T-1},\dot{x}_{T-1},F_{T-1},x_{T},\dot{x}_{T})$, where $x_{t+1}=x_{t}+\dot{x}_{t}dt$ is the predicted next pose. Velocity $\dot{x}_{t}$ is a vector corresponding to translation and angular velocities in the plane for each of the $n$ cells, it is predicted in simulation as $\dot{x}_{t+1}=V(x_{t},\dot{x}_{t},F_{t},\mathcal{M},\mu)$. The goal is to identify mass distribution $\mathcal{M}$ and friction map $\mu$ that result in simulated trajectories that are as close as possible to the real observed ones. Therefore, the objective is to solve the following optimization problem,

Since $x_{t}$ is a vector containing all cells’ positions, the loss is the sum of distances between each cell’s ground-truth pose and its predicted pose, which is equivalent to the average distance (ADD) metric as proposed in . In the following, we explain how velocity function $V$ is computed.

IV Forward Simulation

We adopt here the formulation presented in . We adapt and customize the formulation to exploit the proposed grid-structure representation, and we extend it to include frictional forces between a pushed object and a support surface. The transition function is given as $\dot{x}_{t+1}=x_{t}+\dot{x}_{t}dt$ where $dt$ is the duration of a constant short time-step. Velocity $\dot{x}_{t}$ is a function of force $F_{t}$ and mechanical parameters $\mathcal{M}$ and $\mu$. To find $\dot{x}_{t+1}$, we solve the system of equations of motion that we present in Figure 2, where $x_{t}$ and $\dot{x}_{t}$ are inputs, $[\rho,\xi]$ are slack variables, $[\dot{x}_{t+dt},\lambda_{e},\lambda_{f},\gamma]$ are unknown vectors, and $[\mathcal{M},\mu]$ are hypothesized mass and friction matrices. $\operatorname{diag}(\mathcal{M})$ is a $1\times 3n$ vector corresponding to the main diagonal of $\mathcal{M}$.

$\mathcal{J}_{e}(x_{t})$ is a global Jacobian matrix of all the adjacency constraints in the grid structure. These constraints ensure that the different cells of the object move together with the same velocity. $\mathcal{J}_{e}(x_{t})$ is an $m\times n$ matrix where $n$ is the number of cells, and $m$ is the number of pairs of adjacent cells.

If cell $i$, whose four sides have length $l$, is one of the two adjacent cells in the pair indexed by $k$, then

$\lambda_{e}$ is a $2m\times 1$ variable vector that is multiplied by the Jacobian $\mathcal{J}_{e}(x_{t})$ to generate the vector of impulses $\mathcal{J}^{T}_{e}(x_{t})\lambda_{e}$, which are time-integrals of internal forces that preserve the rigid structure of the object.

$\mathcal{J}_{f}(\dot{x}_{t})$ is an $n\times n$ Jacobian matrix related to the frictional forces between the object’s cells and the support surface, and the corresponding constraints. The main block-diagonal of $\mathcal{J}_{f}(\dot{x}_{t})$ is $[\mathcal{J}^{1}_{f}(\dot{x}_{t}),\mathcal{J}^{2}_{f}(\dot{x}_{t}),\dots,\mathcal{J}^{n}_{f}(\dot{x}_{t})]$, wherein

and the remaining entries of $\mathcal{J}_{f}(\dot{x}_{t})$ are all zeros. $\lambda_{f}$ is a $2n\times 1$ variable vector. $\mathcal{J}_{f}(\dot{x}_{t})$ is multiplied by $\lambda_{f}$ to generate a vector of the frictional forces and torques between the support surface and each cell of the object. $\mathcal{J}_{f}(\dot{x}_{t})$ defines the direction of the frictional forces and torques as the opposite of its current velocity $\dot{x}_{t}=[\dot{\theta}^{i}_{t},\dot{p}^{i}_{x,t},\dot{p}^{i}_{y,t}]_{i=1}^{n}$, whereas $\lambda_{f}$ defines the scalar magnitudes of the frictional forces and torques.

The friction terms have complementary constraints, stated in Fig. 2. These constraints are used to distinguish between the cases when the object is moving and friction magnitudes $\lambda_{f}$ are equal to $\mu\operatorname{diag}(\mathcal{M})$, and the case when the object is stationary and the friction magnitudes $\lambda_{f}$ are smaller than $\mu\operatorname{diag}(\mathcal{M})$. When the object moves, and assuming that the change in the direction of motion happens smoothly, we have $\mathcal{J}_{f}(\dot{x}_{t})\dot{x}_{t+dt}<0$. Therefore, $\gamma>0$ because of the constraints $\rho=\mathcal{J}_{f}(\dot{x}_{t})\dot{x}_{t+dt}+\gamma I$ and $\rho\geq 0$ and $\gamma\geq 0$. Then $\xi=0$ because of the constraint $\gamma\xi=0$. We conclude that $\lambda_{f}=\mu\operatorname{diag}(\mathcal{M})$ from the constraint $\xi+\lambda_{f}=\mu\operatorname{diag}(\mathcal{M})$. Similarly, one can show that $\lambda_{f}<\mu\operatorname{diag}(\mathcal{M})$ if $\dot{x}_{t+dt}=0$.

To simulate a trajectory $(x_{0},\dot{x}_{0},F_{0},x_{1},\dot{x}_{1},F_{1},\dots)$, we iteratively find velocities $\dot{x}_{t+dt}$ by solving the equations in Fig. 2 where $({x}_{t},\dot{x}_{t},F_{t},\mu,\mathcal{M})$ are fixed inputs and the remaining variables are unknown. The solution is obtained, after an initialization step, by iteratively minimizing the residuals from the equations in Fig. 2, using the convex optimizer of .

V Mass and Friction Gradients

To obtain material parameters $[\mathcal{M},\mu]$, a gradient descent on the loss function in Equation 1 is performed. A first approach to compute the gradient is to use the Autograd library for automatic derivation in Python. We propose here a second simpler and faster approach based on deriving analytically the closed forms of the gradients $\frac{\partial loss(\mathcal{M},\mu)}{\partial\mu}$ and $\frac{\partial loss(\mathcal{M},\mu)}{\partial\mathcal{M}}$.

Let us denote by $(\dot{x}^{*}_{t+1},\lambda^{*}_{e},\lambda^{*}_{f},\gamma^{*})$ the solutions for $(\dot{x}_{t+1},\lambda_{e},\lambda_{f},\gamma)$ in the system in Fig. 2. Let us also use $D(x)$ to denote a matrix that contains a vector $x$ as a main diagonal and zeros elsewhere. Finally, let $\operatorname{diag}(\mathcal{M})$ refer to the main diagonal of $\mathcal{M}$. In other terms, $\operatorname{diag}(\mathcal{M})=[\mathcal{I}_{1},\mathcal{M}_{1},\mathcal{M}_{1},\dots,\mathcal{I}_{n},\mathcal{M}_{n},\mathcal{M}_{n}]$. Then,

The differentials of the system are given as

wherein $\partial x_{t}$, $\partial\dot{x}_{t}$, $\partial\mathcal{J}^{T}_{e}(x_{t})$, $\partial\mathcal{J}^{T}_{f}(\dot{x}_{t})$ are all zero matrices and vectors because ${x}_{t}$ and $\dot{x}_{t}$ are fixed and treated as a constant since they are set to ${x}^{g}_{t}$ and $\dot{x}^{g}_{t}$ in Equation 1. Also, $\partial F_{t}=\mathbf{0}$ because the applied force at time $t$ is given as a constant in the identification phase. The differentials can be arranged in the following matrix form: $GX=Y$,

wherein $[\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}]$ is defined as $[\frac{\partial loss}{\partial(-\dot{x}_{t+1})},\mathbf{0},\mathbf{0},\mathbf{0}]G^{-1}$. We use the blockwise matrix inversion to compute $G^{-1}$,

where $A$ and $D$ are square matrices, and $D$ and $(A-BD^{-1}C)$ are invertible. Notice that to compute $[\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}]$, we only need the upper quarter of $G^{-1}$, because the remaining raws will be multiplied by $\mathbf{0}$. Consequently, terms $g(A,B,C,D)$ and $h(A,B,C,D)$ do not matter here, and we only need the terms $(A-BD^{-1}C)^{-1}$ and $-(A-BD^{-1}C)^{-1}BD^{-1}$. The first term $(A-BD^{-1}C)^{-1}$ corresponds to

In the model identification phase, we only utilize data points where the object actually moves when pushed by the robot. Thus, $-\mathcal{J}_{f}(\dot{x}_{t})\dot{x}^{*}_{t+1}-\gamma^{*}=\mathbf{0}$ and $\lambda^{*}_{f}-\mu\operatorname{diag}(\mathcal{M})=\mathbf{0}$, and

Using the blockwise matrix inversion, we find that

We will see in the following that the remaining matrices, $X_{1,2},X_{2,1}$ and $X_{2,2}$, will not be needed. Similarly, the top right of $G^{-1}$ is $-(A-BD^{-1}C)^{-1}BD^{-1}$. It is given as

The first term in $[\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}]$ is then

Since $\frac{\partial loss}{\partial(-\dot{x}_{t+1})}(-\partial\dot{x}_{t+1})=[\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}][\partial\mathcal{M}(\dot{x}^{*}_{t+1}-\dot{x}_{t}),\mathbf{0},\mathbf{0},D(\gamma^{*})\partial(\mu\operatorname{diag}(\mathcal{M}))]^{T}$, then

In the following, we show how to use the equation above to derive $\frac{\partial loss}{\partial\mathcal{M}}$ and $\frac{\partial loss}{\partial\mu}$ and use them in a coordinate descent algorithm to identify $(\mathcal{M}^{*},\mu^{*})$ from data.

We calculate $\frac{\partial loss}{\partial\mathcal{M}}$ while setting $\partial\mu=\mathbf{0}$.

From the definition of the loss function in Equation 1, we can see that \frac{\partial loss}{\partial(\dot{x}_{t+1})}=2dt\sum_{t=1}^{T-1}D\big{(}x_{t+1}-x^{g}_{t+1}\big{)}, wherein $x_{t+1}^{g}$ is the observed ground-truth pose of the object and $x_{t+1}$ is its predicted pose, computed as $x_{t+1}=x^{g}_{t}+\dot{x}^{*}_{t}dt$. Finally,

V-B Friction Gradient

We calculate $\frac{\partial loss}{\partial\mu}$ while setting $\partial\mathcal{M}=\mathbf{0}$.

V-C Mass and Friction Identification Algorithm

VI Policy Gradient

After identifying the mass distribution $\mathcal{M}$ and friction map $\mu$ using Algorithm 1, we search for a new sequence of forces $(F_{t})_{t=0}^{T-1}$ to push the object toward a desired terminal goal configuration $x_{T}^{d}$. Algorithm 2 summarizes the main steps of this process. We start by creating a rapidly exploring random tree ($RRT^{*}$) to find the shortest path from $x_{0}$ to $x_{T}^{d}$. While searching for the shortest path, we eliminate from the tree object poses that are unstable (based on the identified mass distribution $\mathcal{M}$) or that are in collision with other objects. $RRT^{*}$ returns a set of waypoints $\mathcal{X}_{waypoint}$. At each iteration of the main loop of the algorithm, we find the nearest waypoint in $\mathcal{X}_{waypoint}$ and search for actions that would push the object toward it. A pushing force is parameterized by a contact point, a direction and a magnitude, as discussed in Section III. We focus here on optimizing the contact point, and we keep the magnitude constant. The direction of the force is chosen to be always horizontal. It is given as the opposite of the surface normal of the object at the contact point, projected down on the 2D plane of the support surface. This choice is made to avoid slippages and changes in contact points during a push.

A contact point is always located on the outer side of a cuboid (cell). Therefore, we limit the search to the outer cells of the grid. The objective of this search is to select a contact point that reduces the gap between the predicted pose of the object after pushing it, and the nearest waypoint $x_{target}$ that has not been reached yet. We select the initial contact point as the outer cell that is most aligned to the axis $\hat{x}-x_{target}$, where $\hat{x}$ is the estimated center of mass. The gradient of the gap with respect to the contact point is computed by using the finite-difference method. The contact point is moved in the direction that minimizes the gap until a local optimum is reached. Force $F_{t}$ is then defined based on the selected contact point. The pose and velocity of the object are replaced by the predicted ones that result from applying force $F_{t}$. This process is repeated until the object reaches the desired goal configuration. The time duration of each pushing action is also optimized by using finite differences. This part is omitted for simplicity’s sake.

The surface of the object often contains non-differentiable parts where the analytical gradient with respect to the contact point is undefined. Even on smooth parts, there is no clear advantage of computing the analytical gradient here, because the space of contact points is uni-dimensional, in contrast to the high-dimensional space of non-uniform mass and friction distributions. In low-dimensional search spaces, finite-difference methods are computationally efficient.

VII Main Algorithm

Algorithm 3 summarizes the main steps of the proposed approach and the protocol followed in the experiments. In summary, the robot first “plays” with the unknown object by applying random short-lasting horizontal forces for a safe and local exploration. A mass and friction model is then inferred from the gathered data by using Algorithm 1. Based on the inferred mass distribution, a safe goal configuration is sampled from a desired goal region. For example, a pregrasp sliding manipulation can be used to grasp a thin object that cannot be directly grasped from a flat surface. In , a known object is pushed to the edge of a table and then grasped from there. Pushing an unknown object to the edge of a table results often in losing the object. Our method avoids this issue by sampling a goal configuration that allows a sufficient part of the object to be graspable, while keeping the object balanced on the edge thanks to the identified mass distribution. Once a goal is selected, Algorithm 2 is used to generate a sequence of actions to push the object to the goal.

VIII Evaluation

We report here the results of three sets of experiments to evaluate the proposed approach. The most important set is the one related to mass and friction identification.

The experiments are performed on both simulated and real robot and objects. In the real robot setup, a Robotiq 3-finger hand mounted on a Kuka robot is repeatedly moved to collide with a rigid object that is set on a table-top and to push it forward, as shown as Figure 1. The initial, final and intermediate point clouds of the object are recorded. The simulation experiments are performed using the physics engine Bullet and models of the robot and objects. The experiments are performed on five real objects: a book, a hammer, a snack, a toolbox and a spray gun, and eight simulated objects: a box, a hammer, a book, a crimp, a snack, a ranch, a spray gun and a toothpaste. The number of cells per object varies from $28$ to $88$ depending on the size of the object.

VIII-B Tasks

Model Identification. Each real and simulated object is pushed by the robot randomly $10$ times on the table. Half of the recorded trajectories are used for learning a mass matrix $\mathcal{M}$ and friction map $\mu$. The other half is used for testing the identified models. Since the ground-truth values of mass and friction are unknown, the identified models are evaluated in terms of the accuracy of the predicted pose of each cell after applying the sequence of actions provided in the test set. The experiments on the real objects are repeated with $25$ randomized splits into training and testing sets. The simulation experiments are repeated with $10$ different models per object.

Planning and Control. For each one of the eight objects in simulation, we randomly sample $10$ values for their mass and friction matrices, and $10$ random goal configurations in a disk of a radius of $1m$ around the initial configuration. The rotations of the goal configurations are also selected randomly. The number of settings is then $8\times 10\times 10$. The task is to generate, in each setting, a sequence of forces that pushes the object from the initial configuration to the goal. The objective is to asses the computational efficiency of Algorithm 2.

Pre-grasp Sliding Manipulation. Finally, we evaluate the entire system (Algorithm 3) on the task of sliding an object from a random initial pose on a table to a desired goal region at the edge of the table where the object can be grasped. The object cannot be directly grasped from a flat surface, the goal is to push it to the edge where part of it sticks out of the table and becomes graspable. Since the mass distribution of the object is highly heterogeneous and unknown, the object often becomes unbalanced at the edge and falls from the table if the identified model is incorrect. We report here the percentage of experiments where the object is successfully pushed to the edge and grasped without losing it. The experiments on this task are performed using the real Kuka robot and a real hammer. The exploration phase contains only $5$ random pushing actions that are used for model identification. The reported results are averaged over $16$ independent runs, with a different initial pose in each run.

VIII-C Compared Methods

Model Identification. Algorithm 1 is compared against the following methods. Random search is a baseline method that repeatedly samples random values of the mass and friction matrices and returns the best sampled values that minimize $loss(\mathcal{M},\mu)$. Weighted sampling search generates random values uniformly in the first iteration, and then iteratively generates normally distributed random values around the best parameter obtained in the previous iteration. The standard deviation of the random values is gradually reduced over time, to focus the search on the most promising region of the search space. The finite differences gradient is an approximation of the analytical gradient. We add or subtract a small amount to the current parameter values and simulate the trajectories using the neighboring parameter values to approximate the derivatives $\frac{\partial loss}{\partial\mathcal{M}}$ and $\frac{\partial loss}{\partial\mu}$. Because a large number of simulations is required to compute the gradient for all the cells, the parameters of each cell are updated using the coordinate gradient descent. We also compare the proposed method with two black-box optimization methods: CMA-ES and Nelder-Mead, and an automatic differentiation of the LCP solver using the Autograd function of PyTorch . The same minimum and maximum bounds of mass and friction are provided to all methods and are also used for all cells of objects.

Planning and Control. We perform an ablation study where we substitute the finite-difference gradient in Algorithm 2 with an exhaustive search of the optimal contact point.

Pre-grasp Sliding Manipulation. We compare Algorithm 3 to two alternatives. The first one assumes a uniform and homogenous mass and friction values, and uses directly Algorithm 2 to push the object to the goal without model identification. The second alternative is identical to Algorithm 3, except that no upper limit on the friction coefficients is used.

VIII-D Results

Model Identification. Figure 3 shows mass distributions identified using Algorithm 1. The mass distributions of the book, the hammer, the snack, the toolbox and the spray gun are all identified using the real objects and robot. The experiments on the remaining objects were performed in simulation because they were too thin for a safe robotic manipulation. The results show that the identified models quickly converge to the ground-truth models. Figure 4 (a)-(m) shows the average distance between the predicted cell positions and the ground-truth ones in the test data as a function of the number forward physics simulations used by the different identification methods. In our method, the number of physics simulations is the same as the number of steps of the gradient-descent, since one simulation is performed after each update of $\mathcal{M}$ and $\mu$. Note that the physics simulations dominate the computation time, with $1.3448(\pm 0.5604)$ second per a simulation, while the computation of the gradient using the proposed algorithm takes only $0.0069(\pm 0.0024)$ second. The results demonstrate that global optimization methods suffer from the curse of dimensionality due to the combinatorial explosion in the number of possible parameters for all cells. The results also demonstrate that the proposed method can estimate the parameters within a surprisingly small number of gradient-descent steps and a short computation time ($30$ seconds). The average error of the predicted cell positions using the identified mass and friction distributions is less than $1.53cm$ in simulation and $2.27cm$ in the real experiments. Note that the error in real experiments is higher due to sensing and control errors. The finite differences approach also failed to converge to an accurate model due to the high computational cost of the gradient computation, as well as the sensitivity of the computed gradients to the choice the grid size. Figure 4 (o) shows that the number of training actions improves the accuracy of the model learned by Algorithm 1. Increasing the number of training actions allows the robot to uncover properties of different parts of the object more accurately.

Planning and Control. Figure 4 (p) shows the number of physics simulations used to optimize the contact location. The proposed policy gradient algorithm requires a small fraction of the exhaustive search’s computation time, while both methods attained a 100$\%$ success rate in finding a sequence of pushes that reaches the goal in simulation.

Pre-grasp Sliding Manipulation. This task integrates the two previous tasks and evaluates the entire proposed system. Figure 4 (q) shows that in all of the $16$ trials with random initial poses, the robot successfully identified the hammer’s mass and friction distributions using Algorithm 1, selected a physically stable goal configuration in the desired goal region based on the identified model, planned and executed a sequence of pushing actions using Algorithm 2, and grasped the object from the part pushed out of the table. Only $12$ trials resulted in successful grasps when the mass distribution was not learned and was assumed to be uniform. The upper bound on the dynamic coefficient of friction is set to $1.00$ in Algorithm 1. The success rate drops to $13/16$ when no upper bound limit is set on the friction in Algorithm 1.

IX Final Remark

The decline in the success rate in Figure 4 (q) when no upper limit on the friction is used is an important observation. The algorithm attributed most of the rotations to unrealistically high frictions in certain regions, instead of an uneven mass distribution. Despite identifying a wrong model of mass and friction, the predicted motion of the object was accurate as long as the object was entirely on the table. But once pushed to the edge, and its heavy side is not anymore supported by the table’s surface, the object falls. This clearly demonstrates the importance of identifying not only the friction, but also the accurate mass distribution of objects for a safe manipulation. The upper bound can be seen as an inductive bias that is necessary for learning from limited data.

References