A Convex Polynomial Force-Motion Model for Planar Sliding: Identification and Application  Jiaji Zhou, Robert Paolini, J. Andrew Bagnell and Matthew T. Mason The Robotics Institute, Carnegie Mellon University  { jiajiz, rpaolini, dbagnell, matt.mason } @cs.cmu.edu  Abstract — We   propose   a   polynomial   force-motion   model for   planar   sliding.   The   set   of   generalized   friction   loads   is the 1-sublevel set of a polynomial whose gradient directions correspond to generalized velocities. Additionally, the polyno- mial is confined to be convex even-degree homogeneous in order to obey the maximum work inequality, symmetry, shape invariance in scale, and fast invertibility. We present a simple and statistically-efficient model identification procedure using a sum-of-squares convex relaxation. Simulation and robotic ex- periments validate the accuracy and efficiency of our approach. We also show practical applications of our model including stable pushing of objects and free sliding dynamic simulations.  I. INTRODUCTION We develop a data-driven but physics-based method for modeling planar friction. Manipulations employing friction are ubiquitous in tasks including positioning and orienting objects by pushing [1]–[4], controlled slip with dexterous hands [5] and assembly of tight-fitting parts [6]. In the case of planar robot pushing, indeterminacy of the pressure distribution between the object and support surface leads to uncertainty in the resultant velocity given a particular push   action.   Despite   such   inherent   difficulty,   algorithms and analysis have been developed with provable guarantees. Mason [1] derived the voting theorem to determine the sense of rotation of an object pushed by a point contact. Lynch and Mason [2] developed a stable pushing strategy when objects remain fixed to the end effector with two or more contact points. However, minimal assumptions on friction conditions inherently lead to conservative strategies. By explicitly mod- eling and identifying the friction space, we can improve strategies for planning and control. Our contribution lies in developing a precise and statistically-efficient (i.e., requiring only a few collected force-velocity data pairs) model with a computationally efficient identification procedure. Fig. 1 illustrates an outline of the paper. We assume a quasi-static regime [7] where forces and moments are balanced with negligible inertia effects. II. B ACKGROUND ON   P LANAR   F ORCE -M OTION   M ODELS  The classical Coulomb friction law states that for a point contact with instantaneous planar velocity   v   = [ v x ,   v y ] T   , the incurred friction force   f   = [   f x ,   f y ] T   the point applies on the surface is parallel to   v , i.e.,   f / | f |   =   v / | v | . We refer the readers to [1] for details of friction analysis for planar sliding under isotropic Coulomb friction law. In this paper, we build our analysis on a generalized friction law formulated first in [8],  Fig. 1:   The robot randomly pokes the object of known shape with a point finger to collect force-motion data. We then optimize a convex polynomial friction representation with physics-based constraints. Based on the representation, we demonstrate applications of stable pushing and dynamic sliding simulation.  in which   v   and   f   may not be parallel, but only need to obey the maximum work inequality:  ( f   −   f ′ )   ·   v   ≥   0 ,   (1) where   f ′   is an arbitrary element from the set of all possible static and kinetic friction forces. Let   V   = [ V x , V y ,   ω ] T   be the instantaneous generalized ve- locity and   F   = [ F x ,   F y ,   τ ] T   be the generalized friction load for a rigid body sliding on a planar surface with a contact area  R . Both   V   and   F   are in the local body frame 1 .   F   can be computed by integration over   R :  F x   =  ∫  R  f ax   da ,   F y   =  ∫  R  f ay   da ,   τ   =  ∫  R  ( r ax   f ay   −   r ay   f ax )   da .  (2) The maximum work inequality in equation (1) can be ex- tended for generalized friction load   F   and velocity   V :  F   ·   V   =  ∫  R  f ax ( V x   −   ω r ay )   da   +  ∫  R  f ay ( V y   +   ω r ax )   da  =  ∫  R  f ax v ax   +   f ay v ay   da   =  ∫  R  f a   ·   v a   da   ≥   F ′   ·   V ,   (3)  1 Throughout the paper, we use a local coordinate frame with the origin set as the projection of the COM onto the supporting surface. However, the choice of the origin can be any other point of convenience.  arXiv:1602.06056v3 [cs.RO] 16 Jun 2016
among   any   other   possible   friction   load   F ′ .   Due   to   the converse supporting hyperplane theorem [9], the set of all generalized friction loads form a convex set   F   . An important work that inspires us is Goyal et al. [10] who found that all possible generalized friction loads during sliding form a limit surface (LS) constructed from the Minkowsky sum of limit curves at individual support points. Points inside the surface correspond to static friction loads. Points on the surface correspond to friction loads with normals parallel to sliding velocity directions, forming a mapping between generalized friction load and sliding velocity. An ideal LS is always convex due to the maximum work inequality but may not be strictly convex when a single point supports finite pressure. As shown in Fig. 3b, facets can occur since the object can rotate about one of the three support points whose velocity is zero with indeterminate underlying friction. Erdmann [11] proposed a configuration space embedding of friction. In his work, the third component of   F   is   F z   =   τ / ρ  and the third component of   V   is   V z   =   ωρ , where   ρ   is the radius of gyration. In doing so, all three components in   F  and   V   have the same unit. Observe that such normalized representation also obeys maximum work inequality with   ρ  being any characteristic length. In our experiments, we have found that the normalized representation yields better numer- ical condition and different values of   ρ   including radius of gyration, average edge length and minimum enclosing circle radius lead to similar performance. III. RELATED WORK Yoshikawa and Kurisu [12] solved an unconstrained least- squares problem to estimate the center of friction and the pressure distribution over discrete grids on the contact sur- face. With similar set up, Lynch [13] proposed a constrained linear programming procedure to avoid negative pressure assignment. However, methods based on discretization of the support surface introduce two sources of error in both localization of support points and pressure assignment among those points. We do not need to estimate the exact location of support points. Coarse discretization loses accuracy while fine discretization unnecessarily increases the dimensionality of estimation and model complexity. Howe and Cutkosky [14] presented an ellipsoid approximation of the limit sur- face assuming known pressure distribution. The ellipsoid was constructed by computing or measuring the major axis lengths (maximum force during pure translation and maxi- mum torque during pure rotation). Facets can be added by intersecting the ellipsoid with planes determined by each support point. The pressure distribution (except for 3 points support with known center of pressure), nevertheless, is non- trivial to measure or compute. We also show that the ellipsoid approximation, as a special case of our convex polynomial representation, is less accurate due to lack of expressiveness. Recent data-driven attempts [15], [16] collected visual data from random push trials and applied “off-the-shelf” machine learning algorithms to build motion models. We also embrace a data-driven strategy but bear in mind that physics princi- ples should guide the design of the learning algorithm (as constraints and/or priors), hence reducing sample complexity and increasing generalization performance. IV. REPRESENTATION AND IDENTIFICATION In this section, we propose the sublevel set representation of friction with desired properties and show that convex even- degree homogeneous polynomials are valid solutions. Then we formulate an efficient convex optimization procedure to identify such polynomials.  A. Polynomial sublevel set representation  Let   H ( F )   be a differentiable convex function that models the generalized friction load and velocity as follows:  •   The 1-sublevel set   L −  1   ( H ) =   { F   :   H ( F )   ≤   1 }   corresponds to the convex set   F   of all generalized friction loads.  •   The 1-level set   L 1 ( H ) =   { F   :   H ( F ) =   1 }   corresponds to generalized friction loads (during slip) on the boundary surface of   F   .  •   The surface normals given by gradients   { ∇ H ( F )   :   F   ∈  L 1 ( H ) }   represent instantaneous generalized velocity di- rections during slip, i.e.,   V   =   s ∇ H ( F )   where   s   >   0.  Theorem 1:   The set of friction loads represented by the 1-sublevel set of a differentiable convex function follows the maximum work inequality.  Proof:   When the object remains static,   F   belongs to the interior of   L −  1   ( H )   and   V   equals zero, the inequality holds as equality. When the object slips,   F   ∈   L 1 ( H )   and   V   is nonzero, we have for any other generalized friction load   F ′   ∈   L −  1   ( H ) :  V   ·   ( F ′   −   F ) =   s ( ∇ H ( F )   ·   ( F ′   −   F ))   ≤   s ( H ( F ′ )   −   H ( F ))   ≤   0 ,  where the first inequality is due to the convexity of   H ( F ) . In addition to enforcing convexity (discussed in IV-B), we choose   H ( F )   to obey the following properties: 1) Symmetry:   H ( F ) =   H ( − F )   and   ∇ H ( F ) =   − ∇ H ( − F ) . 2) Scale invariance:   ∇ H ( a F ) =   g ( a ) ∇ H ( F ) , where   g ( a )  is a positive scalar function. 3) Efficient   invertibility:   there   exists   efficient   numer- ical   procedure   to   find   a   F   ∈   L 1 ( H )   such   that  ∇ H ( F ) / ‖ ∇ H ( F ) ‖   =   V   for a given query unit velocity  V . We denote such operation as   F   =   H inv ( V ) . Symmetry is based on the assumption that negating the velocity direction would only result in a sign change in the friction load. Scale invariance is desired for two reasons: 1) scaling in mass and surface coefficient of friction could only result in a change of scale but not other geometrical properties of the level-set representation; and 2) predicting directions of generalized velocities (by computing gradients and   normalizing   to   a   unit   vector)   only   depends   on   the direction of generalized force. Such a property is useful in the context of pushing with robot fingers where applied loads are represented by friction cones. The inverse problem of finding the friction load for a given velocity naturally appears in seeking quasi-static balance for stable pushing or computing deceleration during free sliding, as shown in Section VI. In general, efficient numerical solution to the
inverse problem, which our representation enables, is key to planning and simulation. One solution family for   H ( F )   that obeys these properties is the set of   strongly convex even- degree homogeneous polynomials .  Theorem 2:   A strongly convex even degree- d   homoge- neous polynomial   H ( F ;   a ) =   ∑ m i = 1   a i F i 1  x   F i 2  y   F d − i 1 − i 2  z   with   m  monomial terms 2   parametrized by   a   satisfies the properties of symmetry, scale invariance, and efficient invertibility.  Proof:   Proving symmetry and scale invariance are trivial due to the homogeneous and even-degree form of  H ( F ) . Here, we sketch the proof that efficient invertibility can be achieved by first solving a simple non-linear least square problem followed by a rescaling. Construct an objective function   G ( F ) =   1 2   ‖ ∇ H ( F )   −   V ‖ 2  whose gradient   ∂   G  ∂   F   =   ∇ 2 H ( F )( ∇ H ( F )   −   V ) . Note that its stationary point   F ∗ , which iterative methods such as Gauss- Newton or trust-region algorithms will converge to, satisfies  ∇ H ( F ∗ )   −   V   =   0. Hence   F ∗   is globally optimal with value zero. Let   ∆ F t   =   ∇ 2 H ( F t   ) − 1 ( V t   −   V ) , then the update rule for Gauss-Newton algorithm is   F t + 1   =   F t   −   ∆ F t   . Although the final iteration point   F T   may not lie on the 1-level set of   H ( F ) , we can scale   F T   by ˆ F T   =   H ( F T   ) − 1 / d   F T   such that   H (   ˆ F T   ) =   1 and   ∇ H (   ˆ F t   ) / ‖ ∇ H (   ˆ F t   ) ‖   =   V   due to the homogeneous form of   H ( F ) . Therefore   H inv ( V ) =   ˆ F T   .  B. Sum-of-squares Convex Relaxation  Enforcing strong convexity for a degree-2 homogeneous polynomial   H ( F ;   A ) =   F T   A F   has a straightforward set up as solving a semi-definite programming problem with constraint of   A      ε I . Meanwhile, for a polynomial of degree greater than 2 whose hessian matrix   ∇ 2 H ( F ;   a )   is a function of both   F   and   a , certification of positive semi-definiteness is NP-hard. However, recent progress [17], [18] in sum-of- squares programming has given powerful semi-definite re- laxations of global positiveness certification of polynomials. Specifically, let   z   be an arbitrary non-zero vector in   R 3   and  y ( F ,   z ) = [ z 1 F x ,   z 1 F y ,   z 1 F z ,   z 2 F x ,   z 2 F y ,   z 2 F z ,   z 3 F x ,   z 3 F y ,   z 3 F z ] T   . If there exists a positive-definite matrix   Q   such that  z T   ∇ 2 H ( F ;   a ) z   =   y ( F ,   z ) T   Qy ( F ,   z )   >   0 ,   (4) then   ∇ 2 H ( F ;   a )   is positive definite for all non-zero   F   under parameter   a   and   H ( F ;   a )   is called as sos-convex. Further, equation (4) can be written as a set of   K   sparse linear constraints on   Q   and   a . Tr ( A k Q ) =   b T k   a ,   k   ∈ { 1   . . .   K }  Q      ε I ,   (5) where   A k   and   b k   are constant sparse element indicator matrix and vector that only depend on the polynomial degree   d . The number of constraints   K   equals 27 for   d   =   4.  C. Identification  This   section   sets   up   an   efficient   convex   optimization for identifying the coefficient   a   of the polynomial   H ( F ;   a )  given a set of measured noisy generalized force-motion  2 The number of different monomial terms   m   is bounded by   ( d + 2 2  ) .  { F i ∈{ 1 ... N } ,   V i ∈{ 1 ... N } }   pairs.   In   our   experiments,   we   use homogeneous 4th order polynomial. The optimization should find the coefficient   a   such that the measured   forces   F i  are close to the 1-level set surface and the corresponding gradients   are   aligned   well   (up   to   scale)   w.r.t   measured velocities   V i . Let   α i   =   || ∇ H ( F i ;   a )   −   ( ∇ H ( F i ;   a )   ·   V i ) V i || 2 2   be the L2-projection residual of   ∇ H ( F i ;   a )   onto the measured unit velocity vector   V i , and let   β i   = ( H ( F i ;   a )   −   1 ) 2   be a distance measurement of   F i   from the 1-level set of   H ( F i ;   a ) . We set up the optimization as follows: minimize  a , Q   ‖ a ‖ 2 2   +  N  ∑  i = 1  ( η 1 α i   +   η 2 β i )   (6) subject to   Tr ( A k Q ) =   b T k   a ,   k   =   1 , . . . ,   K ,   (7)  Q      ε I .   (8) The first term is for parameter regularization.   η 1   and   η 2   are trade-off parameters determined by cross-validation. Equa- tions (7) and (8) enforce convexity. Note that the objective is quadratic in   a   with sparse linear constraints and a semi- definite constraint on   Q . 3   We would like to point out that the formulation can be adapted online using projected gradient descent so that the importance of historical data is diminish- ing as the object moves, enabling the estimation to adapt to changing surface conditions. Evaluating such online version of the identification algorithm is deferred to future work. V. E XPERIMENTS  We conduct simulation and robotic experiments to demon- strate the accuracy and statistical-efficiency of our proposed representation. The model converges to a good solution with few available data which saves experimental time and design efforts.   We   compare   the   following   four   different   force- motion model representations   H   : 1) degree-4 convex homo- geneous polynomial (poly4-cvx); 2) degree-4 homogeneous polynomial (poly4) with convexity constraints 3) convex quadratic (quad) as degree-2 polynomial, i.e.,   H ( F ) =   F T   A F  with ellipsoid sublevel set; and 4) gaussian process (GP) with squared exponential kernel 4 . Denote by   V i   the ground truth instantaneous generalized velocity direction and   V p ( F i ;   H   )  as the predicted generalized velocity direction based on   H  for the input generalized load   F i , we use the average angle  δ   ( H   ) =   1  N   ∑ N i = 1   arccos ( V p ( F i ;   H   )   ·   V i )   between   V p ( F i ;   H   )  and   V i   as an evaluation criterion.  A. Simulation Study  Two kinds of pressure distribution are studied.  •   “Legged” support: Randomly sampled three support points on a unit circle with randomly assigned pressure.  •   “Uniform” support: Uniformly distributed 360 support points on a unit circle and 400 support points within a unit square. Each point has the same support pressure.  3 Code link:   https://github.com/robinzhoucmu/MLab_EXP/blob/ master/SlidingExpCode/LimitSurfaceFit/Fit4thOrderPolyCVX.m  4 The squared exponential kernel gives better performance over linear and polynomial. Normalizing the input load to a unit vector improves performance by requiring the GP to ignore scale. Every ( F ,   V ) input pair is augmented with ( − F ,   − V ) for training.
For each pressure configuration, we conduct 50 experimen- tal trials. To generate the simulated force-motion data, we assume a Coulomb friction model at each support point with a uniform coefficient of friction. Without loss of generality, sum of pressure over all contact points is normalized to one and the origin is set as the center of pressure. For each trial of “uniform” support, we sampled 150 instantaneous generalized velocities directions   V i   uniformly on the unit sphere and compute the corresponding generalized friction loads   F i . For each trial of “legged” support, 75   ( F i ,   V i )   pairs are uniformly sampled on the facets (same   V i   but different   F i  for each facet) and another 75 pairs are uniformly sampled in the same fashion as “uniform” support. In doing so, the dataset has a diverse coverage. Among the 150 pairs, 50% is used for hold-out testing, 20% is used for cross validation and four different amounts (7, 15, 22, 45) from the rest of 30% are used as training. In order to evaluate the algorithms’ robustness under noise, we additionally corrupt the   training   and   validation   set   using   Gaussian   noise   of standard deviation   σ   =   0 . 1 to each dimension of both   F i   and  V i   (renormalized to unit vector). From Fig. 2 we can reach the following conclusions. 1) Poly4-cvx has the smallest  δ   ( H   )   for different amounts of training data and pressure configurations. 2) Both poly4-cvx and convex quadratic show superior performance when data is scarce and noisy, demon- strating convexity is key to data-efficiency and robustness. Poly4-cvx model additionally shows larger improvement as more data is available due to stronger model expressiveness. 3) Poly4 (without convexity constraint) performs the worst when only few data is available, but gradually improves as more data is available for shaping the surface. 5   GP has similar performance trends as poly4 but worse on average. 4) Polynomial models enjoy significant performance advantages when limit surface is smoother as in uniform point support (approximation of uniform patch contact). Such advantage is smaller for three-points support whose limit surface has large flat facets.  B. Robotic Experiment  We mount three screws at four different sets of locations underneath an alluminium right-angle triangular work ob- ject 6 . Given known mass and COM projection, ideal ground truth pressure for each support point can be computed by solving three linear equations assuming each screw head approximates a point contact. Fig. 3a shows a flipped view of one arrangement whose ideal LS is illustrated in Fig. 3b, constructed by Minkowski addition of generalized friction at each single point support assuming Coulomb friction model with uniform coefficient of friction. Three pairs of symmetric  5 For noise-free experiments shown in Fig. 2b and 2d, when enough training data (more than 22) is presented, poly4 performs slightly better than poly4-convex. We conjecture such difference is due to the gap between sos-convex polynomials and convex polynomials.  6 The triangular object weighs 1.508kg with edge lengths of 150mm, 150mm and 212.1mm. The four different set of support point locations (in mm) with respect to the right angle corner vertex are: [(10,10), (10,130), (130,10)], [(30,30), (30,90), (90,30)], [(10,10), (10,130), (90,30)], [(30,30), (63.33,43.33), (43.33,63.33)].  facets 7   characterize indeterminate friction force when rotat- ing about one of the three support points. Comparison among identified fourth-order homogeneous polynomials with and without convexity constraint is shown in Fig. 3c and 3d. We can see that convex-shape constraint is essential to avoid poor generalization error when little data is available. Fig. 3e and 3f compare the level sets of a convex quadratic (ellip- soid) and a sos-convex degree-4 homogeneous polynomial, demonstrating that the higher degree polynomial captures the facets effect better than quadratic models. We conduct robotic poking (single point pushing) experi- ments on wood and paper board surfaces. In each experiment, we generate 50 pokes (30 for training set, 10 for validation set and 10 for test set) with randomly chosen contact points and pushing velocity directions. 8   Fig. 4 shows model accu- racy (averaged over four different pressure arrangments) with respect to increase in amount of training data for different methods evaluated on both the hold-out test sensor data and samples from ideal LS. We can see similar performance trends as in simulation experiments. Note that both evalu- ations only serve as certain reference criteria. Sensor data is noisy and all possible force measurements from a single point pusher only cover a limited space of the set of friction loads. We also do not intend to treat the idealized limit surface as absolute ground truth as there is no guarantee on uniform coefficient of friction between the support points and the underlying surface. Additionally, point contact and isotropic Coulomb friction model are only approximations of reality. Nevertheless, both evaluations demonstrate performance ad- vantage of our proposed poly4-cvx model. VI. A PPLICATIONS  A. Stable Push Action Generation  The resultant object velocity under a single point push action is hardly fully predictable. However a two-points push action against an edge of the object can be stable such that the object will remain attached to the pusher without slipping or breaking contact [2]. That is, the slider and pusher will move about the same center of rotation (COR) point   p c . Given the level set representation   H ( F ) , the condition of determining whether a two-points push with instantaneous generalized velocity   V p c   is stable or not is equivalent to check if the corresponding generalized friction force   F p c   =  H inv ( V p c   )   lies in the applied composite generalized friction cone   F c . To validate predictions based on the model, we sampled 60 random CORs and execute with the robot for three different pressure arrangements on a novel support surface material (hard poster paper). 9   15 out of the 60 CORs are labelled as stable. The training force-motion data are  7 The third one is in the back not visible from presented view.  8 During   each   pushing   action,   the   robot   moves   at   a   slow   speed   of 2.5mm/s with a total small push-in distance of 15mm. Each generalized velocity direction is approximated as the direction of pose displacement and generalized force is averaged over the action duration.  9 We use the same triangular block in Fig. 3a with two three-points contacts [(10,10), (10,130), (130,10)], [(30,30), (30,90), (90,30)] as well as full patch contact. The 60 CORs are tight rotation centers within a 400mm × 400mm square centered at the COM.
7   15   22   45  0  5  10  15  20  25  30  amount of training data velocity alignment error(degree)  Poly4cvx  Poly4  Quad  GP (a) Three support points with noisy   training   and   validation data. 7   15   22   45  0  5  10  15  20  25  amount of training data velocity alignment error(degree)  Poly4cvx  Poly4  Quad  GP  (b) Three support points with noise-free training and valida- tion data. 7   15   22   45  0  2  4  6  8  10  12  14  16  18  amount of training data velocity alignment error(degree)  Poly4cvx  Poly4  Quad  GP  (c)   Uniform   circular   support points with noisy training and validation data. 7   15   22   45  0  1  2  3  4  5  6  7  8  9  10  11  amount of training data velocity alignment error(degree)  Poly4cvx  Poly4  Quad  GP  (d)   Uniform   circular   support points with noise-free training and validation data.  Fig. 2:   Test error comparison for simulation experiments with 95% confidence bar (50 random evaluations) among different methods as amount of training data increases for three random support points and 360 support points on a ring respectively. Results for uniform pressure distribution within a square are similar to uniform circular support and omitted for space.  (a)   Triangular   block   with three support screws. −4   −2   0   2   4   −4   −2   0   2   4  −3  −2  −1  0  1  2  3  F y F x  F z  (b) Ideal limit surface with facets.  (c) Poly4 fit with 5 training and 5 validation data.  (d) Poly4-cvx fit with 5 training and 5 validation data.  (e) Convex quadratic fit with 10 training and 10 validation data.  (f) Poly4-cvx fit with 10 train- ing and 10 validation data.  Fig. 3:   Level set friction representations for the pressure arrange- ment in Fig. 3a. Red dots and blue arrows are collected generalized forces and velocities from force-torque and motion capture sensor respectively. Fig. 3c and 3d, Fig. 3e and 3f share the same data.  collected from pushing the object on a wood surface. Table I and II summarize the classification accuracy and positive (stable) class recall measurements of three invertible methods with respect to increase in amount of training data. Fig. 5 shows an example (full patch contact) that the stable regions generated from the identified poly4-cvx model is much larger than the conservative analysis as in [2] which misses the tight/closer rotation centers. TABLE I:   Comparison of average accuracy with 95% confidence interval as amount of training data increases.  10   20   30 poly4-cvx   88.13 ± 1.80   91.33 ± 1.61   93.07 ± 1.45 poly4   85.27 ± 2.12   89.40 ± 1.98   93.00 ± 1.62 quadratic   87.93 ± 1.72   87.20 ± 1.65   88.00 ± 1.39  TABLE   II:   Comparison   of   average   positive   recall   with   95% confidence interval as amount of training data increases.  10   20   30 poly4-cvx   90.13 ± 3.54   96.69 ± 1.93   98.18 ± 1.32 poly4   79.96 ± 5.25   92.76 ± 2.90   97.18 ± 1.84 quadratic   73.18 ± 4.61   73.38 ± 4.69   73.87 ± 4.63  B. Free Sliding Dynamics Simulation  Given   H ( F ) , the equation of motion (with respect to the object local coordinate frame) during free sliding assuming a uniform surface can be written as   I   dV dt   =   − H inv ( V ) , where  I   is the moment of inertia. We use the Runge-Kutta method provided in MATLAB ODE45 and demonstrate several ex- ample sliding trajectories in Fig. 6. As studied in [19], given an ideal limit surface, a free sliding object comes at rest with one of several definite generalized velocity directions (in local body frame), termed as eigen-directions. We have empirically found a similar trend that there exists multiple converging sets of initial generalized velocities. An example behavior is shown in Fig. 6d where the final velocity directions (instantaneous rotation centers) remain in the same small region regardless of different initial velocities. VII. CONCLUSION AND FUTURE WORK In this paper, we propose to use the sub-level sets and gradients of a function to represent rigid body planar fric- tion loads and velocities, respectively. The maximum work inequality implies that such a function needs to be convex. We additionally require the properties of symmetry, scale invariance, and efficient invertibility which lead us to choose a convex even-degree homogeneous polynomial representa- tion. We apply the representation to applications including stable pushing and dynamic simulation. For future work, we plan to evaluate the model on a larger dataset with varying object and surface material physical properties. We will also explore methods for online model identification.
5   10   20   30  0  5  10  15  20  25  30  35  40  45  amount of training data velocity alignment error(degree)  Poly4cvx  Poly4  Quad  GP (a) Test on sensor data (wood surface). 5   10   20   30  0  5  10  15  20  25  30  35  40  45  amount of training data velocity alignment error(degree)  Poly4cvx  Poly4  Quad  GP  (b) Test on data sampled from ideal LS (wood surface). 5   10   20   30  0  5  10  15  20  25  30  35  amount of training data velocity alignment error(degree)  Poly4cvx  Poly4  Quad  GP  (c) Test on sensor data (paper board surface). 5   10   20   30  0  5  10  15  20  25  30  35  40  amount of training data velocity alignment error(degree)  Poly4cvx  Poly4  Quad  GP  (d) Test on data sampled from ideal LS (paper board surface).  Fig. 4:   Test error comparison for robotic experiments with 95% confidence bar (50 random evaluations) among different methods as amount of training data increases for three support points on wood and hard paper board surfaces. −500   0   500  −400  −300  −200  −100  0  100  200  300  400  X/mm Y/mm  Fig. 5:   Hatched areas correspond to stable CORs region based on the conservative analysis [2]. Red triangles are stable CORs and gray stars are non-stable CORs based on the poly4-cvx model. The two push points are 50mm in width. The pusher and objects are covered with electrical tape and guffer tape respectively with measured coefficient of friction equals one.  ACKNOWLEDGMENT This work was conducted in part through collaborative par- ticipation in the Robotics Consortium sponsored by the U.S Army Research Laboratory under the Collaborative Technol- ogy Alliance Program, Cooperative Agreement W911NF-10- 2-0016 and National Science Foundation IIS-1409003. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory of the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for   Government   purposes   notwithstanding   any   copyright notation herein. R EFERENCES [1]   M.   T.   Mason,   “Mechanics   and   planning   of   manipulator   pushing operations,”   IJRR , vol. 5, pp. 53–71, Fall 1986. [2]   K. M. Lynch and M. T. Mason, “Stable pushing: Mechanics, control- lability, and planning,”   IJRR , vol. 15, pp. 533–556, Dec. 1996. [3]   S. Akella and M. T. Mason, “Posing polygonal objects in the plane by pushing,”   IJRR , vol. 17, pp. 70–88, January 1998. [4]   M. Dogar and S. Srinivasa, “Push-grasping with dexterous hands: Mechanics and a method,” in   IROS , 2010. [5]   A. Cole, P. Hsu, and S. Sastry, “Dynamic control of sliding by robot hands for regrasping,” in   IEEE Transactions on robotics and automation , vol. 8, 1992. [6]   D. E. Whitney, “Quasi-static assembly of compliantly supported rigid parts,”   ASME Journal of Dynamic Systems, Measurement, and Control , vol. 104, pp. 65–77, Mar. 1983. [7]   M.   T.   Mason,   “On   the   scope   of   quasi-static   pushing,”   in   ISRR , Cambridge, Mass: MIT Press, 1986. −100   −50   0   50   100  −80  −60  −40  −20  0  20  40  60  80  100  X/mm Y/mm  (a) Trajectory with   V ( 0 )   equals (150mm/s, -150mm/s, 2 π rad/s). −100   −50   0   50   100  −80  −60  −40  −20  0  20  40  60  80  100  X/mm Y/mm  (b) Trajectory with   V ( 0 )   equals (250mm/s, -100mm/s, 3 π rad/s). −100   −50   0   50   100   150   200  −50  0  50  100  150  200  X/mm Y/mm  (c) Trajectory with   V ( 0 )   equals (500mm/s, 100mm/s,   π rad/s).  (d) Trajectories of generalized friction loads.  Fig. 6:   Dynamic simulations based on a polynomial model using 10 training and 10 validation data for a triangular object with three support points as in Fig. 3a. Fig. 6a, 6b, and 6c illustrate the sliding trajectories for different initial velocities. The rgb and black triangles are the initial and final poses respectively, with dotted triangles as intermediate poses. Black dots traces correspond to the trajectories of center of mass. Magenta circles are the instantaneous rotation centers at the final time steps. Trajectories of the generalized friction loads are shown in Fig. 6d, where red, black, and blue curves correspond to Fig. 6a, 6b, and 6c respectively.  [8]   J. J. Moreau, “Unilateral contact and dry friction in finite freedom dynamics,” in   Nonsmooth Mechanics and Applications , pp. 1–82, Springer, 1988. [9]   S. Boyd and L. Vandenberghe,   Convex Optimization .   Cambridge university press, 2004. [10]   S. Goyal, A. Ruina, and J. Papadopoulos, “Planar sliding with dry friction. Part 1. Limit surface and moment function,”   Wear , vol. 143, pp. 307–330, 1991. [11]   M. Erdmann, “On a representation of friction in configuration space,”  IJRR , vol. 13, no. 3, pp. 240–271, 1994. [12]   T. Yoshikawa and M. Kurisu, “Identification of the center of friction from pushing an object by a mobile robot,” in   IROS , 1991. [13]   K. M. Lynch, “Estimating the friction parameters of pushed objects,” in   IROS , 1993. [14]   R. D. Howe and M. R. Cutkosky, “Practical force-motion models for sliding manipulation,”   IJRR , vol. 15, no. 6, pp. 557–572, 1996. [15]   M. Kopicki, S. Zurek, R. Stolkin, T. Morwald, and J. Wyatt, “Learning to predict how rigid objects behave under simple manipulation,” in  ICRA , 2011.
[16]   D. Omrcen, C. Boge, T. Asfour, A. Ude, and R. Dillmann, “Au- tonomous acquisition of pushing actions to support object grasping with a humanoid robot,” in   Humanoids 2009 , 2009. [17]   P. A. Parrilo,   Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization .   PhD thesis, California Institute of Technology, 2000. [18]   A. Magnani, S. Lall, and S. Boyd, “Tractable fitting with convex polynomials via sum-of-squares,” in   CDC-ECC , 2005. [19]   S. Goyal, A. Ruina, and J. Papadopoulos, “Planar sliding with dry friction. Part 2. Dynamics of motion,”   Wear , vol. 143, pp. 331–352, 1991.