1 Introduction
As frictional contact is the fundamental driving process by which many robots are able to interact with their surroundings, it is unsurprising that its behavior is central to a large body of robotic locomotion and manipulation research (e.g. [7, 11, 19, 20, 23, 27]). However, dynamical models of these systems are inherently complex and challenging to simulate and analyze. Impacts between rigid bodies induce instantaneous jumps in velocity states and a combinatorial explosion of hybrid modes that in conjunction render application of common tools from control theory and trajectory optimization difficult. While there has been notable progress in planning through unknown contact sequences with full dynamics [19, 15, 13], model complexity has thus far still inhibited realtime usage.
Many applications in robotics involving frictional contact exhibit structure that permits partial or full circumvention of these difficulties. Particularly, we examine planar tabletop manipulation, where a manipulator effects motion of an object that rests upon a flat, frictional surface. Several results in simulation, control, and planning for such systems have been enabled by quasistatic assumptions—that if manipulator accelerations and velocities are low, a forcebalance equation can approximate Netwon’s second law. This assumption enables reducedorder modeling of the object’s movement in response to contact with a manipulator driven at a particular velocity; such models also often circumvent the complexity associated with state discontinuities in dynamical approaches. Furthermore, quasistatic models often eliminate numerical sensitivity induced by the stiff dynamics of manipulators with high feedback gain controllers.
The tractability of these models has enabled impressive results in formal control analysis, task planning, and learning (e.g. [5, 10, 12]). However, the range of motion that these methods are currently able to model is limited. They are often restricted to pushing and nonprehensile motions and are completely unable to usefully express grasping or jamming; in these cases, their associated mathematical programs often yield no solutions or ambiguous behavior. Grasping and jamming objects is crucial for a wide range of robot tasks, and much work has been devoted to planning and controlling action before and after grasping events (e.g.[22, 21, 9, 30, 18]); However, much of this work can only describe static grasp configurations, and is unable to depict grasplike behavior with sliding contacts. The process of acquiring such a grasp itself often involves jamming (e.g. when reorienting an object by pushing it up against a wall), which neither the prehensile nor static grasping models can capture. We therefore find great value in the formulation of a unified quasistatic model that can smoothly capture a complete task involving the acquisition and use of a grasp. The ambiguity in traditional approaches arises from the inconsistent assumptions of rigid bodies and perfect control of manipulator velocity. Our key insight is that by appropriately representing the manipulator’s internal controller, physicallygrounded motion of the objectmanipulator system is guaranteed to exist. We contribute both instantaneous velocity and timestepping position models that formulate this behavior as linear complementarity problems, and prove existence of solutions for each.
2 Related Work
There is a significant body of research that examines manipulation from a quasistatic perspective [14, 27, 11, 23, 17]. For systems in which the object experiences frictional support, relevant research typically examines the pressure distribution supporting the object. Some earlier works provide guaranteed properties of the object’s motion without full knowledge of this distribution. Mason [14] derived the voting theorem to construct a mapping from the center of pressure to the direction of the object’s angular velocity. Lynch and Mason [12] later performed stability and controllability analysis for a manipulator pushing an object with multiple fingers. Other works alternatively contribute models that directly map manipulator joint velocities to object motion. Trinkle [27] characterized vertical planar manipulation using a nonlinear mathematical program that explicitly solved for the contact forces between the object and the manipulator. A similar, more general model for arbitrary 3D rigid multibody systems was proposed in Trinkle et al. [28]. Neither formulation can model detailed pressure distributions between surfaces, as they model friction as acting at a finite set of points.
Efficient and expressive modeling of the complex behavior of these pressure distributions was enabled by Goyal et al. [8], who define the limit surface, the bounded convex set of friction loads that the surfaces in contact might exert on each other. Zhou et al. [31] approximate this set in high detail as a semialgebraic set fitted to experimental data, and from it derive a model for the motion of the manipulatorobject system. This method, however, assumes that the manipulator follows a commanded velocity exactly, and therefore does not return a solution for infeasible commands. As such commands often result in grasping in a full dynamical model, valuable behaviors are not captured. Some planning methods (e.g [6]) introduce manipulator compliance through the generalized dampers of Whitney [29], though they have not been incorporated into models capable of initiating and releasing multiple contacts.
Pang and Tedrake [17] also devise a resolution for nonexistence in velocitycontrolled 3D rigid quasistatic systems. They model deviations from desired velocity as a result of local elastic deformation at point contacts, and preserve realism by minimizing them in a mixedinteger quadratic program (MIQP) formulation. While their work applies to a more general class of systems, our method has three key advantages for planar manipulation: we draw model behavior from a problem class that is far more tractable than MIQPs; our proof of existence makes possible formal guarantees for controller performance as well as simulation reliability; and our inclusion of a limitsurface model allows for more realistic modeling without introducing significant complexity.
3 Background
We now describe the behavior of a frictionallysupported object with center of mass in contact with a manipulator, depicted in Figure 1. For a more detailed treatment of rigidbody contact dynamics, we refer the reader to [1]. Coordinates of the object and manipulator configurations are denoted and , respectively. The manipulator is assumed to be controlled to track commanded velocity . Frictional contact with the manipulator causes the object to move with bodyaxis velocity , which can be converted to the world frame velocity with the fullrank transformation . Contact is modeled as normal and frictional forces and at pairs of points on the boundaries of the object and manipulator. The coefficients of friction at these points are given by the diagonal matrix . Each contact’s tangential force is split into two nonnegative opposing forces and (as in [1]), such that the net tangential force is equal to
. The distances between each pair of points are represented by the vector
. When the points are in contact, their separating velocity can be calculated from the generalized velocities as using Jacobian matrices and . Similarly, there exist Jacobian matrices and such that the velocity at which the bodies slide against each other is given by . For notational convenience, we will often group the normal and tangential terms as , , and . Using the relationship between the contact velocities and the generalized velocities, the contact forces can be used to determine the forces acting on the object and manipulator as and , respectively.We also make use of the constants and the blockdiagonal matrix
(1) 
3.1 Linear Complementarity Problems
Throughout this work, we will make regular use of linear complementarity problems (LCPs). An LCP is a particular type of mathematical program for which solutions can be efficiently computed. LCPs have been widely used by the dynamics and robotics communities for describing the effects of contact (e.g. [25, 1, 31]). Here, we briefly introduce the problem formulation and some useful properties, and we refer the reader to [3] for a more complete description.
The linear complementarity problem for the matrix and vector describes the mathematical program
find  (2)  
subject to  (3)  
(4)  
(5) 
for which the set of solutions is denoted . Constraints (3)–(5), called complementarity constraints, are often abbreviated as .
Note that vector inequalities in the above definition, as well as elsewhere in this work, are taken elementwise. We will find that for LCPs related to frictional behavior, the matrix parameter is often copositive (i.e. for all ). This property is often theoretically useful, as Corollary 4.4.12 of [3], reproduced below, gives a sufficient condition for copositive LCP feasibility. Let , and let be copositive. Suppose that for every , we have . It follows that , and an element of can be discovered by Lemke’s Algorithm in finite time.
3.2 Friction at Point Contacts Between Object and Manipulator
Common to many models of friction is the maximum dissipation principle, which states that if
is the set of all possible forces acting at a contact, then the force at any given moment is one such that the power dissipated at the contact is maximized. For contact at a single point between two rigid bodies with relative velocity
, this condition is realized as(6) 
For point contacts, is often modeled as a cone by Coulomb friction, which enforces the following:

