Modern robots are fast and strong, and, in some situations, their capabilities eclipse those of humans. However, when these robots interact with their environment, whether by manipulating objects or traversing over uneven surfaces, they do so with far less skill than their human counterparts. The speed and accuracy that robots display in static or structured settings disappears when faced with less structured tasks. Critical challenges facing the field lie in modeling, planning, and control of robots in these complex, multi-contact settings, particularly for locomotion  and manipulation .
Rigid-body models of dynamics and contact (see Stewart  or Brogliato  for an overview) are widely used in robotics, as they can lead to far more tractable methods than approaches which explicitly attempt to capture the stiff interaction between objects. These approaches have also led to complementarity-based simulation schemes, such as [32, 1, 17, 13] and others. Recent research, using complementarity models, has also been conducted into multi-contact optimal planning [22, 20, 23] and control [24, 12]. Similar applications have been seen for manipulation (e.g. ), including quasi-static approaches [7, 11]. When impacts occur, rigid-body models approximate the event as an instantaneous change in velocity due to an impulsive force.
The approaches above, now deeply ingrained within the robotics community, universally assume that it is possible to determine a single potential post-impact velocity, even during simultaneous multi-contact. However, as observed in [37, 31, 36] and others, including recent analysis of robot locomotion , the resolution of simultaneous impacts is dependent upon the sequence in which they are resolved. Simulation schemes to this problem (e.g. [31, 36, 17, 19, 14]
) focus on generation of a single solution via a heuristic (symmetry, potential energy , etc.). However, for many practical applications in robotics, it is not possible to create a model detailed enough to reliably disambiguate between the multiple potential solutions. Furthermore, even given such detail, this lack of uniqueness often represents an extreme sensitivity to initial conditions: slight perturbations in the initial state of the system might lead to different impact sequences.
As the motivating examples in III-A will demonstrate, simultaneous impacts are not limited to unlikely, pathological events but are, in fact, regular occurrences in robotics and require careful analysis. From the perspective of planning, learning and control, it is critical to understand the role of this non-uniqueness (alternatively, extreme sensitivity), as some of the broad challenges in executing dynamic, multi-contact motion likely arise from these issues. For example, methods which use a simulator to learn or plan a motion may, unwittingly, be planning for an ambiguous, therefore unstable, outcome due to multi-contact. Furthermore, as the set of these ambiguous outcomes is often non-convex, it is insufficient to try to capture this sensitivity via simple models of uncertainty.
Many methods have been proposed for modeling single, frictional impacts (e.g. [26, 34, 6, 3], and others) along with recent data-driven models [9, 15], experimental validation , and efforts to translate multi-contact simulated motions to real robots . However, comparable results for simultaneous impacts have been limited to restricted models; addressing the complexity of Coulomb friction has been particularly difficult. Seghete and Murphey 
presented a model where solutions were guaranteed to exist, but assumed that impacts are frictionless and that contact normal vectors are linearly independent.Burden et al.  studied discontinuous vector fields, with strong results and applications to robot impacts, but are similarly restricted to frictionless contact. Johnson et al.  treated a limited form of friction, but assumed that the effects of contact forces are inertially decoupled. For a quasi-static model, thus without impact, Halm and Posa  guaranteed existence of solutions for multi-contact motion.
This work extends Routh’s original graphical model 
to address simultaneous, inelastic impacts by permitting impulses to occur in arbitrary sequences. As a result, the model produces a set-valued map that captures the lack of uniqueness inherent in the problem. We believe this is the appropriate description for robotic planning and control, as motions that present as non-unique will, for physical systems, display extreme sensitivity to any errors in estimation or control. In contrast with prior literature, the presented model captures a broad class of frictional systems. InIII, we describe the model and a number of its theoretical properties and in IV we prove the key result that the impact model is guaranteed to terminate. To the best of the authors’ knowledge, this work presents the first known proof of existence of solutions for simultaneous frictional impact.
We now introduce notation for and study the limiting behaviors of the frictional impact dynamics of rigid multibody systems. Denote the interior, closure, and convex hull of a set as , , and . We identify the -norm and unit direction of a vector as and , respectively. We define the open -radius ball in as . We denote as the vectors with strictly positive entries and definite a function to be positive definite if it is strictly positive on . For a single-valued function and a set-valued function , we denote the image of under and as and respectively.
Ii-a Functional Analysis
The results herein are broadly derived from measure theory and functional analysis; for a thorough background, see Rudin [27, 28]. For a set , unless otherwise specified we equip with the standard Euclidean metric and norm, and integrals on are taken with respect to the Lebesgue measure. We denote the -valued Lebesgue spaces and Sobolev space of degree one on a compact interval as and respectively, and particularly note that is a Hilbert space. The total time derivative of an absolutely continuous functions is taken in the Lebesgue sense and defined almost everywhere () such that
Convergence of sequences of functions will be frequently used. Convergence of a sequence of functions to almost everywhere, uniformly, and with respect to the norm of a Hilbert space are denoted , , and respectively. Furthermore, weak convergence in is denoted as . It is well know that for any compact interval , is compactly embedded in , leading to a key result for the derivations in this work:
[Rellich] Let be a bounded sequence in for some compact interval . if (weakly), then (strongly).
Ii-B Differential Inclusions
The dynamics of many robots can be captured accurately with a system of ordinary differential equations (ODEs)which relates , the state of the robot (typically a generalized notion of its position and velocity), to , a set of inputs (such as motor torques) that can be manipulated. However, the dynamics of rigid bodies under frictional contact present two additional complexities that this formulation cannot capture. Impacts between bodies induce instantaneous jumps in velocity that in general cannot described by an ODE (non-smooth behaviors). Additionally, when contact occurs at many points, multiple frictional forces that obey Coulomb’s laws of friction may exist (non-unique behaviors). It is therefore useful to define an object that is more permissive than ODEs by allowing for the derivative at each state to lie in a set of possible values
As the derivative map associated with Coulomb friction may not be differentiable or even continuous in , conditions for a function to be considered a solution to the differential inclusion (2) are weakened from those of an ODE: For a compact interval , is a solution to the differential inclusion if is absolutely continuous and on . Denote the set of such solutions as . Solutions to initial value problems for (2) are defined in a similar manner: For a compact interval , denote the set of functions that are solutions to the initial value problem as . For example, consider the differential inclusion
where is the set-valued unit direction function
For any compact interval , the initial value problem admits the unique solution
We note that is differentiable except at , thus no ODE would admit as a solution. In general, non-emptiness, regularity, and closure of are highly dependent on the structure of ; fortunately, frictional dynamics admit an upper semi-continuous (u.s.c.) structure that moderates the behavior of their solution sets: A function with values closed in is upper semi-continuous if with , , and , we have . [Aubin and Cellina ] Let and be a compact interval. If is uniformly bounded; u.s.c.; and closed, convex, and non-empty at all , is u.s.c. in . Furthermore as well as are non-empty and closed under uniform convergence.
Intuitively, a map is u.s.c. if its value at each is not significantly smaller than its value at any near . for example obeys all requirements of Proposition 5. We have already seen that is closed, non-empty, and convex. It is easy to see that if , then with . An illustration of this system as well as the function can be found in Figure 2.
Ii-C Frictional Impact Dynamics
The dynamics of robotic systems can be modeled as a system of rigid bodies experiencing frictional contact at up to points (for a thorough introduction, see  and ). The state of such a system can be represented by generalized coordinates and velocities . The continuous evolution is governed manipulator equations
where is the generalized inertial matrix; encompasses Coriolis and gravitational forces; projects the velocity onto the contact normals; and projects onto to the contact tangents of the frictional contacts. We identify the behavior with a set of contacts , and identify each contact with it’s related vectors: row of and rows and of , denoted as and , respectively. Denote the collection of potential contact sets as , thus . We furthermore denote , the collection of sets of contacts of which are frictional. The world-frame contact normal and frictional forces and are constrained to lie within the Coulomb friction cone ; that is, for all and ,
where and are identified similarly to and and is the friction coefficient for the th contact. Additionally, we denote the lumped terms
Intuitively, is the set of actively penetrating velocities, where impact is guaranteed to occur. are separating velocities, where no impact can occur. Note that , and velocities in this set may require impacts, as in Painlevé’s Paradox .
In this work, we focus on inelastic impulsive impacts. During such an impact, the velocity changes instantaneously. Letting represent an impulse, the relationship between pre- and post-impact velocities, and , is
Coulomb friction poses challenges in computing , as an impact may cause stick-slip transitions or the direction of slip to change. For a single contact , Routh  proposed a graphical method describing a path in velocity space (equivalently impulse space), from pre- to post-impact velocities, where the path must differentially satisfy Coulomb friction. To briefly summarize Routh’s technique,
Increase the normal impulse with some slope .
Increment the tangential impulse with slope , satisfying to Coulomb friction, identical to (8) and depending on , the velocity after net impulse .
Terminate when the normal contact velocity vanishes111To permit resolutions to Painlevé’s Paradox, terminate only when consistency no longer requires an instantaneous change in velocity., , and take .
To later proceed to the multi-contact case, we observe that this process could be modeled as an upper semi-continuous differential inclusion:
where we retain the use of for ease of notation. For any , we can associate a set of forces such that
Note that for a frictionless contact (), this simplifies to
A diagram depicting the resolution of a potential planar impacts is shown in Figure 2. Solutions may transition between sliding and sticking, and the direction of slip may even reverse as a result of each impact. While the path is piecewise linear in the planar case, this is not true in three dimensions.
From this point forward, in a slight abuse of notation, we will take to be the “time” during the resolution of an impact event. We will also define the net impulse over a sub-interval of an impact resolution:
Implicit in Routh’s method is an assumption that the terminal condition in step 1) will be eventually be reached by any valid choice of increment on ; if it is possible to get “stuck” with , then Routh’s method would be ill-defined and not predict a post impact state. It is easy to see that this does not happen in the frictionless case, as has constant positive derivative . The frictional case requires more careful treatment. Intuitively, the added effect of the frictional impulse will be to dissipate kinetic energy quickly. One may conclude that termination happens eventually as zero velocity is a valid post-impact state:
Let be a matrix with columns that constitute an orthogonal basis of . By equivalence of norms there exists such that
Pick . Let . Assume for .
on and thus . Therefore . ∎
The implications of Lemma II-C is that a priori, one can determine a proportional to the pre-impact velocity such that any solution to the differential inclusion (12) on can be used to construct the post-impact velocity . We will see, however, that the extension of this methodology to multiple concurrent impacts is non-trivial, and that the physicals systems associated with these models often exhibit a high degree of indeterminacy.
Iii Simultaneous Impact Model
Iii-a Motivating Examples
We include, as motivation, two common robotics examples that exhibit simultaneous impacts: one related to legged locomotion and the other to manipulation. Both examples, depending on initial conditions and model properties, can exhibit non-uniqueness. Before describing our model in full detail, we present these examples by considering the outcome of applying Routh’s method to a single contact at a time.
Iii-A1 Rimless Wheel
The rimless wheel is a commonly used description of simple robotic walking . Here, we will analyze the case where two feet contact the ground. This can occur if the robot were to fall on two feet, simultaneously, or when one foot is in sustained ground contact and the other impacts the ground. Note that this example is not limited to a legged robot with locked hip and knee joints; see  for a thorough analysis of similar legged examples.
For a simple example, illustrated in Fig. 3, we assume that both feet strike the ground vertically, with friction sufficient to sustain sticking. In this case, existing simulation schemes ([32, 1] and others) predict that equal impulses are generated on both feet, brining the robot to rest immediately. However, as illustrated in the figure, if the contacts are sequenced one at a time, other post-impact states are possible where one leg separates from the ground. For other configurations of this problem, non-unique solutions exist spanning sticking, sliding, and separation all for a single initial condition.
Iii-A2 Nonprehensile Pushing
In this second example, motivated by nonprehensile pushing of an object, we take a box-like object (Fig. 4) to have one corner sliding along a surface before impacting a frictionless second surface. Here, the impact on the right wall causes a secondary, frictional impact against the lower wall.
If the first impact is taken to termination before activating the contact on the right wall, the solution in Fig. 3(b) is discovered. Here, the bottom contact is separating and the right contact is sliding upward. Instead, if the impact switches prior to termination, shown in Fig. 3(c), a slightly different solution emerges. This example illustrates that, in simple cases, reminiscent of common robotics applications, subtly different non-unique solutions can emerge from multiple contacts.
Iii-B Model Construction
As post-impact velocity is sensitive to the ordering of individual impact resolutions, if we would like to predict as many reasonable post-impact velocities as possible, we must use as relaxed of a notion of impact resolutions as possible. A similar model, without theoretical results or a detailed understanding, was proposed by Posa et al.  where it proved useful for stability analysis of robots undergoing simultaneous impact. We consider a formulation in which at any given instant during the resolution process, the impacts are allowed to concurrently resolve an any relative rate:
Monotonically increase the normal impulse on each non-separating contact at rate such that
Increment the tangential impulse for each frictional at rate such that
Terminate when .
We can understand the constraint (22) on as choosing a net force that comes from a convex combination of the forces that Routh’s method might select for any of the individual contacts . As in the single contact case, we might instead think of the selection of a as picking an element of a set of admissible values for . As before, we construct a u.s.c. differential inclusion to capture this behavior:
The construction of (24) appears quite similar to that of the single contact systems (13) and (15)—in fact, it is equivalent when is a singleton. We now detail several properties of the multi-contact system that are useful for the analysis of its solution set.
Iii-C1 Existence and Closure
For any , is trivially closed, uniformly bounded, and convex as it is constructed from the convex hull of a set of bounded vectors. Therefore, in light of Proposition 5, we obtain the following: For all , velocities , and compact intervals , and are non-empty and closed under uniform convergence.
As each and are conic, and therefore are positively homogeneous in . That is to say, , . Positive homogeneity induces a similar property on : [Solution Homogeneity] For all , , and compact intervals , if , .
Iii-C3 Equivalent Minimal Coordinate Systems
In light of (25), for a set of contacts , we have that for all solutions . As in Lemma II-C, it will be useful to analyze the the evolution of a minimal-coordinate representation of ’s projection onto . Let be a matrix with columns that constitute an orthogonal basis of . Therefore is an orthogonal projector onto and
Therefore, by defining a new set of contacts with equal size to such that , we have that
We denote the collection of contact sets of this size that comply with the full rank condition (30) as
Note that does not imply that all of the rows of are linearly independent; for systems with many contacts, will have many more rows than columns, and in this case implies that for every element of the velocity , there exists a contact that can perturb it.
Iii-C4 Energy Dissipation
A basic behvaior of inelastic impacts is that they dissipate kinetic energy . We now examine the dissipative properties of the model, which function both as a sanity check on its physical realism and as a device to prove critical theoretical properties. On inspection of (7)–(8) and (14), must be non-increasing, and furthermore, unless is stationary, it will decrease by a non-zero amount: [Dissipation] Let and . Then is non-increasing. Let , and let be a compact interval. If and , constant implies constant.
See Appendix -A. ∎
One might then wonder if is strictly decreasing on . This would obviously be untrue if there existed a with as would be an element of for all compact . We will denote the collection of contacts that do not have this property as
Critically, covers most situations in robotics, including grasping and locomotion, with the notable exception being jamming against an immovable surface. Sums-of-squares programming , a form of convex optimization, can be used to certify membership in .
Iv Finite Time Termination
While solutions to the underlying differential inclusion are guaranteed to exist in the multi-contact model, we have yet to show that they are guaranteed to terminate, as in Routh’s single-contact method. There have been termination proofs for other simultaneous impact models (e.g. ), but only under non-trivial assumptions (no friction; limited number of contacts; smaller space of solutions; etc.). We show what we understand to be the most permissive guarantee of termination: For any pre-impact velocity for a contact set , The differential inclusion (24) will resolve the impact by some time proportional to .
We will show that this claim is true as a consequence of kinetic energy decreasing fast enough to force termination—a significant expansion of the claim in Theorem III-C4. Even though must always decrease, Theorem III-C4 does not forbid . In fact, it is not possible to create an instantaneous bound . For example, consider 2 frictionless, axis-aligned contacts such that . For every , we can pick a velocity
and arrive at . However as we take , converges to to boundary of and thus will only be ably to sustain a small for a small amount of time before terminating the impact. It remains possible that the aggregate energy dissipation over an interval of fixed nonzero length can be bounded away from zero. We establish a rigorous characterization of this quality by defining -dissipativity:
[-dissipativity] For a positive definite function , the system is said to be -dissipative, if for all , for all s.t. , if , . Denote the collection of contact sets with this property as
Unsurprisingly, if on and decreases at a known negative rate, we can show that any trajectory of the multi-contact system will exit at a time linearly bounded in : [Bounded Exit] Let be positive definite and let be -dissipative. Then , , .
It is trivial to show that any contact set that complies with the strong assumption of -dissipativity is an element of , as otherwise and could be constant (i.e. ). Far more useful is that we will show Theorem IV arises from the converse: that every exhibits -dissipativity. This is particularly surprising for systems with not full rank, as could be large, yet the projection of onto could be arbitrary small, which in turn allows to be small. We observe that the rank of does not effect whether or not , as all solutions will fall into two categories: either is large, or the related minimal coordinate system will exit very quickly: ,
See Appendix -B. ∎
Finally, we prove the primary claim of this work. Intuitively, if there exists that is not -dissipative, then one could construct a sequence of convergent solutions to that dissipate arbitrarily small amounts of energy. Therefore their limit, also a solution to , as the solution set is closed, dissipates no energy—leading to a contradiction with Theorem III-C4. This argument will be used in an inductive manner, incrementing the size of the contact sets: [Dissipation Inductive Step] Assume for all with or and . Then .
Suppose not. Then by Theorem IV, there is a set of contacts , , and a corresponding sequence of solutions , , all starting with velocity magnitude () and never exiting . We must also have that each dissipates less energy than the last: . As is uniformly bounded, are equicontinuous and bounded in . By Theorem II-A, we may assume that such that . Therefore for all and by Theorem III-C4 is constant. As , by Theorem III-C4, is not an element of (i.e., ). As is full rank, . Let be the corresponding force vector for each .
Case 1: One contact has strictly deactivated (, ). But then as , by taking a subsequence starting from sufficiently high we may assume that never activates (), and therefore at least one of the other contacts is always active (). But then only the forces from determine , and thus . As the shrinking the contact set shrinks the set of possible forces to apply (), and contains only contacts. But then, by assumption, for some , is -dissipative. But . Contradiction!
Case 2: At least one contact always slides . Let be the set of contacts that slide for velocity . Then as is upper semi-continuous, ; that is to say, convergence of the velocity to implies convergence of the direction of sliding on each contact in . Therefore WLOG by taking a subsequence starting from sufficiently high we may assume sufficiently close to and associated new contacts such that
for . Denote and . (37) and (38) in conjunction imply that, for velocities that stay close to , each sliding frictional contact pushes mostly in one direction. Furthermore, the associated frictional force can be generated by three frictionless contacts tilted away from the sliding direction () which never deactivate (). An illustration for this construction is available in Figure 5. As has strictly fewer frictional contacts than and is not -dissipative (), by assumption we must have that . By definition of there must exist some penetrating velocity such that is a permissible net force. We therefore must be able to find individual contact forces with and for each contact such that and . As no combination of the original contacts can create zero net force alone, one of the must strictly activate (). By construction of and and the assumption of Case 2, we have , and thus for each and for each . Thus . But then . Contradiction! ∎
We are now ready to prove the main result of this section.
Non-unique behavior is a pervasive complexity that is present in both real-world robotic systems and common models capturing frictional impacts between rigid bodies—and thus accurate incorporation of such phenomena is an essential component of robust planning, control, and estimation algorithms. Our model presents a state-of-the-art theoretical foundation for the capture of this behavior, because despite the high versatility of allowing impacts to resolve at arbitrary relative rates, it is guaranteed to terminate in finite time under far more modest conditions than shown for previous models.
The logical progression from these theoretical results is to develop a numerical scheme to generate post-impact velocities. Constructing approximate solutions to the differential inclusion poses significant challenges associated with discontinuities in ; tools from time-stepping schemes (e.g. [1, 32]) may circumvent this issue. Another strategy is to precompute a formula for the entire post-impact set as a function of . Sums-of-squares programming presents potential for construction of an outer approximation.
Future generalizations of the model include elastic impacts using Poisson restitution; resolution of Painlevé’s Paradox; and a full rigid body dynamics model that has continuous solutions through impact.
-a Proof of Theorem Iii-C4
Let with non-constant. Let be the associated vector of force variables. As is continuous, we may select such that , is non-constant on . Let