1  A propagative model of simultaneous impact: existence, uniqueness, and design consequences  Vlad Seghete,   Student Member, IEEE,   Todd D. Murphey,   Member, IEEE  Abstract —This   paper   presents   existence   and   uniqueness results for a propagative model of simultaneous impacts that is guaranteed to conserve energy and momentum in the case of   elastic   impacts   with   extensions   to   perfectly   plastic   and inelastic   impacts.   A   corresponding   time-stepping   algorithm that guarantees conservation of continuous energy and discrete momentum is developed, also with extensions to plastic and inelastic impacts. The model is illustrated in simulation using billiard balls and a two-dimensional legged robot as examples; the latter is optimized over geometry and gait parameters to achieve unique simultaneous impacts.  Note to Practitioners —Simultaneous impacts are a common oc- currence in manufacturing and robotic applications. Simulation- based techniques predicting the motion of a mechanical system subject to simultaneous impacts often use numerical routines that make algorithmic assumptions about impact in order to make   the   simulation   more   tractable.   Such   assumptions   have the potential to significantly influence the simulation outcome and may even invalidate results. This paper provides tools for simulating   simultaneous   impacts,   verifying   a   simulation   that deals with simultaneous impacts, and designing a system so that the simulation will be less dynamically sensitive to indeterminacy in simultaneous impact.  Index Terms —Dynamics, Legged Robots, Animation and Sim- ulation, Impact Modeling  I. I NTRODUCTION  I N this paper we investigate a propagative impact model for rigid body systems [1], [2] and the conditions under which it provides a deterministic prediction for the outcome of simultaneous impacts. Our motivating system is pictured in Fig. 1 and represents a simplified model based on the geometry of numerous legged robots [3]–[7] While contact modeling has received significant attention in relation to plastic impacts and established contacts [8], [9], little work has been done to address the issue of simultaneous non-plastic impacts. Such interactions naturally occur in the system of Fig. 1—between a foot and the floor while the tail is still in contact—but also in the gait of robots such as RHex [7] or IMPASS [5]. As simulation has proven a strategic tool in the design of mobile robots [10], we present theoretical results that can be used both as building blocks for physical simulations as well as a validation tool for existing simulation methods. When simulating impacts in physical systems two main approaches have been widespread. The most popular method is to solve unilateral contact problems implicitly—usually by  V.   Seghete   and   T.   Murphey   are   with   the   Department   of   Mechanical Engineering at Northwestern University, 2145 Sheridan Rd., Evanston, IL 60660.  solving a linear complementarity problem (LCP) associated with every time step of the simulation [11]–[18]. The defining feature of this implicit approach is that all the collisions detected over a time step are processed together with the dynamics   for   that   time   step,   which   eases   implementation and scalability. In addition, implicit LCP methods potentially have very fast execution times and are parallelizable [19]– [21]. This makes such methods ideal in most real-time phys- ical simulations—e.g. interactive graphics applications, video games, etc. However, the implicit nature of such algorithms is double-edged: at the cost of a high-performing low-hassle implementation   they   can   allow   for   unrealistic   behavior— especially noticeable when dealing with non-plastic impacts. To illustrate such behavior, we simulated the three sphere system shown in Fig. 2a using the popular Bullet physics engine [21]. We let spheres   B   and   C   be in contact initially and had sphere   A   impact   B   with varying velocities. Figure 2b shows the amount of energy lost through this impact for three different values of the coefficient of restitution. For perfectly elastic collisions the system energy increased by as much as   30%   due to the impact. Furthermore, this error does not scale with the time step; while varying the time step frequency between 60 Hz and 6 kHz the error did not change, indicating we cannot expect this algorithm to converge to a “true” solution as the step size goes to zero. Aditionally, even at velocities for which the energy was conserved, the solution did not come close to the experimental results [22]. It is worth noting here that while most LCP methods are implicit, work has been done that incorporates the LCP approach to impacts into event-driven integration schemes [23], [24] as discussed below. An alternative method to modeling rigid body collisions is to treat impacts as impulsive events where a jump in velocity occurs at the moment of impact. At each impact a reset map determines the post-impact state of the system. The main benefit of this approach is the explicit nature of the reset map which allows one control of the simulation’s physical accuracy  Fig. 1: The running robot model (TREX) that we used in order to illustrate the results of this paper. The arrow depicts the intended direction of movement which results from a counter-clockwise torque applied to the legs at the hips.  arXiv:1709.02296v1 [cs.RO] 7 Sep 2017
x x A  x B  x C  Newton’s Cradle (a) Three spheres restricted to move on a line. The configuration space is three dimensional, with   x A ,   x B   and   x C   as the config- uration variables. This system is a simplification of the Newton’s cradle toy. 1   2   3   4   5 -40 0 40 80 Initial velocity   ( m/s )  Relative energy loss   (%)  COR = 1.0 COR = 0.7 COR=0.0  (b) Energy loss versus initial velocity for the simplified Newton’s cradle simulated in Bullet (dots) and using our method (dotted lines). Each band of points and dotted line represents one of three coefficients of restitution (COR) tested: 1, 0.7 and 0. Fig. 2: Simplified Newton’s cradle showing unrealistic behavior in the Bullet physics engine. Notice the inconsistencies in energy behavior present in the Bullet data points.  and compliance with energy conservation and restitution laws at the time of impact. Using an explicit method also gives us the freedom to employ a   propagative   model for handling simultaneous impacts, a model which has better agreement with experimental data [22] than the implicit methods used in   obtaining   the   results   of   Fig.   2b.   Propagative   collision models solve simultaneous impacts sequentially, splitting one simultaneous impact into as many two body impacts as needed to reach a solution. The method is based on assuming an infinitesimal gap between all contacts and solving for two- body impacts in an order of one’s choice. While there is experimental evidence that better qualitative results can be achieved by using soft-body models [25], [26] and similar approaches [27]–[29], the errors incurred by using a rigid body propagative model are slight at best [22] and certainly not to the degree of the errors presented in Fig. 2b. 1   The rigid body approach we undertake allows us to model more complex systems, like the legged robot model shown in Sec. VII-A. Un- fortunately, the main drawback in using rigid body propagative approaches, as discussed in [1] and [28], is the inadequate modeling of continuum mechanics phenomena. In particular, the coupling effects between simultaneous impacts cannot be generally addressed. In the rigid body limit such coupling effects manifest themselves as a lack of unique solutions to the equations governing simultaneous impact. However, rather than proposing to solve the general simultaneous impact problem, we focus instead on identifying cases in which coupling effects do not generate non-uniqueness in the rigid body limit for simultaneous impacts. Such cases, therefore, can be modeled under the rigid body assumption provided mechanism and gait design goals are set appropriately. As major contributions the present paper provides the fol- lowing: 1) sufficient conditions for solution existence for the prop- agative model of elastic impact, assuming two simulta-  1 Note that the results presented in [22] reflect the behavior of the system when the initial velocities have relatively high values compared to the stiffness of the bodies involved in the collisions.  neous impacts 2) sufficient   conditions   for   solution   uniqueness   for   the propagative model of elastic impact, assuming   n   simul- taneous impacts; when uniqueness is satisfied, existence is guaranteed for arbitrary numbers of simultaneous impacts. 3) a time stepping algorithm which preserves discrete mo- mentum and continuous energy 4) extensions of 1-3 above to the cases of plastic and in- elastic impacts, as well as automated distinction between plastic and non-plastic impacts when using the time- stepping method 5) application to three example systems: Newton’s cradle, billiards, and a two dimensional legged robot model robot which includes both external forcing and Coulomb friction The rest of the paper is structured around the above con- tributions. Section II goes over the collision model in the simple case of one isolated impact and introduces notation that will be used throughout the rest of the paper. Section III expands the basic collision model by extending it to the case of simultaneous collisions and explains its use in the example case of Newton’s cradle. In Sec. IV we examine existence arguments for the solutions of our algorithm, giving an upper bound on the number of iterations needed to achieve a feasible result for the case of two simultaneous contact points. The existence results shown here correspond to and extend those presented in [30], substituting a rigorous mathematical proof in place of a geometric argument. In Sec. V we investigate the uniqueness of the solutions obtained by the propagative method, for both the case of   n   perfectly elastic and perfectly plastic simultaneous impacts. For the plastic case the solution is trivially unique. However, for the elastic case we find that, while in certain cases, including Newton’s cradle, our approach gives unique solutions, in general we have a non-unique, al- though finite, number of outcomes. Similarly, standard implicit methods do not provide uniqueness results for their solutions, as the final solution is usually dependent on both the initial 2
condition and the particular solver. Our approach has the advantage of offering a set of countably many solutions, all of which satisfy the LCP criteria. Moreover, the non-unique solutions that we present have been posited to correspond to uncertainty in the elastic body physical system [2]. This gives rise to an interesting question: could we leverage control of the configuration at impact to always obtain predictable collision outcomes when dealing with simultaneous impacts? In Sec. VII-A we give an affirmative answer to this question— an extension to our work in [31]—and present its use on two example systems: an extension of the simplified system in Fig. 2a and a two-legged and tailed robot intended for locomotion seen in Fig. 1. All the simulation results are based on a time stepping scheme consisting of variational integrators which is briefly described in Sec. VI and the foundations of which were discussed in more detail in [32] and [30]. II. B ASIC   I MPACT   M ODEL  The following presents well known results concerning the equations governing impacts treated with an impulsive ap- proach. We present this both for reference and to introduce the reader to notation. Of particular note is our definition of inner product and norm on the cotangent space—see (5)— which makes use of the dual of the kinetic energy metric on the tangent space. Assume a simple mechanical system described by config- uration   q   and Lagrangian   L ( q,    ̇ q ) =   K ( q,    ̇ q )   −   V   ( q,    ̇ q ) . The equations of motion for such a system can be derived by applying a variational principle, and are the known Euler- Lagrange equations [33]:  d dt ∂L ∂    ̇ q   ( q,    ̇ q )   −   ∂L ∂q   ( q,    ̇ q ) = 0 .   (1) For a collision at time   t ∗   we say that the velocity at that time is   infeasible   if it points “into” the contact manifold—opposite to the contact manifold normal: D φ ( q ∗ )   ·    ̇ q ∗   <   0 ,   (2) where   q ∗   =   q ( t ∗ )   ∈   Q   is the configuration at time of impact and   φ   :   Q   →   R   is the gap function describing the contact manifold   C   ⊂   Q   [16]—it takes positive values in the feasible region of space, negative values in the infeasible region, and is identically zero when   q   is on the contact manifold. D φ   is the first derivative of the gap function with respect to its argument, and hence is a covector that belongs to the cotangent space at   q ∗ : D φ ( q ∗ ) =   ∂φ ∂q   ( q ∗ )   ∈   T   ∗  q ∗   Q , and provides the normal to the contact manifold. Using a variational approach one finds that the equations governing a single elastic collision at time  t ∗   are:  ∂L ∂    ̇ q  ∣ ∣ ∣ t +  ∗  t − ∗  =   λ D φ ( q ∗ ) ,   (3a)  [   ∂L ∂    ̇ q   ·    ̇ q   −   L  ] t +  ∗  t − ∗  = 0 .   (3b) Here   λ   is   a   Lagrange   multiplier   such   that   λ D φ   can   be interpreted as the impulse imparted to the system through impact. The two equations have a classical interpretation: (3a) is the conservation of momentum—tangentially to the impact manifold—and (3b) is the conservation of energy through the impact. From here on, we work under the assumption that we are dealing with a non-degenerate simple mechanical system such that potentials do not depend on the velocity    ̇ q , and the kinetic energy term is quadratic in    ̇ q   . Under this assumption, we can write the Lagrangian as   L ( q,    ̇ q ) =   1 2    ̇ q T   M   ( q )  ̇ q   −   V   ( q ) , where  M   ( q )   is the mass matrix   M   ( q ) =   ∂    ̇ q    ̇ q   L ( q,    ̇ q )   and is positive definite. The system in (3) becomes   ̇ q T   ( t +  ∗   ) M   −    ̇ q T   ( t − ∗   ) M   =   λ ∗ D φ ∗ ,   ̇ q T   ( t +  ∗   ) M    ̇ q ( t +  ∗   ) =  ̇ q T   ( t − ∗   ) M   ( t − ∗   ) ,  where, for ease of notation, we dropped the   q -dependency of  M   and   φ ∗ , implicitly assuming they are evaluated at   q ∗ , the impact configuration. We can rewrite these equations as  p +   =   p −   +   λ u ,   (4a)  p + M   − 1 p +T   =   p − M   − 1 p − T ,   (4b) where   p ±   =    ̇ q T  ∗ ±   M   is the momentum before and after the collision and   u   =   D φ T  ∗   is the normal to the manifold of impact (for the purposes of this paper we use bold notation to denote covectors, which in our case are elements of the cotangent bundle—e.g.   p ,   u   ∈   T   ∗  q ∗   Q —while regular script denotes vectors:   q   ∈   Q ,    ̇ q   ∈   T q ∗   Q , etc.) Also, for the remainder of the paper, the norms and dot products between covectors are assumed to be those defined under the local kinetic energy metric [34]:  〈 u ,   v 〉 g   =   u M   − 1 v T   ,   (5a)  ‖ u ‖ 2  g   =   〈 u ,   u 〉   .   (5b) Using   this   notation,   we   redefine   infeasibility   in   terms   of momentum. We say that   p −   is   infeasible   with respect to the contact manifold represented by the normal   u   if  〈 p − ,   u 〉  g   <   0 .   (6) By the same token,   p −   is   feasible   if   〈 p − ,   u 〉 g   ≥   0 . Assuming an infeasible   p − , we solve for   p +   using (4):  p +   =   p − Γ( u ) ,   (7a)  Γ( u ) =   I   −   2 M   − 1   u T   u  ‖ u ‖ 2  g  =   I   −   2   M   − 1 u T   u u M   − 1 u T   .   (7b) Here   Γ( u )   is a momentum (reset) map [2] that describes an instantaneous change in momentum due to an impact with a manifold normal   u . Equation (7a) is only one of the two solutions to the system in (4), but is the only feasible one. Indeed, after eliminating   p +   we are left with a quadratic equation in   λ :  λ 2   ‖ u ‖ 2  g   + 2 λ   〈 p − ,   u 〉 g   = 0 . The feasible value is given by  λ   =   − 2   〈 p − ,   u 〉 g   /   ‖ u ‖ 2  g   . Note that if we write (2) using the notation just introduced, we have that   λ >   0 . This is consistent with the LCP formulation of contact, as   λ >   0   satisfies the classic complementarity conditions. The mapping   Γ( u ) , as defined in (7a), has several properties that will be useful in Sec. III, Sec. IV, and Sec. V. First,  Γ( u ) 2   = Γ( u )Γ( u ) =   I,   (8) 3
which means that resolving an impact across a manifold twice in a row returns the original momentum. Also,   Γ   is not dependent on the magnitude of   u , only on its direction:  Γ( α u ) = Γ( u ) .   (9) III. S IMULTANEOUS   C OLLISIONS  Simultaneous impact can occur if two contacts are made at the same time or if one contact is already present when a second impact occurs. When an impact is assumed to be plastic the equations governing the interaction are  ∂L ∂    ̇ q  ∣ ∣ ∣ t +  ∗  t − ∗  =   ∑  i ∈   U ∪   V  λ i D φ i ( q ∗ ) ,   (10a) D φ i ( q ∗ )  ̇ q +  ∗   = 0 ,   λ i   ≥   0 ,   ∀ i   ∈   U   ∪   V ,   (10b) where   U   is the set of indices of manifolds already in contact while   V   is the set of indices of new contact manifolds. These equations generate a unique solution, which corresponds to eliminating the portion of    ̇ q ∗   that is orthogonal—under the kinetic metric—to the manifolds of collision   φ   at the time of impact. Using the notation introduced in Sec. II, the outcome of a plastic impact is given by  p +   =   P null(Span { D φ i ( q ∗   ) | i ∈   U ∪   V } )  ( p − )   ,   (11) where   P S   ( p )   represents the projection of covector   p   onto a subspace   S . In essence,   p +   is always tangent to all contact manifolds—or, equivalently, it is orthogonal to all contact manifold normals at the current configuration. The uniqueness of the result is lost when considering elastic impacts:  ∂L ∂    ̇ q  ∣ ∣ ∣ t +  ∗  t − ∗  =   ∑  i ∈   U ∪   V  λ i D φ i ( q ∗ ) ,   (12a)  [   ∂L ∂    ̇ q   ·    ̇ q   −   L  ] t +  ∗  t − ∗  = 0 ,   (12b)  λ i   ≥   0 ,   ∀ i   ∈   U   ∪   V ,   (12c) These equations form a fully determined system only when there is a single term in the summation on the right hand side of (12a)—in which case they reduce to the case treated in Sec. II. Otherwise, the equations governing the impact dynamics become underdetermined (there are more variables than equations). This gives rise to a whole continuum of solutions. For generating trajectories in simulation, an element of this continuum must be chosen. An a priori relation between lambdas could be chosen in order to solve this dilemma. However, it is unclear what physical principle to use. Instead, we expand on a version of the propagative method discussed in   [2].   The   reason   behind   using   a   propagative   model   of simultaneous impact is twofold: the method gives unique and correct results in simple, intuitive cases in which other methods fail—e.g. Newton’s cradle—and it also provides at most a finite number of valid solutions in other cases, as discussed in Sec. V. We investigate the simplest case of simultaneous impacts: that of two manifolds of contact, such that   U   ∪   V   =   { a, b } . Suppose that the impact occurs across two manifolds at the exact same time, such that (12) becomes  ∂L ∂    ̇ q  ∣ ∣ ∣ t +  ∗  t − ∗  =   λ a D φ a ( q ∗ ) +   λ b D φ b ( q ∗ ) ,   (13a)  [   ∂L ∂    ̇ q   ·    ̇ q   −   L  ] t +  ∗  t − ∗  = 0 ,   (13b)  λ a   ≥   0 ,   λ b   ≥   0 ,   (13c) which has the same number of equations as (3), but one extra variable. Instead, a propagative approach consists of applying (7a) repeatedly until a feasible momentum is found—whether this is even operationally valid is something we address in Sec. IV. Using the notation of the previous section, where   p  represents momentum and   u a   and   u b   represent the contact manifold normals, we would have a sequence of operations such as  p 1   =   p − Γ( u a )   with   〈 p 1 ,   u b 〉 g   <   0 ,  p 2   =   p 1 Γ( u b )   with   〈 p 2 ,   u a 〉 g   <   0 , . . .   =   . . .   with   . . .  p +   =   p n − 1 Γ( u b )   with   〈 p + ,   u a  〉  g   ≥   0 ,   〈 p + ,   u b  〉  g   ≥   0 .  The above sequence generates a solution of the form  p +   =   p −  n ∏  i =1  Γ( w i ) ,   (14) where   ∏ n i =1   x i   is short for “the product   x A   ·   x B   · · ·   x n ” and  { w i }   is a sequence that alternates between   u a   and   u b   as in  { w i }   =   { u a ,   u b ,   u a ,   u b , . . .   } . The solution offered by (14) is certainly not unique since either the sequence or the number of terms might vary. For example, a solution where   { w i }   =  { u b ,   u a ,   u b ,   u a , . . .   }   might be equally valid, provided that   p +  is feasible. Similarly, since   Γ 2 ( u ) =   I , prepending the first element an even number of times is equivalent to applying the original sequence. Thus, the previous alternating sequence is equivalent to   { w i }   =   { u b ,   u b ,   u b ,   u a ,   u b ,   u a , . . .   } . Note that we are working with infinite sequences, as we do not want to assume a priori that a feasible momentum can be found in a finite number of steps. The question of existence (and thus finiteness of the length of the sequence) is, contrary to intuition, non-trivial. For a discussion and proof of termination for two surfaces, see sec. IV. In general, for a set of more than two contact manifold normals  S   =  {  u   =   Dφ i ( q ∗ )  ‖ Dφ i ( q ∗ ) ‖ g  ∣ ∣ ∣ ∣   φ i ( q ∗ ) = 0 ,   〈 p − ,   u 〉  g   <   0  }  we will have several choices of   w i   ∈   S   and a number of terms such that   〈 p + ,   u 〉 g   ≥   0   for all   u   ∈   S.   With the above in mind we propose a minimality condition which will reduce the possible mapping sequences and, in certain cases such as discussed in Sec. V, will provide uniqueness.  Definition 1 (Minimality Condition):   We say that a sequence  W   =   { w i }   of contact manifold normals is   minimal   with respect to an infeasible momentum   p   and Lagrangian   L   if the following hold 4
1) the application of the corresponding sequence of reset maps   { Γ( w i ) }   generates a feasible momentum:  〈  p  n ∏  i =1  Γ( w i ) ,   w  〉  g  ≥   0 ,   ∀ w   ∈   S  2) no proper subsequence of   W   can generate a feasible momentum. We proceed to only consider minimal sequences of reset maps, with two direct consequences which help us make stronger statements in regards to existence and uniqueness of feasible solutions to (13). The first, and most important, is that we do not consider solutions obtained when continuing to apply reset maps to an already feasible solution. The second consequence is that none of the mapping sequences we consider will have consequent members that are identical, which is guaranteed by the following lemma:  Lemma 1:   Given a sequence   { w i }   which is minimal with respect to some momentum   p − , we must have that   w j   6   =  w j +1   for all   j . 2  An important consequence of Lemma 1 is that, in the case of two contact manifolds, described by   u   and   v , any minimal sequence of normals has the form of an alternating sequence of   u   and   v . This fact will be central to later results.  Example 1 (Newton’s cradle):   Consider the system shown in Fig. 2a. We can define manifolds of contact by  φ 1 ( q ) =   x B   −   x A   −   2 r, φ 2 ( q ) =   x C   −   x B   −   2 r,  u   =   D φ 1 ( q ) = [ − 1 ,   1 ,   0]   ,  v   =   D φ 2 ( q ) = [0 ,   − 1 ,   1]   .  The configuration space is three dimensional and the contact manifolds are two planes. Under the assumption that all masses are equal   M   =   mI , we have that  Γ( u ) =     0   1   0 1   0   0 0   0   1      ,   Γ( v ) =     1   0   0 0   0   1 0   1   0      .  It takes three iterations to find a feasible solution, making  { u ,   v ,   u }   and   { v ,   u ,   v }   the only minimal sequences with respect to an infeasible momentum   p − :  p +   =   p − Γ( u )Γ( v )Γ( u ) =   p − Γ( v )Γ( u )Γ( v ) =   p −      0   0   1 0   1   0 1   0   0      .   (15) Note that we obtained the same final result regardless of the order of impacts, which is not expected to hold for general systems. For example, if we were to choose unequal masses for the balls, the solution would not be, in general, unique. This indicates that the uniqueness of the solution is related to an interplay between both the geometry of the system and its inertia tensor at the time of impact. We will present and discuss sufficient conditions for uniqueness in Sec. V and use the results of that section to develop the main contribution of the paper, the impact design approach in Sec. VII-A.  2 For clarity, the proof of this and all following lemmas can be found in the appendix.  A. Extension to Inelastic Collisions  In general, collisions models include a coefficient of resti- tution, to account for energy loss during impact. While the coefficient   of   restitution   is   usually   defined   as   a   ratio   of collinear impulses [35], it cannot be used with a propagative model since propagative models deal with multiple impulse exchanges. However, the restitution coefficient can be defined using the energy loss through the impact:  R   =  √  1   −   ∆ E E p  ,   (16) where   E p   is the energy that would have been lost during the collision through a perfectly plastic impact and   ∆ E   is the energy lost when considering the coefficient of restitution   R . A coefficient of restitution of zero will determine a plastic impact while a value of one will generate a perfectly elastic impact. For all other inelastic impacts we represent the solution as a convex combination between the plastic and elastic outcomes, parametrized by   α :  p +   =   α   p +  e   + (1   −   α ) p +  p   ,   α   ∈   [0 ,   1] ,   (17) where   p +  e   is the momentum outcome for a perfectly elastic collision and   p +  p   is the momentum outcome of a perfectly plastic collision. Using (17) in the right hand side of (16), along with the observation that   E p   =   1 2  ∥ ∥ p +  p  ∥ ∥ 2  g   , we obtain:  ∥ ∥ p + ∥ ∥ 2  g   −   ∥ ∥ p +  e  ∥ ∥ 2  g   =   R 2   (∥ ∥ p +  e  ∥ ∥ 2  g   −   ∥ ∥ p +  p  ∥ ∥ 2  g  )  .  The   values   of   p +  e   and   p +  p   are   given   by   (14)   and   (11), respectively. Since   p +  p   is an orthogonal projection of   p − , and  p +  e   is conserved in the direction of   p +  p   , we have that  ∥ ∥ p +  p  ∥ ∥ 2  g   =   〈 p +  p   ,   p +  e  〉  g   .   (18) We solve for   α   and use (18) to obtain  α   =   √ 1   −   R 2 ,   R   ∈   [0 ,   1] .  IV. E XISTENCE  While the previous section presents an overview of the propagative method for solving simultaneous impacts, it also raises two important questions regarding the same method: do solutions always exist and, if they do, are they unique? The current section addresses a special case of existence by proving that, for a simultaneous impact involving two contact manifolds a minimal sequence of mappings exists and it is finite. The uniqueness of the corresponding momentum outcome is investigated in Sec. V. In what follows we assume that we are dealing with a simultaneous elastic impact involving two contact manifolds. A representation of the contact manifolds at the impact con- figuration is given by their normals   u   and   v , which, without any loss of generality, can be assumed as being of unit length and not collinear:  ‖ u ‖ g   =   ‖ v ‖ g   = 1 ,   (19a)  u   6   =   ± v .   (19b) 5
Let   T   be the intersection of the two hyperplanes orthogonal to   u   and   v , defined by  T   = null   { u ,   v }   =  {  p   ∈   T   ∗  q ∗   Q | 〈 p ,   u 〉 g   =   〈 p ,   v 〉 g   = 0  }  .  Additionally, let   N   be the plane defined by the two covectors  N   = span { u ,   v }   =  {  p   ∈   T   ∗  q ∗   Q | 〈 p ,   w 〉 g   = 0 ,   ∀ w   ∈   T  }  .  Notice   that,   by   definition,   T   and   N   are   orthogonal   and complementary. We will make use of these properties in the following lemmas. We start out by presenting two lemmas which allow us to equate—through the use of an orthogonal projection—solutions   p   ∈   T   ∗  q ∗   Q   to the lower dimensional  p N   ∈   N . The lemmas provide the link between solving the problem of existence in the lower dimensional space   N  and the problem of existence in   T   ∗  q ∗   Q : the former becomes sufficient in order to prove the latter. Thus, the substantial part of the proof for Theorem 1 consists of showing existence when   dim( N ) = 2 . We do so by using a geometric argument involving the angle between successive reflection mappings of a covector   r   = ( u   +   v ) /   ‖ u   +   v ‖ g   and a momentum   p 0 , which is related to the pre-impact momentum   p −   through an orthogonal projection. Figure 3 illustrates the geometric argument used in Theorem 1. For the rest of the this section we will use the following notation to denote the orthogonal projection of a covector onto a set   S :  P S ( p ) = arg min  p ∗   ∈ S  ‖ p   −   p ∗ ‖ g   .  This notation, along with standard properties of orthogonal projections, will be heavily used in the following lemmas. First we show that the feasibility of   p   is equivalent to the feasibility of its component that can be represented as a linear combination of   u   and   v .  Lemma 2:   Let   p ,   u ,   v   ∈   T   ∗  q ∗   Q   with   u ,   v   satisfying the restrictions in (19). Then,   p   is feasible   iff   p N   =   P N ( p )   is feasible. Next, we show that any sequence of reflection transformations will only affect the component of momentum that is in the span of   u   and   v , leaving the component orthogonal to that plane unchanged.  Lemma 3:   Given   p   ∈   T   ∗  q ∗   Q   and a sequence of trans- formations   Γ( w i )   with   w i   ∈ { u ,   v }   that map   p   to   p f   =  p   ∏  i   Γ( w i ) , we always have that  p f   =   P T ( p ) +   P N ( p )   ∏  i  Γ( w i ) ,  Next, we present a lemma that equates feasibility in   N   to a trigonometric condition on the inner product between   p   and a covector   r :   p   is feasible   iff   the angle between   p   and   r   is smaller than half the angle between   u   and   v .  Lemma 4:   Let   p ,   u ,   v ,   r   ∈   N   such that   r   =   u + v  ‖ u + v ‖ g  and  u   6   =   ± v . Let   γ   = arcsin   〈 r ,   u 〉 g   . Given the above we have that   p   is feasible   iff  〈 p ,   r 〉 g   ≥ ‖ p ‖ g   cos   γ.   (20) Note that for the backwards implication we did not use the fact that   p   ∈   N , which means that this implication can be easily generalized to more than two contact manifolds. The same cannot be said about the forward implication, as counterexamples can be easily found when   p   / ∈   N : consider the case where   〈 p ,   u 〉 g   =   〈 p ,   v 〉 g   = 0   but   ‖ p ‖ g   6   = 0 . We use the result of Lemma 2 to simplify the problem of existence of solutions in the higher dimensional cotangent space to an equivalent problem in a two-dimensional subspace. In particular, since reflection transformations can only affect the   P N ( p )   component of the momentum, it is enough to show that, given a finite number of transformations   Γ( w i ) , we can transform   P N ( p )   into a valid momentum   p N   ∈   N  such that   〈 p N ,   u 〉 g   ≥   0   ≤   〈 p N ,   v 〉 g   . We can   then   find the corresponding transformation of the original momentum through the same sequence of mappings, which, according to Lemma 3 is guaranteed to be feasible as well. Lemma 4 gives us a way to rewrite the feasibility conditions into a form dependent on the inner product between the momentum and a unit covector   r   in the cotangent space. Before we show the main result, we present one last lemma which gives shows that the mapping   Γ   is conformal under the kinetic metric:  Lemma 5:   Given a momentum map   Γ( w )   and two momenta  p a ,   p b   ∈   T   ∗  q ∗   Q   we have that  〈 p a Γ( w ) ,   p b Γ( w ) 〉 g   =   〈 p a ,   p b 〉 g   .   (21) In   words,   Γ( w )   is   guaranteed   to   be   a   conformal   (angle preserving) mapping under the kinetic metric. This enables us to formulate and prove the main result:  Theorem 1:   Given a momentum   p 0   ∈   T   ∗  q ∗   Q   and two man- ifolds of impact with linearly independent normals   u   6   =   ± v  at the point of simultaneous contact, there exists a minimal sequence   { w i } , as per Def. 1, of length   n   ≤ d π/γ e   such that  p f   =  n ∏  i =1  Γ( w i )  is feasible, where   γ   is the angle defined in Lemma 4.  Proof:   We start by considering a momentum   p N  0   ∈   N . The subspaces   T   and   N   are both orthogonal and complementary, which makes it possible for us to write  p 0   =   P T ( p 0 ) +   P N ( p 0 ) .  We focus our attention on   p N  0   =   P N ( p 0 ) , since Lemma 2 guarantees that any feasibility results on   p N  f   extend to   p f   . This allows us to work in a two dimensional vector space isomorphic to   R 2 . Furthermore, Lemma 3 gives us that we can find   p f   by applying the same sequence of mappings to  p 0   as that used when obtaining   p N  f   from   p N  0   . In what follows we make use of two properties of the momentum map defined in (7). The first is invertibility. In fact, a mapping   Γ( w )   is its own inverse, since   Γ( w )Γ( w ) =   I . In addition, Lemma 5 shows that   Γ( w )   is a conformal momentum map under the kinetic energy metric: it preserves angles— inner products—according to this metric. Now let   r 0   =   u + v  ‖ u + v ‖ g  . Using the fact that   Γ( w i )   are invertible and taking into consideration the result of Lemma 4, it is sufficient to find a sequence of mappings   Γ( w i )   that maps  r   to  r f   =   r 0  n ∏  i =0  Γ( w i ) ,  6
such that  〈 r f   ,   p N  0  〉  g   ≥  √  1   − 〈 u ,   v 〉 g  2  ∥ ∥ p N  0  ∥ ∥ g   .  Applying the sequence of mappings to   p 0   in reverse order will generate a feasible  p f   =   p 0  n ∏  i =0  Γ( w n − i ) .  Given the above, the following has to hold:  Lemma 6:   For any minimal sequence   { w i }   we have that  〈 r i ,   r 0 〉 g   = cos(2 iγ ) ,  where, as in Lemma 4,   γ   = arcsin   〈 r 0 ,   u 〉 g   = arcsin   〈 r 0 ,   v 〉 g  and   r i   =   r 0  ∏ i j =0   Γ( w j   ) . Thus, for any value of   i   we can calculate   〈 r i ,   r 0 〉 g   = cos(2 iγ ) . Knowing this angle, the fact that   r i   ∈   N , and  ‖ r i ‖ g   = 1 , we have only two possible   r i   to satisfy these conditions. One of the two covectors is obtained by setting  w 0   =   u   while the other is obtained by setting   w 0   =   v —see Fig. 3. If we account for both these possibilities and consider the set  R N   =  {  r n  ∣ ∣ ∣ ∣ n   ≤  ⌈   π γ  ⌉ }  ,   (22) we have that, for any   r i   ∈   R N   , the smallest angle between it and any other   r j   ∈   R N   is at most   2 γ . Formally,  min  i   max  j   〈 r i ,   r j   〉 g   ≥   cos(2 γ ) .  The implication of this is that for any vector   p N ,   ∃ i   ≤   N   s.t.  〈 p ,   r i 〉 g   ≥   cos( γ ) . This proves our theorem, since applying the mappings used to get   r i   in reverse order to   p 0   guarantees that  〈  p N  f   ,   r 0  〉  g  ≥   cos( γ ) . The above proof shows that, in the case of two contact manifolds with distinct normals at the point of contact, a feasible solution can be obtained through a finite number of successive applications of the reflection mappings defined by  Γ( u )   and   Γ( v ) . Generalizing to more simultaneous impacts is complicated by the lack of an alternating structure in   w i   in such cases. Considering all this, we leave the proof of such existence an open question to further study. V. U NIQUENESS  In   this   section   we   present   several   results   that   help   us determine uniqueness of simultaneous impacts. In particular, Theorem 2 and Theorem 3 guarantee uniqueness of the impact results in the cases that the inner product between the two con- tact manifold normals takes a value of   0   or   ± 0 . 5 , respectively. Also, as a consequence, if the impact involves a higher number of contact manifolds which are all pairwise orthogonal at the impact configuration, the result of applying the reset maps will be unique, regardless of order of application. We first show that orthogonality between the normals of the contact manifolds at the simultaneous impact configurations is both sufficient and necessary for a unique feasible solution to be found in exactly two steps.  Lemma 7 (Orthogonality Condition):   For two contact man- ifolds described by their normals   u   6   =   ± v , and any infeasible momentum   p   with   〈 u ,   p 〉 g   ≤ 〈 v ,   p 〉 g   <   0 , we have that  p f   =   p Γ ( u ) Γ ( v ) =   p Γ ( v ) Γ ( u )   is feasible   (23)  iff  〈 u ,   v 〉 g   = 0 .   (24) As part of the above result we have that two reset maps commute if and only if the covectors that generate them are orthogonal. This helps us show the following important theorem.  Theorem 2:   Given a set of   n   contact manifolds described by their normals   u i   such that they are all pairwise orthogonal  〈 u i ,   u j   〉 g   = 0 ,   ∀ i, j,   i   6   =   j,  and any infeasible momentum   p   with   〈 u i ,   p 〉 g   <   0 ,   ∀ i , we then have that  p f   =   p  n ∏  i  Γ( u i )   (25) is the same for any ordering of the indices   i   and the result is feasible.  Proof:   It follows directly from Lemma 2 that, if   u i  are pairwise orthogonal then   Γ( u i )   pairwise commute, which gives that the order of terms in the product of (25) is irrelevant. It remains to show that   p f   is feasible. We do this by showing that the application of a given mapping   Γ( u k )   affects the inner product between the mapped momentum and no other normal but   u k . Indeed, assume that we have a   p k   and a   u i   6   =   u k . Then we will have that  〈 p k Γ( u k ) ,   u i 〉 g   =  〈  p k   −   2   〈 p k ,   u k 〉 g   u k ,   u i  〉  g  =   〈 p k ,   u i 〉 g   −   2   〈 p k ,   u k 〉 g   〈 u k ,   u i 〉 g  =   〈 p k ,   u i 〉 g   .  Thus, applying all of the mappings   Γ( u i )   only once, in any order, we make sure the sign of the inner product   〈 p f   ,   u i 〉 g  is positive for all normals   u i . Hence,   p f   is feasible. The orthogonality condition in Lemma 2 was found under the assumption that the system undergoes two impacts before a feasible exit velocity is found. Sometimes this might not be the case—e.g. due to design constraints—and systems where more than two impacts are needed to find a feasible exit velocity need to be considered. Such a system is Newton’s cradle, where, for equal masses, the inner product is   〈 u ,   v 〉 g   =   − 0 . 5 . The following theorem provides sufficient conditions for a three-stage impact to generate a unique outcome:  Theorem 3 (Three Stage Impacts):   Given two contact man- ifolds described by their normals   u   6   =   ± v   and   p   that satisfies  〈 u ,   p 〉 g   ≤ 〈 v ,   p 〉 g   <   0 , we have that  p f   =   p Γ ( u ) Γ ( v ) Γ( u ) =   p Γ ( v ) Γ ( u ) Γ( v )   (26) is feasible   iff  〈 u ,   v 〉 g   =   −   1 2   .   (27) 7
v  u  r 0  N   '   R 2  p 0  γ   γ  r 0  N  r 1  r 1  p 0  2 γ  2 γ   λ 1 v λ 1 u   r 0  N  r 1  r 1  r 2  r 2  ... r n p 0  p f  γ γ γ γ  r 0  N  r 1  r 1  r 2  r 2  p 0  2 γ  2 γ λ 2 v λ 2 u  (a)   (d) (c) (b) Fig. 3: A sketch of the subspace   N   when the kinetic energy metric is identical to the Euclidean metric. The figure gives an interpretation to equations (40) and (22). The dotted circle is the intersection of the unit sphere in   T   ∗  q ∗   Q   with   N   and the dotted lines represent the contact manifolds tangents—they are orthogonal to   u   and   v   and represent the axis across which the reflection transformations operate. Figure (a) shows the two normals   u   and   v , the additional variables   r 0  and   γ   as well as an initial infeasible momentum   p 0 ; (b) shows the two potential applications of a momentum map to   r 0   and the two corresponding values for  r 1 ; a second sequence of mappings is presented in (c)—notice the   2 γ   angle increase between   r 1   and   r 2 ; finally, in (d), we have gone through   n   mappings and obtained an   r n   that is within at most   γ   of   p 0 . Applying the mappings that generated   r n   to   p 0   in reverse order we obtain   p f   which is guaranteed to be within   γ   of   r 0 , and hence feasible.  Proof:   Through the use of (7) and some algebra, the condition in (26) can be shown equivalent to  (  1   −   4   〈 u ,   v 〉 2  g  ) (  〈 p ,   u 〉 g   u   − 〈 p ,   v 〉 g   v  )  = 0 .  Since we assumed that   u   6   =   ± v , it must be that  〈 u ,   v 〉 g   =   ±   1 2   .  In the case that   〈 u ,   v 〉 g   =   1 2   , we substitute back into (27) to obtain:  p f   =   p   −   2( v   −   u )  (  〈 p ,   v 〉 g   − 〈 p ,   u 〉 g  )  ,  〈 p f   ,   u 〉 g   = 2   〈 p ,   u 〉 g   − 〈 p ,   v 〉 g   ,  〈 p f   ,   v 〉 g   = 2   〈 p ,   v 〉 g   − 〈 p ,   u 〉 g   ,  which implies that at least some of the infeasible   p   are mapped to an infeasible   p f   . So   〈 u ,   v 〉 g   =   1 2   does not work, thus showing the forward implication. On the other hand, if (27) holds we obtain  p f   =   p Γ ( u ) Γ ( v ) Γ( u ) =   p   −   2  (  v   −   2   〈 u ,   v 〉 g   u  ) (  〈 p ,   v 〉 g   −   2   〈 p ,   u 〉 g   〈 u ,   v 〉 g  )  =   p   −   2 ( v   +   u )  (  〈 p ,   v 〉 g   +   〈 p ,   u 〉 g  )  =   p Γ ( v ) Γ ( u ) Γ( v )  Finally, to show that   p f   is feasible, we calculate  〈 p f   ,   u 〉 g   =   〈 p ,   u 〉 g   −   2  (  1   − 〈 u ,   v 〉 g  ) (  〈 p ,   u 〉 g   +   〈 p ,   v 〉 g  )  =   − 〈 p ,   v 〉 g   >   0 ,  〈 p f   ,   v 〉 g   =   〈 p ,   v 〉 g   −   2  (  1   − 〈 u ,   v 〉 g  ) (  〈 p ,   u 〉 g   +   〈 p ,   v 〉 g  )  =   − 〈 p ,   u 〉 g   >   0 .  To illustrate this, consider the Newton’s cradle system dis- cussed in Sec. I. Assuming all three spheres have the same mass   m , the normals to the contact manifolds are  u   = [ − 1 ,   1 ,   0] √ m/ 2 ,  v   = [0 ,   − 1 ,   1] √ m/ 2 .  The dot product between these two normalized covectors is  〈 u ,   v 〉 g   =   m − 1 u   I   v T   =   −   1 2   ,  which is consistent with the result of Theorem 3. This prop- erty of Newton’s cradle implies outcome uniqueness for any combination of initial momenta of the spheres at the moment of simultaneous impact. VI. T IME - STEPPING   M ETHOD FOR   S IMULTANEOUS  I MPACT  In the previous sections we have only discussed impact resolution, assuming that we are already given an impact momentum and impact manifolds. However, we are interested in simulating a time interval that contains an impact. Further- more, our methods for impact resolution presented in (3) and (14) are formulated in terms of the momentum before and after the impact, the result and proof of Theorem 1 is also specified in terms of momentum. We want, therefore, an integration scheme that works directly with momentum. The method of variational integrators [36]–[39] is ideal in this situation, since, by design, it will conserve a quantity known as the   discrete momentum —obtained through the discrete Legendre transform from a pair of configurations. The discrete momentum is an approximation of the continuous momentum at a given time, the error between the two vanishing in the limit of the time step going to zero. There are other advantages to using variational integrators: they are known to preserve the symplectic form [40] and conserve the average energy of the system over a large number of time steps [37]. Finally, recent work has been taking advantage of variational integrator methods in order to generate optimal controllers for complex systems [41]–[44]. Variational integrators are obtained by discretizing the ac- tion sum directly:  L d ( q i , t i , q i +1 , t i +1 )   ≈  ∫   t i +1  t i  L ( q,    ̇ q, t )   dt  = ( t i +1   −   t i )   L  (   q i   +   q i +1  2   , q i +1   −   q i  t i +1   −   t i  , t i   +   t i +1  2  )  .  8
Note that the discrete Lagrangian depends only on configura- tion variables, and not on velocity information. There are also several quadrature rules one can apply for the discretization. Here we have used the midpoint rule. The discrete equivalent of the Euler-Lagrange equations is the set of Discrete Euler- Lagrange equations [44]:  ∂ ∂q k  L d ( q k − 1 , t k − 1 , q k , t k ) +   ∂ ∂q k  L d ( q k , t k , q k +1 , t k +1 ) = 0 ,  (28) that can be thought of as a mapping from two known config- urations   q k − 1   at time   t k − 1   and   q k   at time   t k to an unknown configuration   q k +1   at time   t k +1 . A common interpretation of (28) is that they enforce the conservation of   discrete momentum   [37], which is defined through the use of the discrete momentum maps  F − ( t a , t b ) =   ∂ ∂q a  L d ( q a , t a , q b , t b ) ,  F + ( t a , t b ) =   ∂ ∂q b  L d ( q a , t a , q b , t b ) .  Using this notation, (28) becomes  F + ( t k − 1 , t k ) +   F − ( t k , t k +1 ) = 0 ,   (29) which states that the forward momentum   F +   at the end of the  ( t k − 1 , t k )   interval has to equal the backward momentum   F −  at the beginning of the following interval,   ( t k , t k +1 ) . Equation (29) is a discrete equivalent of the Euler-Lagrange equations and is known as the DEL (discrete Euler-Lagrange) set of equations. In case of an impact at time   t ∗ , we apply the update map from (14) to the discrete momentum, and we solve  F + ( t k , t ∗ )  N ∏  i  Γ( w i ) +   F − ( t ∗ , t k   + 1) = 0 .   (30) The equations in (29) and (30) are, for all but the simplest systems, nontrivial. We use [36] in which the terms in these equations are derived using a tree structure and used in a root finding algorithm. This is the method we have used when solving the dynamics away from impact for the running mechanism described in the following section. VII. I MPACT   D ESIGN  In this section we apply the uniqueness results of Sec. V to example systems. We calculate combinations of geometries and configurations that have unique outcomes for two mechan- ical systems: a billiard ball break and a tailed, running biped. These calculations are done analytically for the billiards and numerically for the biped.  A. Billiards  Consider the general billiard break shown in Fig. 4. Billiard  c   acts as the cue ball in this situation, imparting momentum to the other two billiards through a simultaneous collision. The configuration vector for this system is  q   = [ x a , y a , x b , y b , x c , y c ] T   ,  and the two gap functions are  φ a ( q ) =   √ ( x a   −   x c ) 2   + ( y a   −   y c ) 2   −   ( r a   +   r c ) , φ b ( q ) =   √ ( x b   −   x c ) 2   + ( y b   −   y c ) 2   −   ( r b   +   r c ) .  The two normals at a point where both these functions are zero are D φ a ( q ∗ ) = [ x c   −   x a , y c   −   y a ,   0 ,   0 , x a   −   x c , y a   −   y c ]  r a   +   r c  ,  D φ b ( q ∗ ) = [0 ,   0 , x c   −   x b , y c   −   y b , x b   −   x c , y b   −   y c ]  r b   +   r c  .  The mass matrix of this system is diagonal and since we are assuming the impacts to be frictionless and are not expecting any energy in the rotational modes of the objects, we ignored the moments of inertia when calculating   M   . The dot product between the manifolds is  〈 D φ a ,   D φ b 〉 g   = ( x a   −   x c ) ( x b   −   x c ) + ( y a   −   y c ) ( y b   −   y c )  m c   ( r a   +   r c ) ( r b   +   r c )   .  Since   φ a   =   φ b   = 0 , we can rewrite this expression in terms of   θ   (see Fig. 4) using the law of cosines:  〈 D φ a ,   D φ b 〉 g  = ( x a   −   x b ) 2   + ( y a   −   y b ) 2   −   ( r a   +   r c ) 2   −   ( r b   +   r c ) 2  2 m c ( r a   +   r c )( r b   +   r c ) = cos( θ )  m c  .   (31) Thus, requiring that the impact manifolds be orthogonal under the kinetic energy metric is, in this case, equivalent to setting  θ   =   π/ 2 . It is interesting to note that the result in (31) depends neither on the masses of the billiards, nor on their radii. While keeping billiard   c   in contact with   a   and   b , we varied the angle  θ   continuously from the minimum value where   a   and   b   were also in contact—somewhere close to   π/ 6   in our case—up to m a  m b  m c  θ x y r a  r c  r b  ( x a , y a ) ( x c , y c ) ( x b , y b )  φ a  φ b  ~ v c  Fig. 4: The schematic of a planar mechanical system consisting of three billiard balls about to experience simultaneous impact. Ball   a   and   b   are stationary while ball   c   has an initial velocity   v c   such that the contact between  a   and   c   is simultaneous with the contact between   b   and   c . The masses of the billiards are proportional to the volumes of the respective spheres and all friction effects are ignored.  9
π  4  π  2 3 π  4   π  Angle between billiards   θ   (rad) 0 2 4 6 8 10 12 Outcome indeterminacy   ξ   ( % )  ξ   =   ‖ p ab   − p ba ‖ g  ‖ p 0 ‖ g Fig. 5: The indeterminacy of the outcome of the billiard break in Fig. 4 as a function of the angle between the billiards. As a measure of indeterminacy   ξ  we looked at the difference in momentum between the two possible outcomes under the kinetic energy metric and as a percentage of the total energy of the system.  π . The only billiard with initial velocity was   c   and   v c   was chosen to point along the bisector of angle   θ . As a measure of indeterminacy we looked at the difference in momentum between the two possible outcomes under the kinetic energy metric, relative to the total kinetic energy of the system:  ξ   =   ‖ p ab   −   p ba ‖ g  ‖ p 0 ‖ g  ,  where   p 0   is the momentum covector before impact,   p ab   is the momentum after impact when solving the   a   −   c   impact first and   p ba   is the other option, where we solve for the   b   −   c  impact before solving the   a   −   c   impact. Figure 5 shows the values taken by   ξ   as we varied   θ . An expected minimum exists at   θ   =   π , which corresponds to a grazing impact. Somewhat less intuitive is the minimum at   π/ 2 , which tells us that when the impact configuration is such that   θ   =   π/ 2   the simulation has no indeterminacy in its solution.  B. TREX  The second system we investigated is the tailed running mechanism (TREX 3 ) in Fig. 6. The model is inspired by the geometry of several legged locomotors, such as the RHex [7], IMPASS [5] and several others [3], [4], [6]. The system consists of a two part body, two articulated legs and a tail. The two knee joints and the joint between the two parts of the body are modeled as linear torsional springs and dash pots. The actuators are located at the hips. The tail is connected rigidly to the posterior body segment. The configuration variables for this system are  q   = [ x, y, θ, φ L , γ L , φ R , γ R , φ T   ] T  ,  which represent the position and the orientation of the anterior body, the left hip and knee angles, the right hip and knee angles  3   The name was chosen due to the geometric similarity between our mechanism and a commercially available dinosaur toy [45]. r 1  r 2  φ L  − γ L  − φ R  − γ R  − θ  − φ T  x y  Fig. 6: The schematic of the two-dimensional TREX model used in simulation. The articulations are at the hips, knees and tail. The degrees of freedom consist of the Cartesian coordinates of the main body   x   and   y   together with all the angles   φ ,   θ   and   γ , as shown. The floor is assumed to be horizontal and positioned at   y   = 0 .  and the tail angle. We fixed the densities of the body and limbs and the section area of the limbs to reasonable values— the body density is that of water,   1 g/cm 3 , the density of the limbs is that of carbon fiber,   2 g/cm 3 , and the cross section of the limbs was assumed to be   1 . 25 cm 2 —while leaving the leg segment lengths and radii of the two bodies as design variables. In order to enforce a tapering shape and reduce the number of variable parameters, we chose the length of the tail to be a function of the body radii  L tail   =   r 3 1   +   r 2 1   r 2   +   √ 2 √ r 4 1   r 2 2   +   r 3 1   r 3 2  r 2 1   −   r 2 2  −   ( r 1   +   r 2 ) ,  which imposes that the tail length be inversely proportional with the ratio of the two radii. The mechanism drive was generated through external forc- ing at the hips, implemented as described in detail in [36], [41], [42], [44]. The values of the forces were chosen by a standard PD controller with gains   K p   = 10 5   and   K v   = 10 4 . Coulomb friction was also added into the model for this system, in order to facilitate its moving forward. Friction was implemented by adding an external force in each independent tangent direction at the contact point. The coefficients of these forces are found using the   maximum dissipation principle   [46] which gives rise to a constrained extremization problem that we solved using standard derivative-free optimization methods at every time step. In configurations similar to the one shown in Fig. 6 the system   would   undergo   a   simultaneous   impact   across   two manifolds: the tip of the tail and the tip of one of the legs. For these configurations we can calculate the angle between the two manifolds, and, according to Sec. V, if this angle is  π/ 2   the indeterminacy of the outcome will be zero. However, picking a random configuration and set of parameters such that the two contacts are established will, most likely, return a non- orthogonal pair of manifolds. This is not hard to imagine, since the two manifolds themselves depend on the configuration and 10
Initial   Optimized Fig. 7: The running robot model (TREX) that we used in order to illustrate the results of this paper. Both systems are presented in a stance in which double impact is occurring. The model on the right has parameters chosen by hand as being reasonable. Its configuration and design parameters act as an initial condition for an optimization to generate the system and configuration on the left. The exact parameters and configuration variables for the two systems are also shown (see Table I).  their dot product is defined through the kinetic energy metric, which also depends on the inverse of the mass matrix, and hence on the configuration. We assume that most of the energy lost through impact comes from dissipation in the knee joints and that a relatively small amount goes into permanent deformation of either the robot or the ground. Thus we model all impacts as elastic and expect some degree of exponentially decaying chatter- ing. When the chattering becomes faster than the time step frequency, we assume we have reached Zeno behavior and consider that impact plastic, in effect taking away the rest of the energy that would be lost through very high frequency— and probably not modelable—motion in the dampers. In order to test that our results presented in this paper would have usefulness in gait and mechanism design, we performed the following simulation. First, we chose a random configuration of the robot at the time of simultaneous impact, such that both the tail and the right foot were touching the floor simultaneously. We made sure that the contact manifolds were not orthogonal under the kinetic energy metric and labeled the configuration along with the design parameters to be the   unoptimized system . Next, we used a classic root finding algorithm for underconstrained systems in order to find a nearby configuration and set of parameters for which the dot product between the two manifolds of contact is zero: this is the   optimized system . The two systems are presented in Fig. 7 and their parameters can be seen in Table I. Furthermore, we also considered a system with parameters identical to the optimized system but with a   initial stance : both the legs were straight at the knee and   180 ◦   out of phase at the hips. For both the optimized gait and the initial gait we raised the mechanism   17 cm   from the ground and used the PD controller at   the   hips   to   keep   the   robot   in   the   same   configuration until the moment of impact—at which time the controller applied torque, spinning the legs counterclockwise in order to create forward movement. During this first interaction several simultaneous impacts occur, and we have a choice of which manifold to solve for first: the leg or the tail. In the first run through we chose the index order for the contact manifolds based on the   arg min   〈 w i ,   p +  ∗   〉 g   . On the second run through, however, we reversed this order, effectively using   arg max  instead of   arg min   for every impact. The parameters that we Time after first collision ( ms ) 100 0   20   40   60   80  10 − 5  10 − 10  Outcome indeterminacy   ξ   (%)  ξ   =   ‖ p −  ̃ p ‖ g  ‖ p ‖ g  10 0  initial gait optimized gait  Fig. 8: A coordinate-independent measure of the indeterminacy in momentum is plotted for two gaits of the optimized mechanism: the optimized gait, and the initial gait, in which both legs are straight at the knee and 180 ◦   out of phase with each other. The measure of indeterminacy was chosen in an analogous way to that in Sec. VII-A. We solved simultaneous impacts in each copy of the mechanism by choosing a different order of reflections. The normed distance between the momenta of the two mechanism versions at each time point is then normalized by the kinetic energy at that time.  considered and their values for both the unoptimized and optimized system are shown in Table I. The difference in behavior between the optimized and the initial case is presented in Fig. 8 which shows an indetermi- nacy measure for each system based on the normalized uncer- tainty of the momentum over time. The indeterminacy for the optimized gait is more than five orders of magnitude smaller than that of the initial gait, suggestive of an improvement in the modelability of the simultaneous impact. VIII. C ONCLUSION  We have shown how the geometry and configuration of a rigid body system at the time of simultaneous impact affects its sensitivity to initial conditions, when a propagative rigid body impact model is used. The measure of sensitivity is obtained from the inner product between the contact manifold normals at the impact configuration, and is related to the uniqueness of solutions under the propagative model—existence of solutions was also shown for the case of two simultaneous impacts. We optimized two example systems—one analytically, the other numerically—to minimize the sensitivity during a simultane- ous impact. Both optimized and non-optimized versions of the numerical model were simulated using a time stepping scheme based on variational integrators, and significant sensitivity improvement was shown between the two cases. A CKNOWLEDGMENTS  This material is based upon work supported by the National Science Foundation under grant IIS-1018167. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily re- flect the views of the National Science Foundation. We would also like to acknowledge Dr. David Pekarek for the invaluable help and support given through fruitful conversations relating to problems common of all collision methods, as well as Dr. Goce Trajcevski for discussions regarding the existence proof. 11
TABLE I: Design parameters, configuration variables, and their values for both the optimized and unoptimized TREX systems. The largest change occurs in  φ R , which is the right hip angle. Configuration variables + Rest angles for torsional springs   Design parameters System   x   y   θ   φ L   φ R   γ L   γ R   φ T   r 1   r 2   l u   l l   〈 u ,   v 〉 g  ( cm )   ( cm )   ( rad )   ( rad )   ( rad )   ( rad )   ( rad )   ( rad )   ( cm )   ( cm )   ( cm )   ( cm )  Initial   0   12.38   5.78   5.13   1.99   5.24   5.24   0   7   4   8   10   1 . 65   ×   10 − 3  Optimized   0   13.69   5.92   4.79   2.85   5.71   5.67   -0.19   6.69   4.22   7.81   9.86   − 2 . 28   ×   10 − 7  R EFERENCES [1]   D.   Baraff,   “Analytical   methods   for   dynamic   simulation   of   non- penetrating rigid bodies,”   SIGGRAPH Comput. Graph. , vol. 23, no. 3, pp. 223–232, Jul. 1989. [2]   A. Chatterjee and A. Ruina, “A new algebraic rigid-body collision law based on impulse space considerations,”   ASME J. Appl. Mechanics , vol. 65, no. 4, pp. 939–950, 1998. [3]   M. Eich, F. Grimminger, and F. Kirchner, “A versatile stair-climbing robot for search and rescue applications,” in   IEEE International Work- shop on Safety, Security and Rescue Robotics, 2008. SSRR 2008 , Oct. 2008, pp. 35 –40. [4]   R. D. Quinn, G. M. Nelson, R. J. Bachmann, D. A. Kingsley, J. T. Offi, T. J. Allen, and R. E. Ritzmann, “Parallel complementary strategies for implementing biological principles into mobile robots,”   The Inter- national Journal of Robotics Research , vol. 22, no. 3-4, pp. 169–186, Mar. 2003. [5]   J. B. Jeans and D. Hong, “IMPASS: intelligent mobility platform with active spoke system,” in   IEEE International Conference on Robotics and Automation, 2009. ICRA ’09 , May 2009, pp. 1605 –1606. [6]   D. Lyons and K. Pamnany, “Rotational legged locomotion,” in   , 12th International Conference on Advanced Robotics, 2005. ICAR ’05. Pro- ceedings , Jul. 2005, pp. 223 –228. [7]   U. Saranli, M. Buehler, and D. Koditschek, “RHex: a simple and highly mobile hexapod robot,”   Int. J. Robotics Research , vol. 20, no. 7, pp. 616–631, Jul. 2001. [8]   Y. Xiong and X. Xiong, “Algebraic structure and geometric interpreta- tion of rigid complex fixture systems,”   IEEE Trans. Autom. Sci. Eng. , vol. 4, no. 2, pp. 252 –264, Apr. 2007. [9]   E. Staffetti, “Anal. of rigid body interactions for compliant motion tasks using the grassmann-cayley algebra,”   IEEE Trans. Autom. Sci. Eng. , vol. 6, no. 1, pp. 80 –93, Jan. 2009. [10]   R. Primerano, D. Wilkie, and W. Regli, “A case study in system-level physics-based simulation of a biomimetic robot,”   IEEE Trans. Autom. Sci. Eng. , vol. 8, no. 3, pp. 664 –671, Jul. 2011. [11]   S. Berard, J. Trinkle, B. Nguyen, B. Roghani, J. Fink, and V. Kumar, “daVinci code: A multi-model simulation and anal. tool for multi-body systems,” in   Proc. IEEE Int. Conf. Robotics and Automation , Apr. 2007, pp. 2588–2593. [12]   D. Stewart and J. C. Trinkle, “An implicit time-stepping scheme for rigid body dynamics with coulomb friction,” in   Proc. IEEE Int. Conf. Robotics and Automation , vol. 1, 2000, pp. 162–169. [13]   N. Chakraborty, S. Berard, S. Akella, and J. Trinkle, “An implicit time- stepping method for multibody systems with intermittent contact,” in  Robotics: Science and Systems , Jun. 2007. [14]   J. Trinkle, J.-S. Pang, S. Sudarsky, and G. Lo, “On dynamic multi-rigid- body contact problems with Coulomb friction,”   ZAMMJ. Appl. Math. and Mechanics/Zeitschrift f ̈ ur Angewandte Mathematik und Mechanik , vol. 77, no. 4, pp. 267–279, 1997. [15]   J. J. Moreau and M. Jean, “Numerical treatment of contact and friction: the contact dynamics method,”   ESDA 1996: Photomechanics; Contact mechanics and tribology , vol. 4, pp. 201–208, 1996. [16]   J. Moreau, “Numerical aspects of the sweeping process,”   Comput. Methods Appl. Mechanics and Eng. , vol. 177, no. 3-4, pp. 329–349, Jul. 1999. [17]   T. Liu and M. Y. Wang, “Computation of three-dimensional rigid- body dynamics with multiple unilateral contacts using time-stepping and gauss-seidel methods,”   IEEE Autom. Sci. Eng. , vol. 2, no. 1, pp. 19–31, Jan. 2005. [18]   “Newton.” [Online]. Available: http://newtondynamics.com/ [19]   “PhysX.” [Online]. Available: http://www.nvidia.com/ [20]   “ODE: open dynamics engine.” [Online]. Available: http://www.ode.org/ [21]   “Bullet   physics   library,”   Sep.   2011.   [Online].   Available:   http: //bulletphysics.org/ [22]   C. M. Donahue, C. M. Hrenya, A. P. Zelinskaya, and K. J. Nakagawa, “Newton’s cradle undone: Experiments and collision models for the normal collision of three solid spheres,”   Physics of Fluids , vol. 20, no. 11, p. 113301, 2008. [23]   B. V. Mirtich, “Impulse based dynamic simulation of rigid body sys- tems,” Dissertation, University of California at Berkeley, 1996. [24]   J. K. Hahn, “Realistic animation of rigid bodies,”   SIGGRAPH Comput. Graph. , vol. 22, no. 4, p. 299–308, Jun. 1988. [25]   Y.   Jia,   M.   Mason,   and   M.   Erdmann,   “Simultaneous   impacts: a   state   transition   diagram   approach,”   2011.   [Online].   Available: http://www.cs.iastate.edu/ ∼ jia/papers/IJRR11b-submit.pdf [26]   Y.-B. Jia, M. Mason, and M. Erdmann, “A state transition diagram for simultaneous collisions with application in billiard shooting,” in  Algorithmic Found. of Robotics VIII , B. Siciliano, O. Khatib, F. Groen, G. S. Chirikjian, H. Choset, M. Morales, and T. Murphey, Eds.   Springer Berlin Heidelberg, 2009, vol. 57, pp. 135–150. [27]   Z. Zhao, C. Liu, and B. Brogliato, “Energy dissipation and dispersion effects in granular media,”   Physical Review E , vol. 78, no. 3, p. 031307, Sep. 2008. [28]   C. Liu, Z. Zhao, and B. Brogliato, “Frictionless multiple impacts in multibody systems. i. theoretical framework,”   Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science , vol. 464, no. 2100, pp. 3193–3211, Dec. 2008. [29]   ——, “Frictionless multiple impacts in multibody systems. II. numerical algorithm and simulation results,”   Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science , vol. 465, no. 2101, pp. 1–23, Jan. 2009. [30]   V. Seghete and T. Murphey, “Variational solutions to simultaneous collisions between multiple rigid bodies,” in   Proc. IEEE Int. Conf. Robotics and Automation , May 2010, pp. 2731–2738. [31]   ——, “Conditions for uniqueness in simultaneous impact with applica- tion to mechanical design,” presented at the IEEE Int. Conf. for Robotics and Automation, St. Paul, MN, 2012. [32]   V. Seghete and T. D. Murphey, “Multiple instantaneous collisions in a variational framework,” in   Proc. IEEE Conf. Decision and Control , 2009, pp. 5015–5020. [33]   J. E. Marsden and T. S. Ratiu,   Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems , 2nd ed.   Springer Verlag, 1999. [34]   F. Bullo and A. D. Lewis, “Section 4.2 the kinetic energy metric,” in   Geometric Control of Mechanical Systems: Modeling, Analysis, and Design for Simple Mechanical Control Systems , ser. Texts in Appl. Math. New York: Springer Sci.+Bus. Media Inc., 2010, no. 49, pp. 162–172. [35]   B. Blazejczyk-Okolewska, K. Czolczynski, T. Kapitaniak, and J. Wo- jewoda, “Section 1.5. coefficient of restitution,” in   Chaotic mechanics in systems with impacts and friction , ser. World Scientific Series on Nonlinear Sci.   World Scientific Publishing Co. Pte. Ltd., 1999, vol. 36, no. A, pp. 55–59. [36]   E. Johnson and T. Murphey, “Scalable variational integrators for con- strained mechanical systems in generalized coordinates,”   IEEE Trans. Robot. , vol. 25, no. 6, pp. 1249–1261, Dec. 2009. [37]   R. C. Fetecau, J. E. Marsden, M. Ortiz, and M. West, “Nonsmooth lagrangian mechanics and variational collision integrators,”   SIAM J. Appl. Dynamical Syst. , vol. 2, no. 3, pp. 381–416, 2003. [38]   A. Lew, J. Marsden, M. Ortiz, and M. West, “An overview of variational integrators,” in   Finite Element Methods: 1970’s and Beyond , L. Franca, Ed., CIMNE, Barcelona, Spain, 2003, pp. 1–18. [39]   ——, “Variational time integrators,”   Int. J. Numerical Methods Eng. , vol. 60, no. 1, p. 153–212, May 2004. [40]   C. Kane, J. E. Marsden, and M. Ortiz, “Symplectic-energy-momentum preserving variational integrators,”   J. Math. Physics , vol. 40, pp. 3353– 3371, 1999. [41]   O. Junge, J. Marsden, and S. Ober-Bl ̈ obaum, “Discrete mechanics and optimal control,” in   Proc. IFAC World Congress .   Czech Republic: IFAC, Jul. 2005, pp. 744–744.  12
[42]   S. Leyendecker, S. Ober-Bl ̈ obaum, J. Marsden, and M. Ortiz, “Discrete mechanics and optimal control for constrained multibody dynamics,” in  Proc. IDETC/MSNDC .   Las Vegas, USA: ASME, 2007, pp. 4–7. [43]   D. N. Pekarek, “Variational methods for control and design of bipedal robot models,” Thesis, California Inst. of Technology, Pasadena CA, 2010. [44]   D. N. Pekarek and J. E. Marsden, “Variational collision integrators and optimal control,”   Proc. 18th Int. Symp. Math. Theory Networks and Syst. (MTNS) , 2008. [45]   “Wrexx.”   [Online].   Available:   http://www.iloverobots.com/ shop-robotic-toys/mechatars/wrexx-product-set [46]   D. E. Stewart, “Rigid-body dynamics with friction and impact,”   SIAM review , pp. 3–39, 2000.  A PPENDIX  Lemma 1:   Given a sequence   { w i }   which is minimal with respect to some momentum   p − , we must have that   w j   6   =  w j +1   for all   j .  Proof:   We prove this by reductio ad absurdum: suppose there exists a   j   such that   w j   =   w j +1 . We have then that  Γ( w j   )Γ( w j +1 ) = Γ 2 ( w j   ) =   I . This means that we can write:  p +   =   p −  n ∏  i =1  Γ( w i ) =   p −  j − 1 ∏  i =1  Γ( w i )  n ∏  k = j +2  Γ( w k ) ,  which   implies   that   the   subsequence   { w k }   of   { w i }   with elements   w j   and   w j +1   removed also generates a feasible momentum. Thus, the sequence   { w i }   cannot be minimal.  Lemma 2:   Given   p   ∈   T   ∗  q ∗   Q   and a sequence of trans- formations   Γ( w i )   with   w i   ∈ { u ,   v }   that map   p   to   p f   =  p   ∏  i   Γ( w i ) , we always have that  p f   =   P T ( p ) +   P N ( p )   ∏  i  Γ( w i ) ,  Proof:   It follows from their definitions that   T   and   N   are complementary, such that   T   ×   N   =   T   ∗  q ∗   Q.   Thus, we can write  p f   =   P T ( p f   ) +   P N ( p f   ) .  Let us look at the first term  P T ( p f   ) =   P T  [  p   ∏  i  Γ( w i )  ]  =   P T  {  [ P T ( p ) +   P N ( p )]   ∏  i  Γ( w i )  }  .   (32) From the definition of   T   and (7) we have that all transfor- mations that reflect across either the   u   or the   v   plane leave vectors in   T   unchanged. Thus, we can say that  P T ( p )   ∏  i  Γ( w i ) =   P T ( p )   ∈   T .   (33) Similarly, we have that  P N ( p )   ∏  i  Γ( w i ) =   P N  [  p   ∏  i  Γ( w i )  ]  ∈   N .   (34) Substituting (33) and (34) into the right hand side of (32) and using the orthogonality of   T   and   N , we obtain  P T ( p f   ) =   P T ( p )   (35) An identical argument gives the following result for   P N ( p f   ) :  P N ( p f   ) =   P N ( p )   ∏  i  Γ( w i ) .   (36) The statement of the lemma follows directly from (35) and (36).  Lemma 3:   Let   p ,   u ,   v ,   r   ∈   N   such that   r   =   u + v  ‖ u + v ‖ g  and  u   6   =   ± v . Let   γ   = arcsin   〈 r ,   u 〉 g   . Given the above we have that   p   is feasible   iff  〈 p ,   r 〉 g   ≥ ‖ p ‖ g   cos   γ.   (37)  Proof:   We start by noting that when   γ   6   = 0 —this is true since we assume   u   and   v   are not collinear—several useful relations hold:  〈 u ,   v 〉 g   =   −   cos 2 γ,  ‖ u   +   v ‖ g   = 2 sin   γ,  and  〈 p ,   r 〉 g   =   〈 p ,   u 〉 g   +   〈 p ,   v 〉 g  2 sin   γ   .   (38) We   first   prove   the   forward   implication.   Note   that   the feasibility conditions together with (38) already give us a lower bound by guaranteeing that   〈 p ,   r 〉 g   ≥   0 . Since   p   ∈  span { u ,   v }   we can write  p   =   a u   +   b v ,   a, b   ∈   R .  The norm of   p   can be expressed as  ‖ p ‖ g   =   √ a 2   −   2 ab   cos 2 γ   +   b 2 ,  and the feasibility conditions become  a   −   b   cos 2 γ   ≥   0 , b   −   a   cos 2 γ   ≥   0 ,  which, in turn, imply  ab   sin 2   2 γ   −   2   ( a 2   +   b 2 )   cos 2 γ   ≥   0 .  Several trigonometric manipulations show that the previous is equivalent to  ( a   +   b ) 2   sin 2   γ   ≥   cos 2   γ   ( a 2   −   2 ab   cos 2 γ   +   b 2 )   .  The right hand side can be written in terms of the norm of   p :  ( a   +   b ) 2   sin 2   γ   ≥   cos 2   γ   ‖ p ‖ 2  g   .  We can also write  〈 p ,   u 〉 g   +   〈 p ,   v 〉 g   = 2( a   +   b ) sin 2   γ,  which, after squaring, gives  (  〈 p ,   u 〉 g   +   〈 p ,   v 〉 g  ) 2  4 sin 2   γ   = ( a   +   b ) 2   sin 2   γ.  Putting it together, we have  〈 p ,   r 〉 2  g   ≥   cos 2   γ   ‖ p ‖ 2  g   ,  which proves the forward implication in our lemma. 13
Now, for the reverse implication, we assume (37) and want to show that   p   is feasible. Following the last few step of the forward proof, we can show that (37) is equivalent to  〈 p ,   u 〉 g   +   〈 p ,   v 〉 g   ≥ ‖ p ‖ g   sin 2 γ.  We want to show that both   〈 p ,   u 〉 g   and   〈 p ,   v 〉 g   are positive. Let   u t   =   u   − 〈 u ,   v 〉 g   v   be the projection of   u   onto the plane normal to   v . Using this notation, we write  〈 p ,   v 〉 g   +   〈 p ,   u t 〉 g   +   〈 p ,   v 〉 g   〈 u ,   v 〉 g   ≥ ‖ p ‖ g   ‖ u t ‖ g   ,  which, after some algebra, becomes  〈 p ,   v 〉 g   ≥ ‖ p ‖ g   ‖ u t ‖ g   − 〈 p ,   u t 〉 g  2 sin 2   γ   ≥   0 ,  where the last inequality holds due to the fact that an inner product is always less than the product of the vector norms. Using an identical argument we show that   〈 p ,   u 〉 g   ≥   0 . The two statements together are equivalent to feasibility.  Lemma 4:   For any minimal sequence   { w i }   we have that  〈 r i ,   r 0 〉 g   = cos(2 iγ ) ,  where, as in Lemma 4,   γ   = arcsin   〈 r 0 ,   u 〉 g   = arcsin   〈 r 0 ,   v 〉 g  and   r i   =   r 0  ∏ i j =0   Γ( w j   ) .  Proof:   We show this using mathematical induction and proving that, when   w 0   =   u  〈 r k ,   u 〉 g   =  {   sin[(2 k   + 1) γ ]   :   w k − 1   =   v  −   sin[(2 k   −   1) γ ]   :   w k − 1   =   u   (39a)  〈 r k ,   v 〉 g   =  {   −   sin[(2 k   −   1) γ ]   :   w k − 1   =   v  sin[(2 k   + 1) γ ]   :   w k − 1   =   u   (39b)  〈 r k ,   r 0 〉 g   = cos(2 kγ )   (39c) for   ∀ k >   0 . A symmetric result holds when   w 0   =   v , such that (39c) remains unchanged. For   k   = 0 , before any reflections are applied, we have that  〈 r 0 ,   u 〉 g   = sin( γ ) ,  〈 r 0 ,   v 〉 g   = sin( γ ) ,  〈 r 0 ,   r 0 〉 g   =   ‖ r 0 ‖ g   = 1 .  For   k   = 1 , we make use of the consequence of Lemma 1 by observing that   { w i }   must consist of alternating elements. This means that   w 1   =   u   and   r 1   =   r 0 Γ( u ) . We calculate  〈 r 1 ,   u 〉 g   =   〈 r 0 Γ( u ) ,   u 〉 g   =   〈 r 0 ,   u 〉 g   −   2   〈 r 0 ,   u 〉 g   ‖ u ‖ 2  g  = sin( γ )   −   2 sin( γ ) =   −   sin( γ ) ,  〈 r 1 ,   v 〉 g   =   〈 r 0 Γ( u ) ,   v 〉 g   =   〈 r 0 ,   v 〉 g   −   2   〈 r 0 ,   u 〉 g   〈 u ,   v 〉 g  = sin( γ ) + 2 sin( γ ) cos(2 γ ) = sin(3 γ ) ,  〈 r 1 ,   r 0 〉 g   =   〈 r 0 Γ( u ) ,   r 0 〉 g   =   ‖ r 0 ‖ g   −   2   〈 r 0 ,   u 〉 g   〈 u ,   r 0 〉 g  = 1   −   2 sin 2 ( γ ) = cos(2 γ ) ,  thus showing the first step of the induction proof when   w 0   =  u . The argument is symmetrical for the case when   w 0   =   v . Next, we assume that (39) holds for   k   and show that it also holds for   k   + 1   under the assumption that   w k − 1   =   u   . We start with  〈 r k +1 ,   u 〉 g   =   〈 r k ,   u 〉 g   − 2   〈 r k ,   u 〉 g   ‖ u ‖ g   =   −   sin[(2 k   +1) γ ] =   −   sin { [2( k   + 1)   −   1] γ } .  We then show that  〈 r k +1 ,   v 〉 g   =   〈 r k ,   v 〉 g   −   2   〈 r k ,   u 〉 g   〈 u ,   v 〉 g  =   −   sin[(2 k   −   1) γ ] + 2 sin[(2 k   + 1) γ ] cos(2 γ ) = sin { [2( k   + 1) + 1] γ } .  Finally, we have  〈 r k +1 ,   r 0 〉 g   =   〈 r k ,   r 0 〉 g   −   2   〈 r k ,   u 〉 g   〈 u ,   r 0 〉 g  = cos(2 kγ )   −   2 sin[(2 k   + 1) γ ] sin( γ ) = cos[2( k   + 1) γ ] .  A symmetric argument holds under the complementary as- sumption that   w k − 1   =   v . Thus, our inductive proof is finished, and we have showed that  〈 r n ,   r 0 〉 g   = cos(2 nγ ) .   (40)  Lemma 5:   For two contact manifolds described by their normals   u   6   =   ± v , and any infeasible momentum   p   with  〈 u ,   p 〉 g   ≤ 〈 v ,   p 〉 g   <   0 , we have that  p f   =   p Γ ( u ) Γ ( v ) =   p Γ ( v ) Γ ( u )   is feasible   (41)  iff  〈 u ,   v 〉 g   = 0 .   (42)  Proof:   Suppose that (41) holds. We then have that  〈 p ,   u 〉 g   〈 u ,   v 〉 g   v   =   〈 p ,   v 〉 g   〈 u ,   v 〉 g   u ,   ∀ p   ∈   T   ∗  q ∗   Q.  Since we assume that   u   6   =   ± v , the only way in which (A) will hold for any   p   is if  〈 u ,   v 〉 g   = 0 .  Conversely, assuming (42), we can write  p f   =   p Γ( u )Γ( v ) =   p Γ( u )   −   2   〈 p Γ( u ) ,   v 〉 g   v  =   p   −   2   〈 p ,   u 〉 g   u   −   2   〈 p ,   v 〉 g   v   + 4   〈 p ,   u 〉 g   〈 u ,   v 〉 g   v  =   p   −   2  (  〈 p ,   u 〉 g   u   +   〈 p ,   v 〉 g   v  )  =   p Γ( v )Γ( u ) .  The symmetry of the result in   u   and   v   assures us of the commutativity of the   Γ( u )   and   Γ( v ) . The feasibility of   p f  is given by the following:  〈 p f   ,   u 〉 g   =   〈 p ,   u 〉 g   −   2   〈 p ,   u 〉 g   =   − 〈 p ,   u 〉 g   ≥   0 ,  〈 p f   ,   v 〉 g   =   〈 p ,   v 〉 g   −   2   〈 p ,   v 〉 g   =   − 〈 p ,   v 〉 g   ≥   0 ,  Thus, the statement has been proven. 14