Modeling and Analysis of Non-unique Behaviors in Multiple Frictional Impacts

TL;DR

This paper introduces a differential inclusion model for resolving simultaneous frictional impacts in robotics, capturing multiple potential outcomes and guaranteeing termination, addressing a key challenge in impact modeling.

Contribution

It extends Routh's method to multiple contacts and provides the first termination guarantee for set-valued outcomes in simultaneous impacts.

Findings

Model captures multiple impact outcomes accurately.

Proves solutions to the impact model always terminate.

Addresses a fundamental challenge in impact dynamics modeling.

Abstract

Many fundamental challenges in robotics, based in manipulation or locomotion, require making and breaking contact with the environment. To represent the complexity of frictional contact events, impulsive impact models are especially popular, as they often lead to mathematically and computationally tractable approaches. However, when two or more impacts occur simultaneously, the precise sequencing of impact forces is generally unknown, leading to the potential for multiple possible outcomes. This simultaneity is far from pathological, and occurs in many common robotics applications. In this work, we propose an approach for resolving simultaneous frictional impacts, represented as a differential inclusion. Solutions to our model, an extension to multiple contacts of Routh's method, naturally capture the set of potential post-impact velocities.We prove that solutions to the presented model must terminate. This is, to the best of our knowledge, the first such guarantee for set-valued outcomes to simultaneous frictional impacts.

Equations

$$\dot{\mathbf{v}}\in D(\mathbf{v}).$$

$$\dot{\mathbf{v}}\in -\mathrm{Unit}\left(\mathbf{v}\right),$$