For each contact, either the normal velocity is zero and the normal force is nonnegative, or vice versa:
(7) 
The magnitude of the frictional force at the th contact is bounded above by . For sliding contacts, the frictional force has magnitude and is antiparallel to the sliding velocity. This behavior can be captured as complementarity constraints with the addition of a slack variable :
(8) (9)
3.3 Friction at Contact Between Object and Surface
Coulomb friction behavior cannot readily be applied to contacts where the normal force is not concentrated at a finite set of points. Zhou et al. [31] therefore devise and experimentally validate a model that directly approximates the limit surface assuming a constant pressure distribution. They parameterize this behavior with a symmetric, scaleinvariant, and strictly convex function , defined over the space of bodyaxis friction wrenches . The physical meaning of this function is as follows:

The set of possible static friction wrenches is .

If the contact is sliding, then the maximum dissipation principle requires that be on the boundary of , . Furthermore, must lie in the normal cone of at .
As is strictly convex, the latter condition is exactly satisfied by
(10) 
We note that (10) will also hold in the case of static friction with .
3.4 Friction Behavior as an LCP
We examine the quasistatic model of Zhou et al. [31], which composes both point and surface contacts. From (10), if has the ellipsoid form ,
(11) 
(12) 
where . Assuming perfect velocity control (i.e. ), (7)–(9) reduce to , a generalization of (27) in [31], where
(13)  
(14) 
Each is a choice for that complies with both the point and surface contact models. While solutions for and are not computed separately, their product is sufficient to calculate object velocity. For simplicity, we will denote as for the rest of this paper. We also define , the set of feasible point contact forces for command velocity . For , is the set of admissible internal forces [16], i.e. forces that map to zero net force and torque on the object ().
For nonquadratic descriptions of , low accuracy solutions may be computed quickly by approximating as an ellipsoid. If higher accuracy solutions are required, one may solve a sequence of programs such that in the th program is equal to the Hessian of evaluated at a solution of th program.
4 Finite Velocity Feedback Quasistatics
While the above formulation has been successful at simulating pushes [31] and planning grasps under stochasticity [30], the range of applications of this method is significantly limited due to undefined and ambiguous behaviors as displayed in Figure 2. Grasping and jamming commands may result in manipulator velocities that cannot be realized without penetrating (see Lemma 1, [31]). Additionally, when the manipulator is commanded to graze the object (that is, the commanded velocities of all the contact points are parallel to the object boundary), there can be an infinite set of possible solutions which result in wildly different motions.
The source of this ambiguity is a critical assumption of the perfect velocity control model—that there exists a feedback controller internal to the manipulator that has high enough (essentially, infinite) gain to overcome any external disturbance. From this perspective, grasping and jamming are undefined as they prescribe two unstoppable objects to oppose each other (””like behavior). Additionally, while exact execution of a grazing maneuver may produce zero external disturbance for the manipulator controller to balance, a small perturbation of the commanded velocities towards the inward normal of the object surface could transform the command to an infeasible grasp or jam.
In order to resolve these issues, we explicitly interpret the quasistatics of perfect velocity control as the limiting case of the quasistatics of finite, stable, linear feedback. That is, we assume some relative feedback gain matrix and scaling factor such that the generalized force due to contact exerted on the manipulator is balanced by a feedback torque
(15) 
where we note that the gain is inversely proportional to and the scaling of . While (13) and (14) assumed that the manipulator velocity directly tracked the desired velocity, , we instead solve (15) for and construct a new program that accounts for the modified velocity:
(16)  
(17)  
(18) 
While we have eliminated the perfect velocity control assumption, our formulation introduces no additional computational complexity, as our LCP is of equal dimension. We now prove that our model is wellbehaved in the sense that solutions are guaranteed to exist (Theorem 4); the manipulator velocity remains bounded as (Theorem 4); and solutions typically converge the perfect velocity control case when it is feasible (Theorem 4).
For all , if and , . Let . As all the indices of are nonnegative, . We also have that and as and . Therefore,
(19) 
and is copositive. Now, suppose . We must have
(20) 
which implies . Therefore, . The result follows from Proposition 3.1.
If the desired velocity lies within a neighborhood of infeasible velocities (symptomatic of grasping or jamming behavior), then as we take , the force balance in (15) implies that , and therefore some of the individual finger forces, will grow unboundedly. However, the net manipulator velocity error remains bounded due to the boundedness of : Let . There exists such that for all , for all , . Let . From Definition 3.1, we have that
(21)  
(22) 
Let
be the minimum eigenvalue of
. Noting , we have that(23)  
(24)  
(25)  
(26) 
Intuitively, if individual finger forces grow without bound, yet the net force remains bounded, then there must be a “canceling out” effect. More precisely, the portion of that experiences this growth must be an internal force:
For all sequences and such that and , Let . by constructing from . It is therefore sufficient to show .
We first note that , and follow directly from . Multiplying (22) by , we have
(27) 
In the limit, the second and fourth terms in this equation vanish due to the boundedness of , rendering
(28) 
As , , and therefore .
can therefore be considered a basis from which one can generate the errors in the manipulator velocity; for sufficiently close to , there exists a , such that . When there are no nonzero internal forces, the manipulator displacement approaches zero, and the solution of the finite linear feedback model approaches the behavior for perfect velocity control:
Suppose that is chosen such that . For all sequences , such that and , then . Assume the contrary, so that there exists and bounded away from such that . Letting , we have that . , given by Theorem 4, implies . For all , we observe the complementarity condition
(29) 
which implies that is bounded for bounded . As , for all , we have . Therefore by Proposition 3.1, and are nonempty.
If there were a subsequence such that were bounded, then from (29) and the nonemptiness of , would have a limit point in , violating our assumptions. Therefore, we must have . Similar to Theorem 4, by dividing (29) by ,
(30) 
and thus . But , a contradiction.
We note that implies the current configuration is not a forceclosure [16]. In practice, we expect the set of nonforceclosure configurations that have to be small. We also expect perfect velocity control models to behave poorly during forceclosure, as they exhibit grasping behavior. Therefore, this model will behave more realistically in a variety of scenarios, with minimal accuracy loss for prehensile commands that perfect velocity control models handle well already.
5 TimeStepping Scheme
Despite the guarantee (Theorem 4) that (17) and (18) provide feasible instantaneous velocity solutions, embedding the LCP into common ODE schemes is an incomplete approach, as the resulting velocities will be discontinuous whenever contact is initiated between the manipulator and objects. Anitescu and Potra [1] resolved a similar issue in their formulation for 3D multibody simulation by rootfinding the first subtimestep impact. While a similar modification could be applied to our LCP, the ability to resolve subtimestep impacts in a single LCP would be beneficial. To that end, we take inspiration from Stewart and Trinkle [25], and instead formulate an alternative LCP that explicitly models the positions of the manipulator and object at the end of the timestep.
To construct this new program, instead of solving for the force at a particular time , we solve for the the net normal and tangential impulse between times and at each contact, and . Using the superscripts and to denote values calculated at the beginning and end of this interval, we linearize as
(31) 
and we make a firstorder approximation of and :
(32)  
(33) 
and are subject to the complementarity constraints
(34)  
(35)  
(36) 
Arranging into the standard format, we arrive at , where
(37)  
(38) 
With the exception of the added in , this program is identical to an instantiation of the velocity formulation defined in (17) and (18). Noting that for any feasible initial condition, a trivial extension of Theorem 4 guarantees existence of solutions for all feasible initial conditions. This is a significant improvement over existing timestepping schemes for the full dynamic behavior, as it circumvents both the expensive rootfinding subroutine of [1] and the nonexistance issues in [25]. However, it is important to note that if is nonlinear, the linearization in (31) does not guarantee that the true value of will be feasible (positive), even if is feasible. In these cases, similar to [25]
, one may rectify this issue by solving a sequence of problems in a fixedpoint iteration scheme, linearizing the problem about the best current estimate for
. This iteration can be conducted at the same time as the iteration for a nonellipsoidal . We additionally note that setting gives a timestepping reformulation of the program in [31].6 Examples
We provide a few examples to illustrate the capabilities and accuracy of our approach. Three examples from our open source MATLAB library^{1}^{1}1https://github.com/mshalm/quasistaticGrasping are provided in conjunction with a video depiction^{2}^{2}2https://www.youtube.com/watch?v=1wAH5o3OLck. Additionally, we compare the output of our model to a fully dynamic, compliant simulation using Drake [26]. All LCPs associated with our model are solved with PATH [4].
6.1 Pushing with Two Fingers
We first consider a flat disk of radius pushed by two fullyactuated point fingers. The coefficient of friction between the fingers and the disk is set at , and the forcemotion map is set as such that . We consider the configuration which directly represents the coordinates of each finger in a fixed frame. We assume each coordinate is controlled independently with equal gain, such that .
We now empirically evaluate the validity of Theorem 4. Analytically, we expect that as we take , finite feedback simulation should converge to perfect velocity command tracking. We simulate three motions over at for various feedback scaling terms , and plot the corresponding object trajectories in Figure 3. We compare the final configuration of each finite feedback trajectory with , the trajectory resultant from executing the same manipulator command with a perfect velocity control. For sufficiently high gains, we see a linear relationship between the log of the feedback scaling term and that of , the error in the final pose.
However, the lowgain performance, particularly on the semicircular command, showcases an important nuance in the convergent behavior described in Theorem 4. While the method may converge over a single timestep, differences in individual timesteps may accumulate such that there is a change in the contact mode during the motion. In these lowgain cases, the fingers tend to slide off of the sides of the object, leaving the object behind its intended goal.
We now show our formulation’s ability to simulate through jamming motions that perfect velocity control models cannot capture. We attempt to push the object into and roll along a wall. We use identical simulation parameters, and plot the results in Figure 4. We can see that the perfect velocity control simulation terminates at the collision with the wall, while the finite gain formulation not only captures the collision, but also still exhibits convergent behavior thereafter.
6.2 Polygonal Peginhole
We apply our method to a peginhole problem, a more complex task requiring initiation and release of several contacts and gripping motions. We simulate a new system consisting of a thin rectangular peg, two triangular manipulator fingers, and a slot into which the peg is fit. The slot has a tolerance on the width of the peg, the relative feedback gains are set to (in units for linear terms and for rotational terms), and the feedback scaling is set to . In a handdesigned trajectory, the peg catches the corner of the slot in the initial insertion attempt, after which the manipulator reorients and successfully inserts it. For each time step, only bodies close enough to make contact are considered, allowing the LCPs to be kept small. We set our time step to , for which PATH is able to compute solutions in on average. The trajectory is displayed in Figure 4.
6.3 Comparison to Full Dynamics
We now evaluate the similarity of our model to Drake [26], a framework for simulation of rigidbody dynamics with contact. We simulate a square object of mass and side length in contact with four manipulators, each consisting of two long links with a circular contact at the end of the second link. The manipulators pinch the object, which is particularly numerically challenging to simulate. While the velocity commands are symmetric, each manipulator is driven with different feedback gains, causing the object to move in the direction and spin in response to the pinch. The dynamic and quasistatic simulations exhibit qualitatively similar behavior, shown in Figure 5.
7 Discussion and Future Work
We have presented a method for generating motion using a quasistatic model for planar manipulation. By explicitly modeling manipulator feedback behavior, our method is able to synthesize physicallygrounded motion in the face of infeasible velocity commands. Despite this added capability, our method is as computationally efficient as previous methods. We have validated our model theoretically and empirically by proving convergence to an established model.
An interesting property of many quasistatic manipulation models is that they allow decoupling of manipulator dynamics from object motion. However, we assert that Theorem 4 implies that controller choice inherently effects grasping equilibrium and is therefore an unavoidable component of a physicallygrounded grasping model. We also view the decoupling of the relative feedback gains from the overall feedback intensity as essential for the accuracy of our model. As the force scaling term of (10) is not explicitly determined, even given true manipulator gains, one cannot synthesize with complete accuracy. However, in light of Theorems 4 and 4, if the controller has high enough gain, a small enough choice for will produce accurate results. In these cases, is a numerical term rather than something physically meaningful; one should choose as small as possible without degrading numerical precision in the construction of .
One might suspect that embedding a highgain controller into the model induces stiff behaviors, as the corresponding full dynamics are stiff. However, our quasistatic assumptions happen to eliminate this behavior. As Theorem 4 proves convergence to the perfect velocity control model in most cases and Theorem 4 proves that velocities are bounded by otherwise, our feedback terms does not add significant stiffness. Some numerical precision may however be lost for small due to roundoff error or poor conditioning of the associated LCPs. In future work, we will conduct a quantitative analysis of this behavior.
We do not expect solutions to our LCPs to be unique, as nonuniqueness is pervasive in complementaritybased contact models [2, 24]. While our model does construct a unique mapping between and , there are unmet assertions required for the mapping between and to be unique. However, it does disambiguate some grazing cases (such as in Figure 2).
Future extension of this result to 3D motion poses significant challenges. In the 2D case, the contact between the object and the surface below generates a unique map from contact forces to object motion. In 3D motion, the manipulator must instead counteract gravity. Furthermore, quasistatic modeling cannot realistically capture certain actions; for instance, dropping the object may either result in a lack of solution, or in a large enough to make linearization of inaccurate. Possible applications of this model include controller synthesis through sumsofsquares based Lyapunov analysis and model predictive control, as well as planning via trajectory optimization methods.
8 Acknowledgements
This material is based upon work supported by the National Science Foundation under Grant No. CMMI1830218. The authors thank Joury Abdeljaber for her work on numerical experiments.
References
 [1] Anitescu, M., Potra, F.A.: A timestepping method for stiff multibody dynamics with contact and friction. International Journal for Numerical Methods in Engineering 55(7), 753–784
 [2] Brogliato, B.: Nonsmooth Mechanics. Springer London (1999)
 [3] Cottle, R.W., Pang, J.S., Stone, R.E.: The Linear Complementarity Problem. Society for Industrial and Applied Mathematics (Jan 2009)
 [4] Dirkse, S.P., Ferris, M.C.: The path solver: a nommonotone stabilization scheme for mixed complementarity problems. Optimization Methods and Software 5(2), 123–156 (1995)
 [5] Dogar, M.R., Srinivasa, S.S.: A planning framework for nonprehensile manipulation under clutter and uncertainty. Auton. Robots 33, 217–236 (2012)
 [6] Erdmann, M.: Using backprojections for fine motion planning with uncertainty. The International Journal of Robotics Research 5(1), 19–45 (1986)
 [7] F. Hogan, E.R., Rodriguez, A.: Reactive planar manipulation with convex hybrid mpc. In: Proceedings 2018 IEEE International Conference on Robotics and Automation (2018)
 [8] Goyal, S., Ruina, A., Papadopoulos, J.: Planar sliding with dry friction part 1. limit surface and moment function. Wear 143(2), 307 – 330 (1991)
 [9] HaasHeger, M., Iyengar, G., Ciocarlie, M.: Passive reaction analysis for grasp stability. IEEE Transactions on Automation Science and Engineering 15(3), 955–966 (July 2018)
 [10] Kloss, A., Schaal, S., Bohg, J.: Combining learned and analytical models for predicting action effects. CoRR abs/1710.04102 (2017), http://arxiv.org/abs/1710.04102
 [11] Lynch, K.M.: The mechanics of fine manipulation by pushing. In: Proceedings, 1992 IEEE International Conference on Robotics and Automation. pp. 2269–2276 vol.3 (May 1992)
 [12] Lynch, K.M., Mason, M.T.: Stable pushing: Mechanics, controllability, and planning. The International Journal of Robotics Research 15(6), 533–556 (1996)
 [13] Manchester, Z., Kuindersma, S.: Variational contactimplicit trajectory optimization. In: International Symposium on Robotics Research (ISRR). Puerto Varas, Chile (2017)
 [14] Mason, M.T.: Mechanics and planning of manipulator pushing operations. The International Journal of Robotics Research 5(3), 53–71 (1986)
 [15] Mordatch, I., Lowrey, K., Todorov, E.: Ensemblecio: Fullbody dynamic motion planning that transfers to physical humanoids. In: 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). pp. 5307–5314 (Sept 2015)
 [16] Murray, R.M., Sastry, S.S., Zexiang, L.: A Mathematical Introduction to Robotic Manipulation. CRC Press, Inc., Boca Raton, FL, USA, 1st edn. (1994)
 [17] Pang, T., Tedrake, R.: A robust timestepping scheme for quasistatic rigid multibody systems (2018)
 [18] Paolini, R., Rodriguez, A., Srinivasa, S., Mason, M.T.: A datadriven statistical framework for postgrasp manipulation. International Journal of Robotics Research (IJRR) 33(4), 600–615 (April 2014)
 [19] Posa, M., Cantu, C., Tedrake, R.: A direct method for trajectory optimization of rigid bodies through contact. The International Journal of Robotics Research 33(1), 69–81 (2014)
 [20] Posa, M., Tobenkin, M., Tedrake, R.: Stability analysis and control of rigidbody systems with impacts and friction. IEEE Transactions on Automatic Control (TAC) 61(6), 1423–1437 (2016)
 [21] Rimon, E., Burdick, J.: On force and form closure for multiple finger grasps. In: Proceedings of IEEE International Conference on Robotics and Automation. vol. 2, pp. 1795–1800 vol.2 (Apr 1996)
 [22] Rodriguez, A., Mason, M.T., Ferry, S.: From caging to grasping. Int. J. Rob. Res. 31(7), 886–900 (Jun 2012)
 [23] Song, P., Kumar, V., Trinkle, J.C., Pang, J.S.: A family of models for manipulation planning. In: (ISATP 2005). The 6th IEEE International Symposium on Assembly and Task Planning: From Nano to Macro Assembly and Manufacturing, 2005. pp. 236–241 (July 2005)
 [24] Stewart, D.E.: Rigidbody dynamics with friction and impact. SIAM Rev. 42(1), 3–39 (Mar 2000)
 [25] Stewart, D., Trinkle, J.: An implicit timestepping scheme for rigid body dynamics with inelastic collisions and coulomb friction. International Journal of Numerical Methods in Engineering 39, 2673–2691 (1996)
 [26] Tedrake, R., the Drake Development Team: Drake: A planning, control, and analysis toolbox for nonlinear dynamical systems (2016), https://drake.mit.edu
 [27] Trinkle, J.C.: A quasistatic analysis of dextrous manipulation with sliding and rolling contacts. In: Proceedings, 1989 International Conference on Robotics and Automation. pp. 788–793 vol.2 (1989)
 [28] Trinkle, J.C., Berard, S., Pang, J.S.: A timestepping scheme for quasistatic multibody systems. In: (ISATP 2005). The 6th IEEE International Symposium on Assembly and Task Planning: From Nano to Macro Assembly and Manufacturing, 2005. pp. 174–181 (July 2005)
 [29] Whitney, D.E.: Force feedback control of manipulator fine motions. Journal of Dynamic Systems, Measurement, and Control 99(2), 91 (1977)
 [30] Zhou, J., Paolini, R., Johnson, A.M., Bagnell, J.A., Mason, M.T.: A probabilistic planning framework for planar grasping under uncertainty. IEEE Robotics and Automation Letters 2(4), 2111–2118 (Oct 2017)
 [31] Zhou, J., Mason, M.T., Paolini, R., Bagnell, D.: A convex polynomial model for planar sliding mechanics: theory, application, and experimental validation. The International Journal of Robotics Research 37(23), 249–265 (2018)