The performance and properties of (asymptotically) optimal grasp planning algorithms heavily depend on the type of grasp quality metrics. A summary of these metrics can be found in . Usual requirements for a good grasp include: force closure, small contact force magnitudes, the preference of normal forces over frictional forces, higher resilience to external wrenches. The type of a metric not only reflects the requirements of an application, but also incurs different planning algorithms. For example, the metrics are submodular, which allows fast discrete grasp point selection . The metric has an optimizable lower-bound, which allows an optimization-based grasp planning algorithm  to jointly search for grasp points and grasp poses.
However, all the metrics considered so far take a common assumption about the target object: these objects are of infinite stiffness and will never be broken. As a result, we can assume that the target object is a rigid body and all the forces and torques are applied on the center-of-mass, which greatly simplify the computation and analysis of metrics. However, this assumption does not hold when grasping fragile objects where certain weak parts of an object should not be touched to avoid possible fractures. The grasp of fragile objects have been considered in prior works [16, 1], which attempt to avoid breaking objects by developing safer grippers. Instead, we argue that, in addition to better robot hardware, a new metric is needed to guide the grasp planning algorithm, so that we can find a grasp that is most unlikely to break the object.
Main Results: We present a novel metric that relaxes the infinite stiffness assumption and takes the material of the target object into consideration. Based on the theory of brittle fracture , we formulate the set of external wrenches that can be resisted without causing fractures in the object. Similar to the metric [7, 19], our metric then measures the size of the space of resistible wrenches. We call this new metric stress-minimizing (SM) metric . We show that, using boundary element method (BEM) , can be computed efficiently, given a closed surface triangle mesh of the object and a set of contact points. The cost of computing is on par with that of computing . In addition, we show that conventional optimization-based grasp planning algorithms  can be modified to search for globally optimal grasps maximizing .
The rest of the paper is organized as follows. We review related work in sec:results and formulation the problem of grasp planning in sec:problem. The formulation of and its properties are summarized in sec:metric. Grasp planning algorithms using as metric are formulated in sec:planning. Finally, we compare with conventional metric in sec:results.
Ii Related Work
We review related work in grasp quality metric, material and fracture modeling, and grasp planning.
Grasp Quality Metric: Although some grasp planners only consider external wrenches along some certain directions , more prominent characterization of robust grasp requires resilience to external wrenches along all directions, which is known as force closure . However, infinitely many grasps can have force closure, of which good grasps are characterized by different quality metrics . Most of the metrics are designed such that implies force closure. A closely related metric to our is , which measures the maximal radius of origin-centered wrench-space circle contained in the convex set of resistible wrenches, under contact force magnitude constraints.
Material and Fracture Modeling: Real world solid objects will undergo either brittle or ductile fractures depending on their materials. But modeling them are both theoretically and computationally difficult . Fortunately, for grasp planning, we do not need to model the deformation of objects after fractures, but only need to detect where fractures might happen. In this case, the theory of linear elasticity suffices, which is efficient to compute using finite element method (FEM)  or boundary element method (BEM) . We use BEM as our computational tool because it only requires a surface triangle mesh which is more amenable to robot grasp applications.
Grasp Planning: Given a grasp quality metric, an (asymptotically) optimal grasp planning algorithm finds a grasp that maximizes the grasp quality metric. Early algorithms  use sampling-based methods for planning. These algorithms are very general and they are agnostic to the type of grasp quality metrics. However, more efficient algorithms such as [6, 19] can be designed if quality metrics have certain properties. Recent work 
uses a set of precomputed quality metrics to train a grasp quality function represented by a deep neural network and then uses the function as reward to optimize a grasp planner via deep reinforcement learning.
Iii Problem Statement
In this section, we formulate the problem of stress-minimizing grasp planning. Throughout the paper, we keep a 3D target object in its reference space, where the origin coincides with its center-of-mass. The object takes up a volume which is a closed subset . In addition, the object is under an external 6D wrench and a set of external contact forces at contact points with unit contact normals . If a grasp is valid, we have the following wrench balance condition:
where is the frictional coefficient. To compare the quality of two different grasps, a well-known method is to compare their metric , which is the maximal radius of the origin-centered inscribed sphere in the convex hull of all possible resistible external wrenches when the magnitude of is bounded. Mathematically, this is:
where is the
positive semi-definite metric tensor in the wrench space.
However, this conventional formulation assumes that the object will never be broken however large the external forces are. To relax this condition, we have to make use of the numerical models of brittle fracture, e.g. . In these formulations, we assume that the object is made of homogeneous isotropic elastic material with being its Lamé material parameters. When under external force fields , an infinitesimal displacement and an stress field will occur . can be computed from using the force balance condition:
where we assume the boundary of is almost everywhere smooth with unit outward normal defined as and is the Dirac’s delta operator. Classical theory of brittle fracture further assumes that there exists a tensile stress and brittle fractures will not happen if the following condition holds:
Propose a grasp quality metric that measures the quality of grasps with eq:safe_condition as precondition (sec:metric).
Propose grasp planning algorithms that generate grasps maximizing (sec:planning).
Iv The Stress-Minimizing Metric
Our construction is illustrated in fig:pipeline. The basic idea behind the construction of is very similar with that of the metric. Intuitively, we first define a convex subset of the 6D wrench space, which contains resistible wrenches that does not violate eq:safe_condition. We then define as the maximal radius of the origin-centered sphere contained in .
Iv-a Definition of and
Given a certain wrench , we need to determine whether . This can be done by first computing and then testing whether eq:safe_condition holds. However, is computed from but not , so that we need to find a relationship between and . In other words, we need to find a body force distribution such that the net effect of is equivalent to applying on the center-of-mass. Obviously, infinitely many will satisfy this relationship and different choices of will lead to different variants of metrics. In this paper, we propose to choose as a linear function in . The most important reason behind this choice is that the computation of can be accomplished using BEM if is a harmonic function of . Under this choice, we have: , where is the constant term, is the constant spatial derivative tensor. Clearly
has 12 degrees of freedom and we can solve forand to equate the effect of and as follows:
where are the first, second, and third column of , respectively. However, there are 9 degrees of freedom in but only 3 constraints, so that we have to solve for in a least square sense:
the solution of which can be computed analytically. In summary, we have:
where denotes Kronecker product. The matrices are constants and can be precomputed from the shape of the target object, or more specifically, from the inertia tensor. Given these definitions, we can now define as follows:
Finally, we are ready to give a mathematical definition of as the following optimization problem:
From the mathematical definition of , we immediately have the following properties of :
is a convex set.
eq:force_balance is a set of quadratic cone constraints, which defines a convex set. eq:elasticity is a set of infinite-dimensional linear constraints, which defines a convex set. eq:safe_condition is an infinite-dimensional PSD-cone constraint, which defines a convex set. Finally, eq:gg_from_w is a linear constraint, which also defines a convex set. As the intersection of convex sets, is convex.
is a compact set so that is finite.
For any , that satisfies eq:force_balance, eq:elasticity, and eq:gg_from_w cannot be uniformed zero. In other words, there exists and such that . If we multiply by , eq:safe_condition will be violated so that . Therefore, is bounded and is obviously closed, so that is compact and is finite.
In addition, the following property of is obvious:
implies force closure.
And the following property has been proved in  for and also holds for by a similar argument:
Iv-B Discretization of
The computation of exact is impossible because it involves infinite dimensional tensor fields: , so that we have to discretize them using conventional techniques such as FEM  or BEM . In comparison, FEM is mathematically simpler but requires a volumetric mesh of the target object, while BEM only requires a surface triangle mesh. We provide the detailed derivation of BEM in appen:BEM and summarize the main results here. Our BEM implementation approximates the stress field to be piecewise constant on each triangular patch of the surface. Assuming that the target object has surface triangles whose centroids are: , we have different stress values:
where are dense coefficient matrices defined from BEM discretization, the definitions of which can be found from appen:BEM. Note that computing the coefficients of these two matrices are very computationally costly, where a naive implementation of BEM requires operations and acceleration techniques such as the H-matrix  can reduce this cost to operations. However, these two matrices are constant and can be precomputed for a given target object shape, so that the the cost of BEM computation is not a part of grasp planning. After discretization, we arrive at the finite-dimensional version of fracture condition:
finite-dimensional version of :
and finite-dimensional version of defined as:
All the properties of the infinite-dimensional and hold for the finite-dimensional version and by a similar argument.
Iv-C Computation of
Even after discretization, computing is costly and non-trivial. According to Lem:opt_formulation, the equivalent optimization problem for is:
which is non-convex optimization, so that direct optimization does not give the global optimum. In this section, we modify three existing algorithms to (approximately) compute . Two discrete algorithms have been proposed in [23, 19] to approximate . Due to the close relationships between and , we show that these two algorithms can be modified to compute .
, the space of unit vectors is discretized into a finite set ofdirections: . As a result, we can compute an upper bound for as:
The global optimum of pb:lower_bound is easy to compute from an H-representation of by computing the distance between the origin and each face of . This lower bound can be iteratively tightened by first computing the blocking face normal of and then expanding via:
These two algorithms can be realized if we can find the supporting point of , which amounts to the following conic programming problem:
The conic programming reformulation in pb:support can be solved using the interior point method . Given this solution procedure, we summarize the modified version of  in Alg:UpperBound and modified version of  in Alg:LowerBound. Note that Alg:LowerBound is advantageous over Alg:UpperBound in that it can approximate up to arbitrary precision , so we always use Alg:LowerBound in the rest of the paper.
Compared with metric, a major limitation of using metric is that the computational cost is much higher. Note that the computational cost of solving pb:support is at least linear in and can be superlinear depending on the type of conic programming solver used. This is the number of surface triangles on the target object, which can easily reach several thousands. Fortunately, we can drastically reduce this cost by using a progressive approach.
Iv-D Performance Optimization
The naive execution of Alg:LowerBound can be prohibitively slow due to the repeated solve of pb:support. The conic programming problem has PSD-cone constraints with being several thousands. Solving pb:support using interior point method  involves repeated solving a sparse linear system with size proportional to . We propose a method that can greatly improve the performance when solving pb:support. Our idea is that when the global optimum of pb:support is reached, more than of the PSD-cone constraints are inactive, so that removing these constraints do not alter the solution. This idea is inspired by  which shows that, empirically, maximal stress only happens on a few sparse points on the surface of the target object. However, we do not know the active constraints as a prior. Therefore, we propose to progressively detect these active constraints.
To do so, we first select a subset such that and are the stresses that are most likely to be violated. To select this set , we use a precomputation step and solve pb:support for times using random , and record which PSD-cones are active. For each PSD-cone, we maintain how many times they become active during the solves of pb:support. We then select the most frequent PSD-cones to form . After selecting , we maintain an active set which initializes to and we solve pb:support using constraints only in , which is denoted by:
After we solve for the global optimum of pb:support_cons, we check the stress of remaining constraint points and we pick the most violated constraint:
If we have , then pb:support_cons and pb:support will return the same solution. Otherwise, we add to . This method is summarized in Alg:progressive and is guaranteed to return the same global optimum of pb:support. In practice, Alg:progressive is orders of magnitude more efficient than throwing all constraints to the interior point method at once.
V Grasp Planning Under The SM Metric
Built on top of the computational procedure of , we can solve the optimal grasp point selection problem. Given a set of potential grasp points sampled on , we want to select points. Mathematically, we want to solve the following mixed-integer optimization problem using branch-and-bound (BB) algorithm :
where is the big-M constant that can be arbitrarily large. Note that the algorithm to solve pb:planning can be different depending on what algorithm we use to (approximately) compute . If approximate Alg:UpperBound is used, then the problem can be efficiently solved using the sub-modular coverage algorithm . However, in order to highlight the advantage of our new metric, we would like to compute accurately using Alg:LowerBound. In sec:opt_planning, we show that pb:planning can be solved alot more efficiently using a series of reformulations, without changing the global optimum.
V-a Optimized BB Algorithm
Our key innovation is via the following reformulation:
where is the special-ordered-set-of-type-1 , which constraints that only one number in a set can take non-zero value. According to , we know that pb:planning requires binary variables while pb:planning_relaxed requires binary variables, which is much fewer as . Intuitively, pb:planning_relaxed build a binary bounding volume hierarchy for the set of contact points and introduce one binary decision variable for each internal level of the tree to select whether the left or the right child is selected. After a leaf node is reached, a contact point is selected. Finally, the last constraint in pb:planning_relaxed, , reflects the order-independence of contact points.
Unfortunately, no optimization tools can solve pb:planning_relaxed in an off-the-shelf manner due to the special objective function, so that we develope a special implementation of BB outlined in Alg:BB. In this algorithm, our binary bounding volume hierarchy is a KD-tree, as illustrated in fig:pipelinede. The most important component for the efficiency of BB algorithm is the problem relaxation (ln:relax) as follows:
which can be solving using Alg:LowerBound by excluding the contact points in . We summarize our main results below:
Alg:BB returns the global optimum of pb:planning_relaxed, which is also the global optimum of pb:planning by setting: .
Finally, note that all the results in this section hold if we replace with .
We implement our algorithms for computing and perform grasp planning using C++. The accuracy of BEM heavily relies on the quality of the surface triangle mesh, so that we first optimize the mesh quality using CGAL . We implement the BEM using kernel independent numerical integration scheme . The most computationally costly step in BEM is the inversion of system matrices, for which we use LU-factorization accelerated by H-matrices . Finally, we use CGAL  to construct convex hulls with exact arithmetics. All the experiments are performed on a single desktop machine with two Xeon E5-2697 CPU and 256Gb memory.
Parameter Choices: Computing requires more parameters than . Specifically, there are three additional variables: tensile stress and Lamé material parameters: . However, if we transform to an equivalent set of parameters: Young’s modulus and Poisson ratio , it is obvious to show that is proportional to and inversely proportional to . Since the absolute value of a grasp metric is meaningless for grasp planning and only the relative value matters, we can always set and choose only according to the material type of the target object, and then setting:
We have for copper and for rubber. In all experiments, we set . Finally, when running Alg:LowerBound, we set .
Shape-Awareness: The most remarkable advantage of over is shape awareness. In fig:shapeAwareness, our target object is a U-shaped tuning fork and we use a 3-point grasp. The shape of the tuning fork is asymmetric along the X-axis, and according to fig:shapeAwarenessabc, is not aware of the asymmetry, while correctly reflects the fact that grasping the leftmost point is better than grasping the mid-left point because it is less likely to break the object. However, the best grasp under and are the same, i.e., grasping the centroid point. If we change the metric and emphasize force resistance over torque resistance, then the difference between and is more advocated.
Robustness to Mesh Resolution: The change of is not sensitive to the resolution of surface meshes as shown in fig:shapeAwarenessabc, which makes robust to target objects discretized using small, low-resolution meshes. As we increase from to and finally to , the change of against is almost intact, with very small fluctuations around .
Computational Cost: does incur a higher computational cost than . The most computational cost lies in the assembly of matrices which involves the direct factorization of a large dense matrix. But this assembly is precomputation and required only once for each target object before grasp planning. In fig:shapeAwarenessabc, this step takes s when , s when , and s when . After precomputation, the cost of evaluating and are very similar, as shown in fig:cost. This implies that using does not incur a higher cost in grasp planning. This is largely due to the progressive Alg:progressive, which greatly reduce the number of constraints in solving pb:support. Without this method, solving pb:support is prohibitively costly by requiring the solve of a sparse linear system of size proportional to .
Grasp Planning: In fig:plan_results, we show globally optimal grasps for 8 different target objects under both the and metric. To generate these results, we choose contact points from potential contact points by running Alg:BB. These contact points are generated using Poisson disk sampling. The computational cost of Alg:BB is hr under and hr under on average. This result is surprising as the cost of computing is comparable to that of computing . We found that tends to create more local minima so that BB needs to create a larger search tree under . In the third row of fig:plan_results, we show the maximal stress configuration in and the corresponding stress configuration under side-by-side. The advantage of is quite clear which suppresses the stress to resist the same external wrench. For some target objects, the high stress is concentrated in a very small region and we indicate them using black circles.
Vii Conclusion and Limitation
We present SM metric, which reflects the tendency to break a target object. As a result, a grasp maximizing
will minimize the probability of breaking a fragile object. We show thatcan be computed using previous methods and its computational costly can be drastically reduced by progressively detecting the active set. Finally, we show that grasp planning under can be performed using BB algorithms. Our experiments show that is aware of geometric fragility while is not. We also show that using does not increase computationally cost in grasp planning.
The major limitation of our work is that computing requires a costly precomputation step for solving the BEM problem. In addition, the BEM problem requires high quality watertight surface meshes of target objects. An avenue of future research is to infer the value of
for given unknown objects using machine learning, as is done in.
-  A. Al-Ibadi, S. Nefti-Meziani, and S. Davis, “Active soft end effectors for efficient grasping and safe handling,” IEEE Access, vol. 6, pp. 23 591–23 601, 2018.
-  P. Alliez, C. Jamin, L. Rineau, S. Tayeb, J. Tournois, and M. Yvinec, “3D mesh generation,” in CGAL User and Reference Manual, 4.14 ed. CGAL Editorial Board, 2019.
E. D. Andersen and K. D. Andersen, “The mosek interior point optimizer for linear programming: an implementation of the homogeneous algorithm,” inHigh performance optimization. Springer, 2000, pp. 197–232.
-  J. Clausen, “Branch and bound algorithms-principles and examples,” 1999.
-  T. A. Cruse, D. Snow, and R. Wilson, “Numerical solutions in axisymmetric elasticity,” Computers & Structures, vol. 7, no. 3, pp. 445–451, 1977.
-  H. Dai, A. Majumdar, and R. Tedrake, Synthesis and Optimization of Force Closure Grasps via Sequential Semidefinite Programming. Cham: Springer International Publishing, 2018, pp. 285–305.
-  C. Ferrari and J. Canny, “Planning optimal grasps,” in Proceedings 1992 IEEE International Conference on Robotics and Automation, May 1992, pp. 2290–2295 vol.3.
-  P. Fiala and P. Rucz, “Nihu: An open source c++ bem library.” Advances in Engineering Software, vol. 75, pp. 101–112, 2014.
-  G. A. Francfort and J.-J. Marigo, “Revisiting brittle fracture as an energy minimization problem,” Journal of the Mechanics and Physics of Solids, vol. 46, no. 8, pp. 1319–1342, 1998.
-  B. Grünbaum and G. C. Shephard, “Convex polytopes,” Bulletin of the London Mathematical Society, vol. 1, no. 3, pp. 257–300, 1969.
-  W. Hackbusch, “A sparse matrix arithmetic based on h-matrices. part i: Introduction to h-matrices,” Computing, vol. 62, no. 2, pp. 89–108, 1999.
-  S. Hert and M. Seel, “dD convex hulls and delaunay triangulations,” in CGAL User and Reference Manual, 4.14 ed. CGAL Editorial Board, 2019.
-  T. J. Hughes, The finite element method: linear static and dynamic finite element analysis. Courier Corporation, 2012.
J. Mahler, J. Liang, S. Niyaz, M. Laskey, R. Doan, X. Liu, J. Aparicio, and K. Goldberg, “Dex-net 2.0: Deep learning to plan robust grasps with synthetic point clouds and analytic grasp metrics,” 07 2017.
-  J. Mahler, F. T. Pokorny, Z. McCarthy, A. F. van der Stappen, and K. Goldberg, “Energy-bounded caging: Formal definition and 2-d energy lower bound algorithm based on weighted alpha shapes,” IEEE Robotics and Automation Letters, vol. 1, no. 1, pp. 508–515, 2016.
-  B. Matulevich, G. E. Loeb, and J. A. Fishel, “Utility of contact detection reflexes in prosthetic hand control,” in 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, Nov 2013, pp. 4741–4746.
-  V. . Nguyen, “Constructing force-closure grasps,” in Proceedings. 1986 IEEE International Conference on Robotics and Automation, vol. 3, April 1986, pp. 1368–1373.
-  M. A. Roa and R. Suárez, “Grasp quality measures: review and performance,” Autonomous robots, vol. 38, no. 1, pp. 65–88, 2015.
-  J. D. Schulman, K. Goldberg, and P. Abbeel, “Grasping and fixturing as submodular coverage problems,” in Robotics Research. Springer, 2017, pp. 571–583.
-  O. Steinbach, Numerical approximation methods for elliptic boundary value problems: finite and boundary elements. Springer Science & Business Media, 2007.
-  J. P. Vielma and G. L. Nemhauser, “Modeling disjunctive constraints with a logarithmic number of binary variables and constraints,” Mathematical Programming, vol. 128, no. 1-2, pp. 49–72, 2011.
Z. Zhang, J. Gu, and J. Luo, “Evaluation of genetic algorithm on grasp planning optimization for 3d object: A comparison with simulated annealing algorithm,” in2013 IEEE International Symposium on Industrial Electronics, May 2013, pp. 1–8.
-  Y. Zheng, “An efficient algorithm for a grasp quality measure,” IEEE Transactions on Robotics, vol. 29, no. 2, pp. 579–585, 2012.
-  Q. Zhou, J. Panetta, and D. Zorin, “Worst-case structural analysis,” ACM Trans. Graph., vol. 32, no. 4, pp. 137:1–137:12, July 2013.