$$\mathrm{Unit}\left(\mathbf{v}\right)=\{\begin{array}{ll}\left\{\hat{\mathbf{v}}\right\}& \mathbf{v}\ne 0,\\ \mathrm{cl}{B}_1& \mathbf{v}=0.\end{array}$$

$${\mathbf{s}}_{{\mathbf{v}}_0}(t)=\{\begin{array}{ll}\left({\Vert {\mathbf{v}}_0\Vert }_2-t\right){\hat{\mathbf{v}}}_0& t\le {\Vert {\mathbf{v}}_0\Vert }_2,\\ 0& t\ge {\Vert {\mathbf{v}}_0\Vert }_2.\end{array}$$

$$\mathbf{M}(\mathbf{q})\dot{\mathbf{v}}+\mathbf{C}(\mathbf{q},\mathbf{v})={\mathbf{J}}_{n,C}(\mathbf{q}{)}^T{\mathbf{\lambda }}_{n,C}+{\mathbf{J}}_{t,C}(\mathbf{q}{)}^T{\mathbf{\lambda }}_{t,C},$$

$${\mathbf{\lambda }}_{n,C}$$

$${\mathbf{\lambda }}_{t,c_j}$$

$${\mathbf{J}}_C$$

$${\mathcal{P}}_C$$

$${\mathcal{A}}_C$$

$$\mathbf{M}(\mathbf{q})({\mathbf{v}}_{+}-{\mathbf{v}}_{-})={\mathbf{J}}_C^T{\mathbf{\Lambda }}_C.$$

$$\dot{\mathbf{v}}\in D_c(\mathbf{v})=\{\begin{array}{ll}\left\{0\right\}& \mathbf{v}\in {\mathcal{A}}_c,\\ F_c(\mathbf{v})& \mathbf{v}\in {\mathcal{P}}_c,\\ \mathrm{co}\left(\left\{0\right\}\cup F_c(\mathbf{v})\right)& \mathrm{otherwise}.\end{array}$$

$$F_c\left(\mathbf{v}\right)={\mathbf{J}}_{n,c}^T-{\mathbf{\mu }}_c{\mathbf{J}}_{t,c}^T\mathrm{Unit}\left({\mathbf{J}}_{t,c}\mathbf{v}\right),$$

$$\dot{\mathbf{v}}={\mathbf{J}}_c^T{\mathbf{\lambda }}_c,{\mathbf{\lambda }}_{n,c}=1,{\mathbf{\lambda }}_c\in {\mathrm{FC}}_c(\mathbf{q},\mathbf{v}).$$

$$F_c\left(\mathbf{v}\right)=\left\{{\mathbf{J}}_{n,c}^T\right\}.$$

$$\sum _{c\in C}{\mathbf{\lambda }}_{n,c}={\Vert {\mathbf{\lambda }}_{n,C}\Vert }_1=1.$$

$$F_C\left(\mathbf{v}\right)=\mathrm{co}\left(\left\{f_c:c\in C,f_c\in F_c(\mathbf{v}),\mathbf{v}\in \mathrm{cl}{\mathcal{P}}_c\right\}\right),$$

$$D_C(\mathbf{v})=\{\begin{array}{ll}\left\{0\right\}& \mathbf{v}\in {\mathcal{A}}_C,\\ F_C\left(\mathbf{v}\right)& \mathbf{v}\in {\mathcal{P}}_C,\\ \mathrm{co}\left(\left\{0\right\}\cup F_C(\mathbf{v})\right)& \mathrm{otherwise}.\end{array}$$

$$\dot{\mathbf{v}}={\mathbf{J}}_C^T{\mathbf{\lambda }}_C,{\Vert {\mathbf{\lambda }}_{n,C}\Vert }_1=1,{\mathbf{\lambda }}_C\in {\mathrm{FC}}_C(\mathbf{q},\mathbf{v}).$$

$${\mathbf{J}}_C\mathbf{v}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left({\mathbf{R}}^T\mathbf{v}\right)$$

$$\mathbf{v}(s)\in {\mathrm{SOL}}_{D_C}\left(I\right)$$

$$\mathbf{v}\in {\mathcal{P}}_C$$

$${\mathbf{J}}_Q\mathbf{v}=0$$

$${\mathcal{F}}_{\left(m,k\right)}=\left\{Q\in {\mathcal{C}}_{\left(m,k\right)}:{\mathbf{J}}_Q\text{ full rank}\right\}.$$

$$\mathcal{N}=\left\{N\in \mathcal{C}:0\notin F_N\left({\mathcal{P}}_N\right)\right\}.$$

$${\mathbf{v}}_{ϵ}=(1+ϵ)\left[\begin{array}{c}-1\\ -ϵ\end{array}\right]\in {\mathcal{P}}_C,$$

$${\dot{\mathbf{v}}}_{ϵ}={\mathbf{J}}_C^T\left[\begin{array}{c}ϵ\\ 1\end{array}\right]\frac{1}{1+ϵ}=\left[\begin{array}{c}ϵ\\ 1\end{array}\right]\frac{1}{1+ϵ}\in D_C({\mathbf{v}}_{ϵ}),$$

$$\mathcal{D}=\left\{D\in \mathcal{C}:\exists {\alpha }_D(s),D\text{ is }{\alpha }_D(s)\text{-dissipative}\right\}.$$

$${\mathbf{J}}_{n,{\bar{w}}_i}$$

$${\mathbf{J}}_{n,{\bar{w}}_i}{\mathbf{v}}^j(s)$$

$${\mathbf{J}}_w^T{\mathbf{\lambda }}_w^j(s)$$

$$D_C(\mathbf{v})=\{\begin{array}{ll}\left\{0\right\}& {\mathcal{A}}_C=\{j\mathbf{v}>0\}\\ \left\{j\right\}& {\mathcal{P}}_C=\{j\mathbf{v}<0\}\\ \mathrm{co}\left(\left\{0,j\right\}\right)& j\mathbf{v}=0\end{array}$$

$${\alpha }_C(s)=\min \left\{\Vert j\Vert t,\frac{1}{2}\right\}.$$

Full text

I Introduction

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. Critical challenges facing the field lie in modeling, planning, and control of robots in these complex, multi-contact settings, particularly for locomotion [43] and manipulation [23].

Rigid-body models of dynamics and contact (see Stewart [38] or Brogliato [4] 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 [1, 10, 16, 22, 36, 37] and others. Recent research, using complementarity models, has also been conducted into multi-contact optimal planning [25, 27, 28] and control [15, 29]. Similar applications have been seen for manipulation (e.g. [34]), including quasi-static approaches [7, 14]. 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 [17, 18, 36, 41, 42] and others, including recent analysis of robot locomotion [30], the resolution of simultaneous impacts is dependent upon the sequence in which they are resolved. Simulation schemes to this problem (e.g. [9, 11, 19, 22, 24, 36, 41] and many others) focus on generation of a single solution via a heuristic (symmetry [22], potential energy [41], 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; essentially, the disambiguation performed by common simulation schemes is not grounded in physical principles. Even were we to be 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 a result, rather than focus on producing a single potential solution, here we consider the set of such solutions.

As the motivating examples in 3 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 impacts (e.g. [3, 6, 31, 39], and others) along with recent data-driven models [12, 20], experimental validation [13], and efforts to translate multi-contact simulated motions to real robots [40]. Comparable results for simultaneous impacts have largely focused on simulation, with the intent to produce a single, reasonable solution (e.g [1, 10, 36]), where Anitescu and Potra [1] and Drumwright and Shell [10] guarantee termination of their numerical methods. Other related work addresses specialized, restricted settings. Seghete and Murphey [35] developed a model where solutions were guaranteed to exist, but assumed that contact normal vectors are linearly independent. Burden et al. [5] studied discontinuous vector fields, with strong results and applications to robot impacts, but are similarly restricted to frictionless contact. Johnson et al. [21] treated a limited form of friction, but assumed that contact occurs only at massless limbs. For a quasi-static model, thus without impact, Halm and Posa [14] guaranteed existence of solutions for multi-contact motion.

This work extends Routh’s graphical model [31] 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 inherent lack of uniqueness. 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. In III, 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 formal result for set-valued solutions to simultaneous frictional impact.

II Background

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 $A$ as $\mathrm{int}(A)$, $\mathrm{cl}(A)$, and $\mathrm{co}(A)$. We identify the $l_{p}$-norm and unit direction of a vector $\bm{v}\in\mathbb{R}^{n}$ as $\left\lVert\bm{v}\right\rVert_{p}$ and $\widehat{\bm{v}}=\frac{\bm{v}}{\left\lVert\bm{v}\right\rVert_{2}}$, respectively. We define the open $r$-radius ball in $\mathbb{R}^{n}$ as $B_{r}$. We denote $\mathbb{R}^{n+}\subseteq\mathbb{R}^{n}$ as the vectors with strictly positive entries and define a function $f:\Omega\subseteq\mathbb{R}^{n}\to\mathrm{cl}\mathbb{R}^{+}$ to be positive definite if it is strictly positive on $\Omega\setminus\left\{\bm{0}\right\}$. For a single-valued function $f:A\to B$ and a set-valued function $D:A\to\mathbb{P}\left(B\right)$, we denote the image of $A^{\prime}\subseteq A$ under $f$ and $D$ as $f(A^{\prime})\subseteq B$ and $D(A^{\prime})\subseteq B$ respectively.

II-A Functional Analysis

The results herein are broadly derived from measure theory and functional analysis; for a thorough background, see Rudin [32, 33]. For a set $\Omega\subseteq\mathbb{R}^{n}$, we equip $\Omega$ with the standard Euclidean metric and norm, and integrals on $\Omega$ are with respect to the Lebesgue measure by default. The total time derivative $\dot{\bm{f}}(t)$ of an absolutely continuous function $\bm{f}(t)$ is taken in the Lebesgue sense (i.e. ${\bm{f}}(t)$ is the anti-derivative of $\dot{\bm{f}}(t)$, which is defined almost everywhere ($a.e.$)). Convergence of a sequence of functions $f_{n}$ to $f$ almost everywhere and uniformly are denoted $f_{n}\xrightarrow{a.e.}f$ and $f_{n}\xrightarrow{u}f$, respectively. A key result for the derivations in this work is the Arzelà-Ascoli Theorem [33]:

Theorem 1** (Arzelà-Ascoli).**

Let $\left({\bm{f}}_{n}\right)_{n\in\mathbb{N}}$ be a uniformly bounded equicontinuous sequence of $\mathbb{R}^{n}$-valued functions on some compact interval $I$. Then there exists a function $\bm{f}$ and subsequence $\left({\bm{f}}_{n_{k}}\right)_{k\in\mathbb{N}}$ such that $\bm{f}_{n_{k}}\xrightarrow{u}\bm{f}$.

II-B Differential Inclusions

The dynamics of many robots can be captured accurately with a system of ordinary differential equations (ODEs) $\dot{\bm{x}}=\bm{f}(\bm{x},\bm{u}),$ which relates $\bm{x}\in\mathbb{R}^{n}$, the state of the robot (typically some notion of position and velocity), to $\bm{u}\in\mathbb{R}^{m}$, a set of inputs (such as motor torques) that can be manipulated. However, the dynamics of rigid bodies under frictional contact present 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, unlike ODEs, allows for the derivative at each state to lie in a set of possible values

[TABLE]

As the map $D(\bm{v})$ associated with friction may not be continuous, conditions for a function $\bm{v}(t)$ to be a solution to the differential inclusion (1) are weakened from those of an ODE:

Definition 1**.**

For a compact interval $I$, $\bm{v}:I\to\mathbb{R}^{n}$ is a solution to the differential inclusion $\dot{\bm{v}}\in D(\bm{v})$ if $\bm{v}$ is absolutely continuous and $\dot{\bm{v}}(t)\in D(\bm{v}(t))$ $a.e.$ on $I$. Denote the set of such solutions as $\mathrm{SOL}_{D}\left(I\right)$.

Solutions to initial value problems for (1) are defined similarly:

Definition 2**.**

For $I=[a,b]$ compact, denote the set of functions $\bm{v}(t)\in\mathrm{SOL}_{D}\left(I\right)$ with $\bm{v}(a)=\bm{v}_{0}$ as $\mathrm{IVP}_{D}\left(\bm{v}_{0},I\right)$.

For example, consider the differential inclusion

[TABLE]

where $\mathrm{Unit}\left(\bm{v}\right)$ is the set-valued unit direction function

[TABLE]

For any compact interval $I=[0,T]$, the initial value problem $\mathrm{IVP}_{-\mathrm{Unit}}\left(\bm{v}_{0},I\right)$ admits the unique solution

[TABLE]

$\bm{s}_{\bm{v}_{0}}(t)$ is non-differentiable at $t=\left\lVert\bm{v}_{0}\right\rVert_{2}$ and thus is not a solution of any ODE. In general, non-emptiness, regularity, and closure of $\mathrm{IVP}_{D}\left(\bm{v}_{0},I\right)$ depend on the structure of $D(\bm{v})$; fortunately, solution sets for frictional dynamics are well-behaved due to their upper semi-continuous (u.s.c.) structure:

Definition 3**.**

A function $D:A\to\mathbb{P}\left(B\right)$ with values closed in $B$ is upper semi-continuous if $\forall\left({a}_{n}\right)_{n\in\mathbb{N}}\in A,\left({b}_{n}\right)_{n\in\mathbb{N}}\in B$ with $a_{n}\to a$, $b_{n}\to b$, and $b_{n}\in D(a_{n})$, we have $b\in D(a)$.

Proposition 1** (Aubin and Cellina [2]).**

Let $\bm{v}_{0}\in\mathbb{R}^{n}$ and $I$ be a compact interval. If $D(\bm{v})$ is uniformly bounded; u.s.c.; and closed, convex, and non-empty at all $\bm{v}$, $\mathrm{IVP}_{D}\left(\bm{v}_{0},I\right)$ is u.s.c. in $\bm{v}_{0}$. Furthermore $\mathrm{SOL}_{D}\left(I\right)$ as well as $\mathrm{IVP}_{D}\left(\bm{v}_{0},I\right)$ are non-empty and closed under uniform convergence.

Intuitively, a map is u.s.c. if its value at each $\bm{v}$ is not significantly smaller than its value at any $\bm{v}^{\prime}$ near $\bm{v}$. $\mathrm{Unit}\left(\bm{v}\right)$, for example, obeys all requirements of Proposition 1. As it is a singleton, $\mathrm{IVP}_{-\mathrm{Unit}}\left(\bm{v}_{0},I\right)$ is closed, non-empty, and convex; furthermore, if $\bm{v}_{n}\to\bm{v}_{\infty}$, then $\bm{s}_{\bm{v}_{n}}\xrightarrow{u}\bm{s}_{\bm{v}_{\infty}}$ with $\bm{s}_{\bm{v}_{\infty}}\in\mathrm{IVP}_{-\mathrm{Unit}}\left(\bm{v}_{\infty},I\right)$. An illustration of this system as well as the function $\mathrm{Unit}(\bm{v})$ can be found in Figure 2.

II-C Frictional Impact Dynamics

Many robots’ dynamics can be modeled as a system of rigid bodies experiencing contact at up to $m$ points (for a thorough introduction, see [38] and [4]). The state of such a system can be represented by configuration $\bm{q}(t)$ and velocities $\bm{v}(t)\in\mathbb{R}^{n}$. The continuous evolution is governed by

[TABLE]

where $\bm{M}(\bm{q})$ is the generalized inertial matrix; $\bm{C}(\bm{q},\bm{v})$ encompasses Coriolis and gravitational forces; $\bm{J}_{n,C}\in\mathbb{R}^{m\times n}$ projects the velocity $\bm{v}$ onto the contact normals; and $\bm{J}_{t,C}\in\mathbb{R}^{2k\times n}$ projects $\bm{v}$ onto the contact tangents of the $k\leq m$ frictional contacts. We identify the behavior with a set of contacts $C=\left\{c_{1},\dots,c_{m}\right\}$, and identify each contact $c_{i}$ with its related vectors: row $i$ of $\bm{J}_{n,C}$ and rows $2i-1$ and $2i$ of $\bm{J}_{t,C}$, denoted as $\bm{J}_{n,c_{i}}$ and $\bm{J}_{t,c_{i}}$, respectively. Denote the collection of potential contact sets as $\mathcal{C}$, thus $C\in\mathcal{C}$. We furthermore define $\mathcal{C}_{\left(m,k\right)}\subseteq\mathcal{C}$ to be the collection of sets of $m$ contacts of which $k\leq m$ are frictional. The world-frame contact normal and frictional forces $\bm{\lambda}_{n,C}(t)\in\mathbb{R}^{m}$ and $\bm{\lambda}_{t,C}(t)\in\mathbb{R}^{2k}$ must lie within the Coulomb friction cone $\mathrm{FC}_{C}\left(\bm{q},\bm{v}\right)$; that is, for all $i\in\left\{1,\dots,m\right\}$ and $j\in\left\{1,\dots,k\right\}$,

[TABLE]

where $\bm{\lambda}_{n,c_{i}}$ and $\bm{\lambda}_{t,c_{j}}$ are identified similarly to $\bm{J}_{n,c_{i}}$ and $\bm{J}_{t,c_{i}}$ and $\bm{\mu}_{c_{j}}>0$ is the friction coefficient for the $j$th contact. Additionally, we denote the lumped terms

[TABLE]

$\mathcal{P}_{C}$ is the set of actively penetrating velocities, where impact is guaranteed to occur. $\mathcal{A}_{C}$ are separating velocities, where no impact can occur. Note that $\mathbb{R}^{n}\setminus(\mathcal{P}_{C}\cup\mathcal{A}_{C})\neq\emptyset$, and velocities in this set may require impacts, as in Painlevé’s Paradox [38].

In this work, we focus on inelastic impulsive impacts, during which velocities change instantaneously. Letting $\bm{\Lambda}_{C}$ represent an impulse, pre- and post-impact velocities, $\bm{v}_{-}$ and $\bm{v}_{+}$ obey

[TABLE]

Coulomb friction poses challenges in computing $\bm{\Lambda}_{C}$, as an impact may cause stick-slip transitions or change in slip direction. For a single contact $C=\left\{c\right\}$, Routh [31] proposed a graphical method describing a path in velocity space (equivalently impulse space) from $\bm{v}_{-}$ to $\bm{v}_{+}$ which satisfies Coulomb friction differentially. To briefly summarize this technique,

Increase the normal impulse $\bm{\Lambda}_{n,c}$ with slope $\bm{\lambda}_{n,c}$. 2. 2.

Increment the tangential impulse $\bm{\Lambda}_{t,c}$ with slope $\bm{\lambda}_{t,c}$, satisfying to Coulomb friction, identical to (7) for the mid-impact velocity $\bar{\bm{v}}=\bm{v}_{-}+\bm{M}(\bm{q})^{-1}\bm{J}_{c}^{T}\bm{\Lambda}_{c}$, the velocity after net impulse $\bm{\Lambda}_{c}$. 3. 3.

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. (i.e. $\bm{J}_{N,c}\bar{\bm{v}}=0$) and take $\bm{v}_{+}=\bar{\bm{v}}$.

To later proceed to the multi-contact case, we observe that this process could be modeled as a u.s.c. differential inclusion:

[TABLE]

where $F_{c}(\bm{v})$ is equal to the net increment in velocity due to the “force” applied in steps 1) and 2) of Routh’s method. Since $\bm{q}$ is constant during an impact, we will apply the transformation $\bm{M}(\bm{q})^{.5}$ to $\bm{v}$ in (11), leaving

[TABLE]

where we retain the use of $\bm{v}$ for ease of notation. For any $\dot{\bm{v}}\in F_{c}(\bm{v})$, we can associate a set of forces $\bm{\lambda}_{C}$ such that

[TABLE]

Note that for a frictionless contact ($\bm{\mu}_{c}=0$), this simplifies to

[TABLE]

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, we will take $s$ to be the “simulation time” during the resolution of an impact event; we note that evolution of $s$ does not correspond to evolution of time, but rather measures the accumulation of impact impulse over an instantaneous collision. In a slight abuse of notation, and we will consider total derivatives such as $\dot{\bm{v}}(s)$ to be taken with respect to $s$. We will also denote the impulse (i.e. the integrated force) on a contact $c$ over a sub-interval $[s_{1},s_{2}]$ of an impact resolution as $\bm{\Lambda}_{c}\left(s_{1},s_{2}\right)$. Implicit in Routh’s method is an assumption that the terminal condition in step 1) will eventually be reached by any valid choice of increment on $\bm{\Lambda}_{c}$; if it is possible to get “stuck” with $\bm{J}_{n,c}\bm{v}<0$, then Routh’s method would be ill-defined and not predict a post impact state. This does not happen in the frictionless case, as $\bm{J}_{n,c}\bm{v}$ has constant positive derivative $\bm{J}_{n,c}\dot{\bm{v}}=\left\lVert\bm{J}_{n,c}\right\rVert_{2}^{2}$. 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:

Lemma 2**.**

$\exists S>0$* such that for any solution $\bm{v}(s)\in\mathrm{SOL}_{D_{c}}\left(\left[0,\left\lVert\bm{v}(0)\right\rVert_{2}S\right]\right)$ of the single frictional contact system defined in (11) and (12), $\exists s^{*}\in\left[0,\left\lVert\bm{v}(0)\right\rVert_{2}S\right]$, $\bm{J}_{n,c}\bm{v}(s^{*})\geq 0$.*

Proof.

See Appendix -A. ∎

The implication of Lemma 2 is that a priori, one can determine an $S>0$ proportional to the pre-impact velocity $\bm{v}_{-}$ such that any solution to the differential inclusion (11) on $\left[0,S\right]$ can be used to construct the post-impact velocity $\bm{v}_{+}$. 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 [8]. 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 [30] 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 ([1, 37] 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. 4(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. 4(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. [29] 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 at any relative rate:

Monotonically increase the normal impulse on each non-separating contact $c$ at rate $\bm{\lambda}_{n,c}\geq 0$ such that

[TABLE] 2. 2.

Increment the tangential impulse for each frictional contact $c$ at rate $\bm{\lambda}_{t,c}$ such that $\bm{\lambda}_{C}\in\mathrm{FC}_{C}(\bm{v})\,.$ 3. 3.

Terminate when $\bm{v}\not\in\mathcal{P}_{C}$.

We can understand the constraint (15) on $\bm{\lambda}_{C}$ 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 $c\in C$. As in the single contact case, we might instead think of the selection of a $\bm{\lambda}_{C}$ as picking an element of a set of admissible values for $\dot{\bm{v}}$. As before, we construct a u.s.c. differential inclusion to capture this behavior:

[TABLE]

[TABLE]

We denote total impulse over an interval $[s_{1},s_{2}]$, $\bm{\Lambda}_{C}\left(s_{1},s_{2}\right)$, as before. Similar to (13), one can extract $\bm{\lambda}_{C}(s)$ from a solution $\bm{v}(s)$ such that

[TABLE]

We illustrate the behavior of this model on the rimless wheel in Figure 3(b). While $\bm{v}$ remains in the third quadrant, the direction of $\dot{\bm{v}}$ is permitted to take any value in the convex cone outlined by the sold blue and dotted red solutions, including the simultaneous impact solution. This results in the velocity terminating at least one impact on the dashed orange set, after which behavior is identical to the single-contact system. The final velocities achievable, shown in solid green, are a superset of those given by sequential and simultaneous resolution.

III-C Properties

The construction of (17) is similar to that of the single contact system (11); it is furthermore equivalent when $C$ is a singleton. We now detail properties of the multi-contact system that are useful for analyzing its solution set.

III-C1 Existence and Closure

For any $C\in\mathcal{C}$, $D_{C}(\bm{v})$ is closed, uniformly bounded, and convex as it is constructed from the convex hull of a set of bounded vectors. Therefore by Proposition 1, we obtain the following:

Lemma 3**.**

For all $C\in\mathcal{C}$, velocities $\bm{v}_{0}$, and compact intervals $I$, $\mathrm{SOL}_{D_{C}}\left(I\right)$ and $\mathrm{IVP}_{D_{C}}\left(\bm{v}_{0},I\right)$ are non-empty and closed under uniform convergence.

III-C2 Homogeneity

As each $\mathcal{P}_{C}$ and $\mathcal{A}_{C}$ are conic, $F_{C}(\bm{v})$ and therefore $D_{C}(\bm{v})$ are positively homogeneous in $\bm{v}$. That is to say, $\forall k>0,\bm{v}\in\mathbb{R}^{n}$, $D_{C}(\bm{v})=D_{C}(k\bm{v})$. Positive homogeneity induces a similar property on $\mathrm{SOL}_{D_{C}}\left(I\right)$:

Lemma 4** (Solution Homogeneity).**

For all $C\in\mathcal{C}$, $k>0$, and compact intervals $I$, if $\bm{v}(s)\in\mathrm{SOL}_{D_{C}}\left(I\right)$, $k\bm{v}(\frac{s}{k})\in\mathrm{SOL}_{D_{C}}\left(kI\right)$.

III-C3 Equivalent Minimal Coordinate Systems

In light of (18), we have that $\bm{v}(s)-\bm{v}(s_{0})\in\mathrm{Range}\left(\bm{J}_{C}^{T}\right)$ for all solutions $\bm{v}(s)\in\mathrm{SOL}_{D_{C}}\left(I\right)$. It will be useful to analyze the the evolution of a minimal-coordinate representation of $\bm{v}$’s projection onto $\mathrm{Range}\left(\bm{J}_{C}^{T}\right)$. Let $\bm{R}$ be a matrix with columns that constitute an orthogonal basis of $\mathrm{Range}\left(\bm{J}_{C}^{T}\right)$. Therefore, $\bm{R}\bm{R}^{T}$ is an orthogonal projector onto $\mathrm{Range}\left(\bm{J}_{C}^{T}\right)$ and

[TABLE]

Therefore, by defining a new set of contacts $Q$ with equal size to $C$ such that $\bm{J}_{Q}=\bm{J}_{C}\bm{R}$, we have that

[TABLE]

We denote the collection of contact sets of this size that comply with the full rank condition (23) as

[TABLE]

Note that $Q\in\mathcal{F}_{\left(k,m\right)}$ does not require $\bm{J}_{Q}$ to have linearly independent rows; $\bm{J}_{Q}$ may have more rows than columns with enough contacts, and $Q\in\mathcal{F}_{\left(k,m\right)}$ would then imply that every perturbation of $\bm{v}$ would perturb at least one contact velocity.

III-C4 Energy Dissipation

A basic behavior of inelastic impacts is that they dissipate kinetic energy $K(\bm{v})=\frac{1}{2}\left\lVert\bm{v}\right\rVert_{2}^{2}$. We now examine the dissipative properties of the model, which function both as a physical realism sanity check and as a device to prove critical theoretical properties. On inspection of (6), (7) and (13), $K(\bm{v}(s))$ must be non-increasing, and furthermore, unless $\bm{v}(s)$ is constant, it will strictly decrease:

Lemma 5** (Dissipation).**

Let $C\in\mathcal{C}$, and let $I$ be a compact interval. If $\bm{v}(s)\in\mathrm{SOL}_{D_{C}}\left(I\right)$, then $\left\lVert\bm{v}(s)\right\rVert_{2}$ is non-increasing.

Theorem 6**.**

Let $C\in\mathcal{C}$, and let $I$ be a compact interval. If $\bm{v}(s)\in\mathrm{SOL}_{D_{C}}\left(I\right)$ and $\bm{v}\left(I\right)\subseteq\mathcal{P}_{C}$, $\left\lVert\bm{v}(s)\right\rVert_{2}$ constant implies $\bm{v}(s)$ constant.

Proof.

See Appendix -B. ∎

One might then wonder if $K(\bm{v})$ is strictly decreasing on $\mathcal{P}_{C}$. One necessary condition would be $\bm{0}\not\in D_{C}(\bm{v}^{*})$ for every $\bm{v}^{*}\in\mathcal{P}_{C}$, as otherwise $\bm{v}(s)=\bm{v}^{*}$ would be a solution to the differential inclusion. We will denote the collection of contacts that have this property as

[TABLE]

Critically, $\mathcal{N}$ covers most situations in robotics, including grasping and locomotion, with the notable exception being jamming between immovable surfaces. Sums-of-squares programming [26], a form of convex optimization, can be used to certify membership in $\mathcal{N}$.

Theorem 6 and Lemma 5 have the immediate implication that $K(\bm{v})$ strictly decreases on on $\mathcal{P}_{N}$ for $N\in\mathcal{N}$:

Theorem 7** (Strict Dissipation).**

Let $N\in\mathcal{N}$ and $I$ be a compact interval. If $\bm{v}(s)\in\mathrm{SOL}_{D_{N}}\left(I\right)$ and $\bm{v}(I)\subseteq\mathcal{P}_{N}$, $\left\lVert\bm{v}(s)\right\rVert_{2}$ is strictly decreasing.

IV Finite Time Termination

While solutions to the underlying differential inclusion are guaranteed to exist in the multi-contact model, we have yet to prove that they terminate, as in Routh’s single-contact method. Termination proofs for other simultaneous impact models (e.g. [1, 10, 35] and others) exist, but these approaches rely on comparatively limited impulsive behaviors, and thus cannot capture essential non-unique post-impact velocities. We now show that our model exhibits what we understand to be the most permissive guaranteed termination behavior:

Theorem 8**.**

For any pre-impact velocity $\bm{v}(0)$ for a contact set $N\in\mathcal{N}$, The differential inclusion (17) will resolve the impact by some $S$ proportional to $\left\lVert\bm{v}(0)\right\rVert_{2}$.

We will prove this claim as a consequence of kinetic energy decreasing fast enough to force termination—a significant expansion of Theorem 7. Even though $K$ must always decrease, Theorem 7 does not forbid $\frac{{\mathrm{d}}}{{\mathrm{d}}s}K(\bm{v})\to 0$. In fact, it is not possible to create an instantaneous bound $\frac{{\mathrm{d}}}{{\mathrm{d}}s}K(\bm{v})\leq-\epsilon<0$. For example, consider 2 frictionless, axis-aligned contacts $C$ such that $\bm{J}_{C}=\bm{I}_{2}$. For every $\epsilon>0$, we can pick a velocity

[TABLE]

[TABLE]

and arrive at $\dot{K}>-2\epsilon$. However as we take $\epsilon\to 0$, $\bm{v}_{\epsilon}$ converges to to boundary of $\mathcal{P}_{C}$ and thus will only be ably to sustain a small $\dot{K}$ 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 $\alpha(s)$-dissipativity:

Definition 4** ($\alpha(s)$-dissipativity).**

For a positive definite function $\alpha(s):\mathrm{cl}\mathbb{R}^{+}\rightarrow[0,1)$, the system $\dot{\bm{v}}\in D_{C}(\bm{v})$ is said to be $\boldmath{\alpha}(s)$-dissipative if for all $s>0$, for all $\bm{v}\in\mathrm{SOL}_{D_{C}}\left([0,s]\right)$ s.t. $\bm{v}\left([0,s]\right)\subseteq\mathcal{P}_{C}$, if $\left\lVert\bm{v}(0)\right\rVert_{2}=1$, $\left\lVert v(s)\right\rVert_{2}\leq 1-\alpha(s)$.

Denote the collection of contact sets with this property as

[TABLE]

Intuitively, if $K>0$ on $\mathcal{P}_{C}$ and $K$ decreases at a known nonzero rate, we can show that any trajectory $\bm{v}(s)$ of the multi-contact system will exit $\mathcal{P}_{C}$ at a time linearly bounded in $\left\lVert\bm{v}(0)\right\rVert_{2}$:

Lemma 9** (Bounded Exit).**

Let $\alpha_{C}(s):\mathrm{cl}\mathbb{R}^{+}\rightarrow[0,1)$ be positive definite and let $C\in\mathcal{C}$ be $\alpha_{C}(s)$-dissipative. Then $\forall S>0$, $\forall\bm{v}(s)\in\mathrm{SOL}_{D_{C}}\left(\left[0,\left\lVert\bm{v}(0)\right\rVert_{2}\frac{S}{\alpha_{C}(S)}\right]\right)$, $\bm{v}\left(\left[0,\left\lVert\bm{v}(0)\right\rVert_{2}\frac{S}{\alpha_{C}(S)}\right]\right)\not\subseteq\mathcal{P}_{C}$.

Proof.

See Appendix -C. ∎

Any contact set $C$ that complies with the strong assumption of $\alpha(s)$-dissipativity is an element of $\mathcal{N}$, as otherwise $\bm{v}$ and $K$ could be constant (i.e. $\mathcal{D}\subseteq\mathcal{N}$). Far more useful is that we will show Theorem 8 arises from the converse: that every $C\in\mathcal{N}$ exhibits $\alpha(s)$-dissipativity. This is particularly surprising for systems $C$ with $\bm{J}_{C}$ not full rank, as $\bm{v}(0)$ could be large, yet $\bm{v}_{R}(s)$, the projection of $\bm{v}(s)$ onto $\mathrm{Range}\left(\bm{J}_{C}^{T}\right)$, could be arbitrary small, permitting small $\dot{K}$. We observe that the rank of $\bm{J}_{C}$ does not effect whether or not $C\in\mathcal{D}$, as all solutions will fall into two categories: either $\bm{v}_{R}(s)$ is large, or the related minimal coordinate system will exit $\mathcal{P}_{C}$ very quickly:

Theorem 10**.**

$\forall m>0,m\geq k\geq 0$, $\mathcal{N}\cap\mathcal{F}_{\left(m,k\right)}\subseteq\mathcal{D}\implies\mathcal{N}\cap\mathcal{C}_{\left(m,k\right)}\subseteq\mathcal{D}$

Proof.

See Appendix -D. ∎

Finally, we prove the primary claim of this work. Intuitively, if there exists $C\in\mathcal{N}$ that is not $\alpha(s)$-dissipative, then one could construct a sequence of convergent solutions to $\dot{\bm{v}}\in D_{C}$ that dissipate arbitrarily small amounts of energy. Therefore their limit, also a solution to $\dot{\bm{v}}\in D_{C}$ as the solution set is closed, dissipates no energy—leading to a contradiction with Theorem 7. This argument will be used in an inductive manner, incrementing the size of the contact sets:

Theorem 11** (Dissipation Inductive Step).**

Assume $\mathcal{N}\cap\mathcal{C}_{\left(m^{\prime},k^{\prime}\right)}\subseteq\mathcal{D}$ for all $m^{\prime}\geq k^{\prime}\geq 0$ with $k^{\prime}<k$ or $k^{\prime}=k$ and $m^{\prime}<m$. Then $\mathcal{N}\cap\mathcal{C}_{\left(m,k\right)}\subseteq\mathcal{D}$.

Proof.

Suppose not. Then by Theorem 10, there is a set of contacts $C\in\mathcal{N}\cap\mathcal{F}_{\left(m,k\right)}$, $S>0$, and a corresponding sequence of solutions $\left(\bm{v}^{j}(s)\right)_{j\in\mathbb{N}}$, $\bm{v}^{j}(s)\in\mathrm{SOL}_{D_{C}}\left([0,S]\right)$, all starting with velocity magnitude $1$ ($\left\lVert\bm{v}^{j}(0)\right\rVert_{2}=1$) and never exiting $\mathcal{P}_{C}$. We must also have that each dissipates less energy than the last: $\left\lVert\bm{v}^{j}(s)\right\rVert_{2}>1-\frac{1}{j}$. As $D_{C}$ is uniformly bounded, $\bm{v}^{j}$ are unformly bounded and equicontinuous. By Theorem 1 and Lemma 3, we may assume that $\exists\bm{v}^{\infty}(s)\in\mathrm{SOL}_{D_{C}}\left([0,S]\right)$ such that $\bm{v}^{j}\xrightarrow{u}\bm{v}^{\infty}$. Therefore $\left\lVert\bm{v}^{\infty}(s)\right\rVert_{2}=1$ for all $s$ and by Theorem 6 $\bm{v}^{\infty}$ is constant. As $C\in\mathcal{N}$, by Theorem 7, $\bm{v}^{\infty}$ is not an element of $\mathcal{P}_{C}$ (i.e., $\bm{J}_{n,C}\bm{v}^{\infty}\geq\bm{0}$). As $\bm{J}_{C}$ is full rank, $\bm{J}_{C}\bm{v}^{\infty}\neq\bm{0}$. Let $\bm{\lambda}_{C}^{j}(s)$ be the corresponding force vector for each $\bm{v}^{j}(s)$.

Case 1: One contact has strictly deactivated ($\exists c\in C$, $\bm{J}_{n,c}\bm{v}^{\infty}>\bm{0}$). But then as $\bm{v}^{j}\xrightarrow{u}\bm{v}^{\infty}$, by taking a subsequence starting from sufficiently high $j$ we may assume that $c$ never activates ($\forall j,t,\bm{J}_{n,c}\bm{v}^{j}(s)>\bm{0}$), and therefore at least one of the other contacts is always active ($\bm{v}^{j}([0,S])\subseteq\mathcal{P}_{C\setminus\left\{c\right\}}$). But then only the forces from $C\setminus\left\{c\right\}$ determine $\dot{\bm{v}}^{j}$, and thus $\bm{v}^{j}\in\mathrm{SOL}_{D_{C\setminus\left\{c\right\}}}\left([0,S]\right)$. As removing a contact shrinks the set of possible forces to apply ($F_{C\setminus\{c\}}\subseteq F_{C}$), $C\setminus\left\{c\right\}\in\mathcal{N}$ and contains $m-1$ contacts. Then by assumption, for some $\alpha(s)$, $C\setminus\left\{c\right\}$ is $\alpha(s)$-dissipative. But $\left\lVert\bm{v}^{j}(s)\right\rVert_{2}\to 1$. Contradiction!

Case 2: At least one contact always slides $\left(\bm{J}_{t,C}\bm{v}^{\infty}\neq\bm{0},\,\bm{J}_{n,C}\bm{v}^{\infty}=\bm{0}\right)$. Let $W=\left\{w\in C:\left\lVert\bm{J}_{t,w}\bm{v}^{\infty}\right\rVert_{2}>\bm{0}\right\}\neq\emptyset$ be the set of contacts that slide for velocity $\bm{v}^{\infty}$. Then as $\mathrm{Unit}$ is u.s.c., $\forall w\in W,\,\mathrm{Unit}\left(\bm{J}_{t,w}\bm{v}^{j}\right)\xrightarrow{u}\mathrm{Unit}\left(\bm{J}_{t,w}\bm{v}^{\infty}\right)$ (i.e. convergence of the velocity to $\bm{v}^{\infty}$ implies convergence of the direction of sliding on each contact in $W$). Therefore WLOG by taking a subsequence starting from sufficiently high $j$ we may assume $\forall w\in W,\,\exists\bm{d}_{1,w},\bm{d}_{2,w},\bm{d}_{3,w}$ sufficiently close to $-\bm{\mu}_{w}\widehat{\bm{J}_{t,w}\bm{v}^{\infty}}$ and associated new contacts $\bar{w}_{1},\bar{w}_{2},\bar{w}_{3}$ such that

[TABLE]

for $i\in\left\{1,2,3\right\}$. Denote $\bar{W}=\bigcup_{i,w}\bar{w}_{i}$ and $\bar{C}=(C\cup\bar{W})\setminus W$. (30) and (31) in conjunction imply that, for velocities $\bm{v}^{j}(s)\approx\bm{v}^{\infty}$, 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 ($\bm{v}^{j}\in\mathrm{SOL}_{D_{\bar{C}}}\left([0,S]\right)$) which never deactivate ($\bm{v}^{j}([0,S])\subseteq\mathcal{P}_{\bar{C}}$). Figure 5 illustrates this construction. As $\bar{C}$ has strictly fewer frictional contacts than $C$ and is not $\alpha(s)$-dissipative ($\left\lVert\bm{v}^{j}(s)\right\rVert_{2}\to 1$), by assumption we must have that $\bar{C}\not\in\mathcal{N}$. By definition of $\mathcal{N}$ there must exist some penetrating velocity $\bm{v}\in\mathcal{P}_{\bar{C}}$ such that $\bm{0}\in F_{\bar{C}}(\bm{v})$ is a permissible net force. We therefore must be able to find individual contact forces $\bm{\lambda}_{n,\bar{c}}f_{\bar{c}}$ with $\bm{\lambda}_{n,\bar{c}}\geq 0$ and $f_{\bar{c}}\in F_{\bar{c}}\left(\bm{v}\right)$ for each contact $\bar{c}\in\bar{C}$ such that $\sum_{\bar{c}\in\bar{C}}\bm{\lambda}_{n,\bar{C}}=1$ and $\sum_{\bar{c}\in\bar{C}}\bm{\lambda}_{n,\bar{c}}f_{\bar{c}}=\bm{0}$. As no combination of the original contacts $C$ can create zero net force alone, one of the $\bar{w}\in\bar{W}$ must strictly activate ($\bm{\lambda}_{n,\bar{w}}f_{\bar{w}}\neq 0$). By construction of $W$ and $\bar{W}$ and the assumption of Case 2, we have $\bm{J}_{C\setminus W}\bm{v}^{\infty}=\bm{0}$, and thus $f_{\bar{c}}^{T}\bm{v}^{\infty}=0$ for each $\bar{c}\in C\setminus W$ and $f_{\bar{w}}^{T}\bm{v}^{\infty}<0$ for each $\bar{w}\in\bar{W}$. Thus $\sum_{\bar{c}\in\bar{C}}\bm{\lambda}_{n,\bar{c}}f_{\bar{c}}^{T}\bm{v}^{\infty}<0$. But then $\sum_{\bar{c}\in\bar{C}}\bm{\lambda}_{n,\bar{c}}f_{\bar{c}}\neq\bm{0}$. Contradiction! ∎

We are now ready to prove the main result of this section.

Proof of Theorem 8.

We will reach the claim by showing $\mathcal{N}=\mathcal{D}$. $\mathcal{N}\supseteq\mathcal{D}$ trivially. Any $C\in\mathcal{F}_{\left(1,0\right)}$ is of the form

[TABLE]

with $\bm{v},j\in\mathbb{R}$, $j\neq 0$. Such a system is $\alpha_{C}$-dissipative with

[TABLE]

$\mathcal{C}_{\left(1,0\right)}\subseteq\mathcal{D}$ by Theorem 10. $\mathcal{N}\subseteq\mathcal{D}$ follows from nested induction on $(m,k)$ via Theorem 11. Therefore, $\mathcal{N}=\mathcal{D}$. ∎

V Conclusion

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 the post-impact velocity set. Constructing approximate solutions to the differential inclusion poses significant challenges associated with discontinuities in $\dot{\bm{v}}$. While simple Euler schemes will converge to the true solution set [2], the convergence rate is unknown, and simulation time and therefore computational complexity would scale linearly with the scale of $\left\lVert\bm{v}_{-}\right\rVert_{2}$ given Theorem 8. Tools from time-stepping schemes (e.g. [1, 37]) may circumvent these issues. Another strategy is to precompute a formula for the entire post-impact set as a function of $\bm{v}_{-}$. 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 Lemma 2

Let $\bm{R}$ be a matrix with columns that constitute an orthogonal basis of $\mathrm{Range}\left(\bm{J}_{c}^{T}\right)$. By equivalence of norms there exists $\epsilon>0$ such that

[TABLE]

Pick $S=\left(\epsilon\min\left(\bm{\mu}_{c},1\right)\right)^{-1}$. Let $V(s)=\left\lVert\bm{R}^{T}\bm{v}(s)\right\rVert_{2}^{2}$. Assume $\bm{v}(s)\in\mathcal{P}_{c}$ for $s<s^{*}=\left\lVert\bm{R}^{T}\bm{v}(0)\right\rVert_{2}S\leq\left\lVert\bm{v}(0)\right\rVert_{2}S$.

[TABLE]

on $[0,s^{*}]$ and thus $V\left(s^{*}\right)\leq\left(\sqrt{V(0)}-\epsilon\min\left(\bm{\mu}_{c},1\right)s^{*}\right)^{2}=0$. Therefore $\bm{J}_{n,c}\bm{v}\left(s^{*}\right)=0$.

-B Proof of Theorem 6

Let $\bm{v}(s)\in\mathrm{SOL}_{D_{C}}\left(I\right)$ with $\bm{v}(s)$ non-constant. Let $\bm{\lambda}_{C}(s)$ be the associated vector of force variables. As $\bm{v}(s)$ is continuous, we may select $s^{*}\in\mathrm{int}I$ such that $\forall\delta>0$, $\bm{v}(s)$ is non-constant on $[s^{*},s^{*}+\delta]$. Let $A=\left\{a\in S:\bm{J}_{n,a}\bm{v}(s^{*})\leq 0\right\}$ be the set of active contacts at $s=s^{*}$. Let $B$ the the largest subset of $A$ such that $\bm{J}_{n,B}\bm{v}=\bm{0}$ and $\bm{J}_{t,B}\bm{v}=\bm{0}$. As $\bm{v}$ is continuous, $\exists\delta_{\epsilon}>0$ and $\epsilon>0$ such that $\forall s\in[s^{*},s^{*}+\delta_{\epsilon}]\subseteq I$,

• •

$\bm{J}_{n,C\setminus A}\bm{v}(s)>\epsilon$

• •

$\bm{J}_{n,c}\bm{v}(s)<-\epsilon$ for $c\in A\setminus B$ frictionless

• •

$\bm{J}_{n,c}\bm{v}(s)<-\epsilon$ or $\left\lVert\bm{J}_{t,c}\bm{v}(s)\right\rVert_{2}>\frac{1}{\bm{\mu}_{i}}\epsilon$ for $c\in A\setminus B$ frictional.

Therefore no new contacts activate before $s^{*}+\delta_{\epsilon}$, and

[TABLE]

on $[s^{*},s^{*}+\delta_{\epsilon}]$. Select one such $s$ with $\bm{v}(s)\neq\bm{v}(s^{*})$. By Lemma 5,

[TABLE]

Therefore, we must have $\left\lVert\bm{\Lambda}_{A\setminus B}(s^{*},t)\right\rVert_{1}>0$. Finally,

[TABLE]

Therefore $\left\lVert\bm{v}\right\rVert_{2}$ is non-constant.

-C Proof of Lemma 9

Assume WLOG by Lemma 4 that $\left\lVert\bm{v}(0)\right\rVert_{2}=1$ and that $\bm{v}(s)\in\mathcal{P}_{C}$ on $0\leq s<\frac{S}{\alpha_{C}(S)}$. As $C$ is $\alpha_{C}(s)$-dissipative, $\exists s_{1}\in\left[0,S\right]$ such that $\left\lVert\bm{v}\left(s_{1}\right)\right\rVert_{2}\leq 1-\alpha_{C}(S)$. A sequence $\left({s}_{k}\right)_{k\in\mathbb{N}}$ can be iteratively constructed by Lemma 4 such that

• •

$s_{k}\in\left[s_{k-1},s_{k-1}+S\left(1-\alpha_{C}(S)\right)^{k-1}\right]\subseteq\left[0,\frac{S}{\alpha_{C}(S)}\right]$

• •

$\left\lVert\bm{v}\left(s_{k}\right)\right\rVert_{2}\leq\left(1-\alpha_{C}(S)\right)\left\lVert\bm{v}\left(s_{k-1}\right)\right\rVert_{2}\leq\left(1-\alpha_{C}(S)\right)^{k}$

Therefore $\exists s_{\infty}\in\left[0,\frac{S}{\alpha_{C}(S)}\right]$ with $s_{n}\to s_{\infty}$ and by continuity of $\bm{v}$, $\bm{v}\left(s_{\infty}\right)=\bm{0}\not\in\mathcal{P}_{C}$.

-D Proof of Theorem 10

Let $C\in\mathcal{N}\cap\mathcal{C}_{\left(m,k\right)}$. Let $\bm{R}$ and $\bm{N}$ be matrices with columns that constitute orthogonal bases of $\mathrm{Range}\left(\bm{J}_{C}^{T}\right)$ and $\mathrm{Null}\left(\bm{J}_{C}\right)$, respectively. Therefore there exists contact set $Q$ of size $(m,k)$ and a positive definite function $\alpha_{Q}(s)$ such that $\bm{J}_{Q}=\bm{J}_{C}\bm{R}$ is full column rank, $\mathcal{P}_{Q}=\bm{R}^{T}\mathcal{P}_{C}$, and $Q$ is $\alpha_{Q}(s)$-dissipative. Let $s>0$, $\bm{v}\in\mathrm{SOL}_{D_{C}}\left([0,s]\right)$, $\left\lVert\bm{v}(0)\right\rVert_{2}=1$, and $\bm{v}\left(\left[0,s\right]\right)\subseteq\mathcal{P}_{C}$. Decompose $\bm{v}(s)=\bm{v}_{R}(s)+\bm{v}_{N}(s)=\bm{RR}^{T}\bm{v}(s)+\bm{NN}^{T}\bm{v}(0)$. We must have $\bm{R}^{T}\bm{v}\in\mathrm{SOL}_{D_{Q}}\left([0,s]\right)$. Therefore as $\bm{R}^{T}\bm{v}([0,s])\subseteq\bm{R}^{T}\mathcal{P}_{C}=\mathcal{P}_{Q}$, by Lemma 9, $s<\left\lVert\bm{R}^{T}\bm{v}(0)\right\rVert_{2}\frac{s}{\alpha_{Q}(s)}$. Thus $\left\lVert\bm{v}_{R}(0)\right\rVert_{2}>\alpha_{Q}(s)\left\lVert\bm{v}(0)\right\rVert_{2}$ and

[TABLE]

Therefore $C$ is $\left(1-\sqrt{1-\alpha_{Q}^{3}\left(1-\alpha_{Q}\right)}\right)$-dissipative.