On Dynamic Multi-Rigid-Body Contact Problems

with Coulomb Friction

\Lambda 

Jeff Trinkley Jong-Shi Pangz Sandra Sudarskyx Grace Lo-

October 3, 1995

Abstract. This paper is summary of a comprehensive study of the problem of predicting the possible acceleration(s) of a set of rigid, three-dimensional bodies in contact in the presence of Coulomb friction. We begin with a brief introduction to this problem and a survey of related work and previous approaches. This is followed by the introduction of two novel complementarity formulations for the contact problem under two friction laws: Coulomb's Law and an analogous law in which Coulomb's quadratic friction cone is approximated by a pyramid. Under a full column rank assumption on the system Jacobian matrix, we establish the existence and uniqueness of a solution to our new models in the case where the friction coefficients are nonnegative and sufficiently small. For the model based on the friction pyramid law, we also show that the classical Lemke almost-complementary pivot algorithm and our new feasible interior point method are guaranteed to compute a solution. Extensive computational results are presented to demonstrate the effectiveness and characteristics of these solution methods for solving the model with friction pyramids. Finally, results are summarized and conclusions are drawn.

Key Words. Multi-rigid-body contact problem, Coulomb friction, complementarity problem, quasi-variational inequality, Lemke's algorithm, interior point algorithm.

\Lambda Any findings, conclusions, or recommendations expressed herein are those of the authors and do not necessarily reflect the views of the funding agencies.

yThe research of this author was supported by the National Science Foundation under grant IRI-9304734, the

Texas Advanced Research Program under grant 999903-078, and the Texas Advanced Technology Program under grant 999903-095. Department of Computer Science, Texas A&M University, College Station, Texas 77843-3112, U.S.A. Email: trink@cs.tamu.edu.

zThe research of this author was based on work supported by the National Science Foundation under grant CCR9213739 and the Office of Naval Research under grant N00014-93-1-0228. Department of Mathematical Sciences, The Johns Hopkins University, Baltimore, Maryland 21218-2689, U.S.A. Email: jsp@vicp.mts.jhu.edu

xThe research of this author was supported by the Texas Advanced Research Program under grant 999903-

078. Department of Computer Science, Texas A&M University, College Station, Texas 77843-3112, U.S.A. Email: sudarsky@cs.tamu.edu

-Department of Mathematical Sciences, The Johns Hopkins University, Baltimore, Maryland 21218-2689, U.S.A.

Email: lo@brutus.mts.jhu.edu

1

1 Introduction The dynamic multi-rigid-body contact problem is concerned with predicting the accelerations and contact forces of several rigid bodies in contact; see Figure 1. Armed with a mathematical model and a solution method for this problem, one can study a wide variety of complex systems. For example, one could simulate the flow of granular material in a hopper or develop a system to automatically plan mechanical assembly tasks. Our new complementarity formulations of the dynamic multi-rigidbody contact problem are applicable to three-dimensional systems analogous to the two-dimensional one shown in Figure 1. There are a number of passive, three-dimensional, bodies referred to as

Figure 1: Multi-Rigid-Body System objects (shown grey) that move in response to external forces (e.g., gravitational forces) and contact with a number of active bodies referred to as manipulator links (shown white). The hatched lines on the bottom link indicate that it is fixed with respect to an inertial frame.

The relevant equations and constraints are the Newton-Euler equations, unilateral and bilateral kinematic constraints, and a friction law. The bilateral (equality) constraints arise from the joint connections between the links of the active mechanism. These constraints reduce the number of degrees of freedom of the mechanism and allow one to define a small set of generalized coordinates to describe its configuration. The unilateral (inequality) constraints disallow interpenetration of the bodies at their contact points. One other type of constraint that appears in such systems is known as a complementarity constraint. Such a constraint arises between the normal component of a contact force and the normal component of relative acceleration at the same contact. These complementarity constraints ensure that either the normal component of the relative acceleration at the contact or the normal component of the contact force is zero. Physically, when the normal component of the relative acceleration is positive (i.e., the contact is breaking), the normal component of the contact force must be zero. Conversely, when the normal component of the contact force is positive, the normal component of the relative acceleration must be zero.

A common approach to solving the multi-rigid-body contact problem has been formalized and implemented in the DADS software package by Haug [10]. In this approach, one formulates the equations of motion as a system of differential algebraic equations (DAE's). To simulate the system,

2

the differential equations are integrated on the constraint "surface" represented by the (algebraic) contact constraints, while contact forces and velocities are monitored to determine when to change the assumed contact interactions (e.g., rolling, sliding, and breaking). Collision detection is used to determine when to apply impact laws.

In more detail, the basic idea is to assume that the contact interactions are known, so that the unilateral kinematic constraints can be replaced by a set of bilateral constraints. In other words, if a pair of bodies touch and the normal component of the contact force is compressive at the current time, then the bodies are assumed to remain in contact during the next integration time step. This is enforced by replacing the fundamental inequality constraint that allows the possibility of positive relative acceleration in the normal direction (indicating the breaking of the contact) by an equation requiring the normal acceleration to be zero. Similar arguments are applied to generate equality constraints for contacts which are expected to roll.

While the DAE approach to solving the multi-rigid-body contact problem seems to work well in practice, it has two fundamental flaws. It is based on an implicit assumption that the contact forces are continuous functions of time (between impact events) and an assumption that one can always deduce the impending contact interactions (e.g., rolling, sliding, and breaking). While the first assumption is reasonable for real systems, it is not justified mathematically when the bodies are assumed to be rigid. The second assumption is more problematic, because immediately following an impact, the finite force (history) information needed to deduce the contact interactions is not available.

For example, consider a rod that falls onto a table top and assume that one of the tips of the rod touches the table first. The equations of motion can easily be integrated until the time of impact. Then an impact model can be applied to yield a post-impact velocity for the rod. From this result, one can determine initial conditions for the first integration step after the impact. However, it is not necessarily possible to determine the correct set of kinematic constraints (needed by the DAE method) for this time step.

Continuing with the example, let us assume that just after impact, the tip of the rod is rolling. The DAE method would formulate the equations of motion under the assumption that the tip of the rod continues to roll, and integration would continue. The problem with this is that L"otstedt determined initial conditions under which the rod's ensuing motion cannot be determined uniquely [15]. Indeed, the rod's tip could begin to slide or separate, or it could continue to roll. Thus, unbeknownst to the DAE user, one of three possible evolutions of the rod has been determined. While it is certainly a deficiency of the model that it can become ambiguous, perhaps a bigger problem is that an unaware user of the DAE approach would take the result as the unique system evolution.

As will be discussed below, multi-body contact problems arise in many engineering applications in which it would be desirable to ignore deformations of the bodies. Thus despite its deficiencies, the multi-rigid-body contact problem is quite important. The primary aim of this paper is to improve our understanding by formulating and studying a mathematical model which retains all the fundamental unilateral and complementarity constraints.

3

1.1 Further background and motivation Henceforth, "model" will refer to the fundamental multi-rigid-body model without the assumptions used in the DAE approach.

In simple situations, such as a block on an inclined plane, the multi-rigid-body model often predicts system motions that are intuitive. However, it can also produce quite unexpected results. For example, consider a rigid object sliding at one point on another. It has been known for at least 100 years that if Coulomb friction exists at the contact, then the rigid-body model will not always possess a solution [18]. To make matters worse, it is possible that the model admits several distinctly different solutions for one set of initial conditions and input. As mentioned earlier,

L"otstedt constructed an example in which the model admitted three solutions [15].

Despite this difficulty, multi-rigid-body models have recently received great interest in the virtual reality (VR), graphics, and robotics research communities. The VR community is motivated by the desire to enhance the believability of the scenes in which the VR user is submersed. Currently, systems cannot handle situations involving the loss of contact and multiple contacts between objects. For example, in some of the most advanced systems, if the user releases his virtual grasp of an object, it will hover as though weightless [19]. A thorough understanding of the multi-rigid-body contact problem could lead to efficient solution algorithms that could be embedded in VR systems, alleviating this problem and related ones.

The robotics and graphics communities are interested in the inverse problem. More specifically, roboticists and animators want to discover input trajectories that accomplish a task or action described at a very high level. The problems to be solved can be formulated as optimal control problems in which the initial and final states of the system are specified along with an objective to optimize. This objective may be to minimize the energy used by an assembly robot or to make an animated creature move gracefully. These optimal control problems are difficult to solve even using the DAE approach, because every time a new contact is formed or the interaction at a contact changes, the equations of motion change structurally. Currently, the most capable animation systems require the animator to tediously specify times at which contacts should occur, slide, roll, and break [13]. Analogous robot manipulation planning systems exist in which a contact mode1 or sequence of contact modes is specified as input, after which, the planner attempts to find an appropriate trajectory [3, 16, 23]. Since these constraints are specified by the user artificially, such planning systems cannot always find solutions when they exist.

Simple robot task planners and animators' tools could be developed without the need for user specification of contact modes, and some have been developed that can successfully plan simple contact tasks [4, 8, 26, 28]. However, their development has been significantly hampered by a lack of understanding of rigid body contact mechanics.

In an effort to improve our understanding, we recently formulated the three-dimensional, multirigid-body contact problem with Coulomb friction as a nonlinear complementarity problem (NCP) [6, 21]. This formulation allowed us to apply existing complementarity theory and the theory of quasi-variational inequalities to prove existence results that extended and subsumed all previous results known to us. In addition, our study provided new insight into why the model sometimes

1A contact mode is a set of contacts and their interactions, i.e., rolling, sliding, or breaking.

4

failed to have a solution, among other things. This led us to propose and study a new model based on a relaxed friction law that always has a solution. Our results yielded strict theoretical guidelines for the applicability of the models, and thus are of great importance in the fields of robotics, graphics, and virtual reality, and more generally in the field of rigid body mechanics.

The specific objectives of this paper are to introduce to the engineering community a set of mathematical tools for the study, analysis, and numerical solution of the three-dimensional, multirigid-body contact problem with Coulomb friction. We will motivate and develop the models, extend the theoretical results established in [21] and discuss their practical significance, present applicable numerical methods, and discuss our computational experience with these methods.

1.2 Related work A derivation of the equations of motion and the relevant contact constraints for planar multirigid-body systems under Coulomb's friction law was given by L"otstedt in [14]. While L"otstedt recognized and exploited the complementarity relationships between the normal components of the contact forces and the relative contact accelerations, he did not formally cast the Coulomb friction constraints at the rolling contacts as standard complementarity relationships, and therefore, he did not achieve a linear complementarity formulation for any problem with one or more rolling contacts. He also did not treat the three-dimensional problem, which turns out to be considerably more complicated, because it no longer corresponds to a linear complementarity problem. Nevertheless, by approximating Coulomb's quadratic friction cone as a four-sided pyramid and introducing appropriate slack variables, we are able to formulate the model for three-dimensional systems as a standard linear complementarity problem (LCP). We also show that the problem with quadratic friction cones cannot be conveniently formulated as a standard NCP, but that it can be formulated as a special type of nonlinear complementarity problem that has been previously studied.

Rajan, Burridge and Schwartz [24] discussed the planar motion of an object initially at rest and in frictional contact with one or two immobile rigid walls. The authors partitioned the force/torque space into several regions, each of which corresponded to a possible contact mode. Any situation in which two or more of these regions overlapped, indicated that several different solutions were simultaneously possible. For the case of a single contact point, the authors derived a condition for which three different contact modes are possible. These conditions can be shown to be equivalent to those derived by L"otstedt for the rod-and-table problem discussed earlier [14] and, once the problem is formulated as an LCP, the conditions can be derived from LCP theory.

Mason and Wang [17], studied L"otstedt's sliding rigid rod problem and derived a condition for inconsistent configurations which is a special case of the condition stated in [14]. In particular, when the rod leans opposite its sliding direction and the contact force is restricted to have finite magnitude (as suggested by Kilmister and Reeve [11]), it can happen that the only possible solution is one in which the rod would penetrate the surface on which it is sliding; thus the inconsistency. Mason and Wang showed that by allowing the force to be impulsive (i.e., of infinite magnitude), the problem could be solved. They argued that there is no physical basis for not allowing impulsive solutions when bodies are rigid. While this stance allowed the reclassification of the "inconsistent" sliding rod problem as consistent, Mason and Wang noted the possibility that problems with unique finite force solutions may have multiple solutions when impulsive forces are also allowed.

5

Motivated by animation applications, Baraff has developed solution algorithms for multi-rigidbody contact problems [2]. From L"otstedt's work it was clear that finding finite force solutions (if they exist) to these problems can be computationally complex. Therefore, Baraff [2] developed an algorithm similar to the Cottle-Dantzig principal pivoting algorithm [6] that could efficiently find a finite force or impulsive solution for planar multi-rigid-body problems; he also extended his algorithm to handle three-dimensional systems. But there is a subtle flaw in the handling of the forces at rolling contacts that can cause Baraff's algorithm to return contact forces and accelerations that do not satisfy the model at rolling contacts. Nonetheless, animations based on his algorithm are intuitively pleasing and incorporating his algorithm into VR systems would be likely to yield "believable" results.

Erdmann used a rigid-body model in his development of a manipulation planner for robotic systems that operated in the space of generalized coordinates referred to as configuration space (C-space). In C-space, the multi-body system was represented as a point, which was restricted to lie on a lower-dimensional "surface" if bodies were maintaining contact. Since Erdmann's planner operated in C-space, he studied the contact problem in that space. He developed the notion of the C-space friction cone and studied its properties [9]. This friction cone provides a geometric method to determine the possible motions of the rigid object in contact with immobile walls, given its initial velocity and an applied force and torque. Erdmann noticed that when the C-space friction cone corresponding to a single contact point dips below the C-space tangent plane, the model of the system may have zero, one, or more solutions. For the sliding rod problem, Erdmann's conditions are equivalent to L"otstedt's and ours.

Finally, we note that with the difficulties discussed above, one may wonder whether it would be beneficial to include a model of contact compliance. Abel, Kumar and Wang [29] showed that multibody contact models admit unique solutions when compliance is lumped at the contact points. This fact makes compliant models desirable in simulation environments, but the addition of compliance leads to stiff differential equations (assuming that the bodies are stiff as many every-day items typically are). However, there is currently no efficient and accurate way to choose the compliance parameters. In addition, for motion planning problems, the uncertainty in the compliance parameters is a serious problem, because the planner uses the model solutions to determine the response of the system to its proposed actions. If the solution of the model is sensitive to the compliance values, then planned tasks may fail for no apparent reason when executed by the real system. If one attempts to combat this problem by expanding the configuration space to include unknown compliance parameters, the growth in the dimension of this space would make all but the simplest planning problems intractable.

2 The Multi-Rigid-Body Contact Model The objective of this section is to introduce the fundamental model for the multi-rigid-body contact problem. This model, and an ability to solve it, is fundamental to any dynamic simulator or automatic task planner for multi-rigid-body systems. During a typical integration step, the equations must be formulated and solved to determine the accelerations of the bodies. The accelerations are then used to update the velocities and positions of the bodies. After this, the bodies would be tested for collisions. Upon the detection of a collision, an impulse model would be applied to

6

prevent interpenetration. Then the process repeats. For the purposes of this paper, we will assume that an impulse model has already been applied if necessary. We will study only the formulation of the equations of motion and the applicable solution techniques for a given time step. The time instant for which the equations of motion are formulated will be referred to as the current time.

The assumptions that will apply for the remainder of this paper are:

1. The bodies are rigid. 2. The normal direction at each contact is well-defined. 3. Dry friction exists at each contact point. 4. Each manipulator joint has one degree of freedom. 5. The manipulator has no closed loops formed by the links and joints. Closed loops involving

(unilateral) contacts are allowed.

6. All links are connected (at least indirectly) to a grounded link as indicated by the hatching

in Figure 1.

The fourth assumption is not limiting, because joints with more than one degree of freedom can be modeled as a set of one degree-of-freedom joints. The reason for the fifth assumption is to avoid having to consider load distribution in the manipulator joints. This assumption could be removed by applying known techniques. Assumption 6 is applicable to typical Earth-bound systems, (e.g., manufacturing systems). Removing this assumption would fundamentally alter the formulation of the model. Notice that there are no restrictions on the shapes of the bodies.

2.1 The Newton-Euler equations While the derivation of the basic equations has appeared in the literature, we derive them briefly here in a manner that facilitates their formulation as complementarity problems later.

Let us begin by numbering the objects from 1 to nobj and manipulator links from 1 to nman. All grounded links are considered to be a single manipulator link. When two bodies i and k (with i 6= k) come in contact, we label the contact point as j and consider this point j as associated uniquely with the pair (i; k). Let nc denote the total number of contact points at the current time. Each contact point j of bodies i and k defines the origin of the contact frame Cj. Without loss of generality, define the contact normal, ^nj, to point inward with respect to body i, if i ! k, and outward otherwise. The other two axes of frame Cj; ^tj and ^oj , span the contact tangent plane and satisfy ^nj \Theta ^tj = ^oj. Further, let ci;j and ck;j denote, respectively, the contact force acting on body i and k through contact j and expressed in frame Cj. For ` = i; k, we may write

c`;j j 2664

(c`;j)n

(c`;j)t (c`;j)o

377 5 (1)

7

where (c`;j)n is the component of the contact force c`;j in the direction of the contact normal, ^nj, and (c`;j)t and (c`;j)o are the force components in the contact tangent plane spanned by ^tj and ^oj. By Newton's third law, we have ck;j = \Gamma ci;j. Next, define W `;j to be the (6 \Theta  3) wrench matrix that transforms contact force c`;j into the equivalent wrench (generalized force) in the frame B`. The wrench matrices contain all geometric information relevant to contact j (definitions can be found in [27]).

Assuming that object i is acted upon by an external generalized force vector gobj;i of length 6, and a set of contact forces, the Newton-Euler equation governing the motion of object i isX

j2Bi

W i;j ci;j + gobj;i + hobj;i = M obj;i"qi (2)

where Bi is the index set of contact points that involve body i; hobj;i is the vector of velocity product terms, "qi is the 6-vector of linear and angular accelerations of the center of mass, and M obj;i is the (6 \Theta  6) positive definite and symmetric mass matrix of object i. This matrix is block diagonal with the upper left (3 \Theta  3) block given by the mass mobj;i times the (3 \Theta  3) identity matrix and the lower right block given by the (3 \Theta  3) inertia tensor of body i expressed in frame Bi.

The Newton-Euler equations of all the objects can be combined and cast in matrix form. For this purpose, it would be convenient to introduce the global wrench intensity vector c of length 3nc which is composed of nc subvectors cfl;j , each of dimension 3, where j runs from 1 to nc and for each j that corresponds to contact of bodies i and k; fl is the smaller of i and k. Note that fl is uniquely determined by j. Concatenating the equations (2) for i = 1; : : : ; nobj, we obtain

W c + gobj + hobj = M obj"q (3) where gobj; hobj, and "q are vectors of length 6nobj formed by the vertical concatenation of the 6-vectors gobj;i, hobj;i, and "qi respectively; the mass matrix of the objects, M obj, is block diagonal with dimension (6nobj \Theta  6nobj) and is given by M obj j diagfM obj;1; : : : ; Mobj;nobjg; and W is the

(6nobj \Theta  3nc) global wrench matrix composed of (nobj \Theta  nc) blocks, each of which is of dimension 6 \Theta  3 and most of which are identically zero. Given that bodies i and k (with i ! k) are both passive objects and they touch at contact j, then the (i; j) and (k; j) (block) elements of W are nonzero. Block (i; j) is W i;j and block (k; j) is \Gamma W k;j, where the minus sign arises as a result of eliminating ck;j for k ? i by applying the relation ci;j = \Gamma ck;j. If body i or k is a link of the manipulator, then correspondingly, block (i; j) or (k; j) of W is zero.

Consistent with the composition (1) of each contact vector ci;j, we denote the three columns of W corresponding to the j-th contact as (wi;j)n, (wi;j)t, and (wi;j)o. Next, define the n; t, and o partitions of W to be the (6nobj \Theta  nc) matrices W n; W t, and W o, with their j-th columns equal to (wi;j)n, (wi;j)t, and (wi;j)o, respectively. Finally, let cn, ct, and co be the nc-vectors whose j-th components are (cfl;j)n; (cfl;j)t, and (cfl;j)o, respectively, where fl is as defined above. Then equation (3) can be rewritten as follows:

W ncn + W tct + W oco + gobj + hobj = M obj"q: (4)

The development of the dynamic equations of the manipulator parallels the development of equation (3). Let n` be the number of joints in the manipulator and let o/ be the n`-vector of

8

joint efforts 2. When a link of the manipulator, body i, is involved in contact j, the effect of the contact force, ci;j, on the joint effort vector o/ can be expressed as the product of JTi;j (the first three columns of the (n` \Theta  6) transposed manipulator Jacobian matrix) times the contact force, ci;j. The definition of this matrix can be found in most introductory text books on robotics (e.g., [7]).

Summing the effects of all the contacts on the manipulator yields the following dynamic equation of the manipulator:

o/ \Gamma  (JT c + gman + hman) = M man"`

where J T is the (n` \Theta  3nc) global (transposed) Jacobian matrix, M man is the (n` \Theta  n`) positive definite and symmetric inertia matrix, gman and hman, are the n`-vectors of joint efforts induced by external wrenches and velocity product wrenches, respectively, and "` is the n`-vector of joint accelerations. The matrix JT is formed from the individual transposed Jacobian matrices as follows. Assume that contact j is between bodies i and k with i ! k. If both bodies are links of the manipulator, then the j-th set of three columns in JT is given by JTi;j \Gamma  JTk;j. If body i or k is an object, then the corresponding individual transposed Jacobian matrix is removed from the definition.

Partitioning the wrench intensity vector, c, into its n, t, and o components yields:

o/ \Gamma  (JTn cn + JTt ct + J To co + gman + hman) = M man"` (5) where J Tn , JTt , and JTo all have dimension (n` \Theta  nc).

2.2 Kinematic constraints According to our convention, a contact point j uniquely determines the two bodies i and k that touch at this point, and vice versa. Throughout the rest of the paper, we shall refer only to the contact points and not to the bodies in contact; in particular, the pair (i; j) which refers to body i in contact with another body k ? i at contact point j is abbreviated as j (with i omitted). We will use the notation jff to mean contact point j along the direction ff 2 fn; t; og.

At the current time, there are nc = nR + nS contacts, where nR and nS are the number of rolling and sliding contacts, and R and S are the index sets corresponding to the rolling and sliding contacts, respectively. Let vj = [vjn vjt vjo]T and aj = [ajn ajt ajo]T , for j = 1; : : :; nc, be the relative linear velocity and acceleration vectors, expressed in frame Cj of bodies i and k at contact point j. Next, partition vn, vt, and vo and an, at, and ao into their n, t, and o component vectors, vn, vt, vo, an, at, and ao, each of length n`. The vectors vff and aff, for ff 2 fn; t; og, can be expressed in the following convenient form:

vff = W Tff .q \Gamma  Jff .`; ff 2 fn; t; og;

aff = W Tff "q \Gamma  Jff"` + .W Tff .q \Gamma  .Jff .`; ff 2 fn; t; og; (6) 2By "joint effort" we mean the force in the direction of motion of a prismatic joint or the moment in the direction of motion of a revolute joint.

9

By definition, the normal component of relative velocity at a contact must be zero and rolling is indicated by both the t and o components being zero also. To prevent interpenetration, the normal component of relative acceleration at must be nonnegative. If it is zero, the contact will be maintained. Otherwise, the contact will break. Specifically, we have: j 2 R if v2jt + v2jo = 0 and j 2 S otherwise. The following equations and inequalities apply at the current time:

vn = 0 (7) vjt = vjo = 0; for all j 2 R (8)

an * 0: (9)

Notice that equations (7) and (8) constrain the initial conditions of the problem at the current time, but inequality (9) is a constraint on the unknown accelerations of the bodies.

2.3 Friction constraints The remaining constraints enforce certain friction laws. Two will be considered: first, Coulomb's Law which requires each contact force to lie within a quadratic cone, and second, an approximation of Coulomb's Law in which the quadratic cone at each rolling contact is replaced by a four-sided pyramid.

2.3.1 Coulomb's friction law Coulomb's Law stipulates that the j-th contact force cj lies within or on the boundary of its corresponding friction cone represented as follows:

c2jt + c2jo ^ _2j c2jn; for j = 1; : : :; nc (10) where _j is the coefficient of friction at the j-th contact point. Since the contact forces must be nontensile, we have:

cjn * 0; for j = 1; : : :; nc: (11)

If a contact is rolling, the contact force may have any direction and magnitude, provided it lies in the cone defined above. However, in the case of a sliding contact, Coulomb's law eliminates the two degrees of freedom associated with the direction of the force. For each sliding contact, the force must lie on the boundary of the friction cone with its friction component directly opposing relative sliding velocity:

_j cjnvjt + cjtqv2jt + v2jo = 0

_jcjnvjo + cjoqv2jt + v2jo = 0

9?=? ; for all j 2 S: (12) Similarly, at each rolling contact we have:

_j cjnajt + cjtqa2jt + a2jo = 0 _j cjnajo + cjoqa2jt + a2jo = 0

9?= ?; for all j 2 R: (13)

10

Physically, these two equations implement the correct sliding constraints if sliding is imminent (i.e., when a2jt + a2jo 6= 0), and provide no constraint otherwise.

The two sets of sliding constraints differ in one significant way. The first pair is linear in the unknowns, cjn, cjt, and cjo, while the second pair is not, because both ajt and ajo are unknown. After the formulation of the nonlinear model based on the quadratic friction cones, this observation will motivate us to approximate the friction cones at the rolling contacts with friction pyramids. As will be seen, the approximation compromises Coulomb's law slightly, but allows the application of linear complementarity theory [6].

Finally, the model of the multi-rigid-body system must allow for four different types of contact transitions. These are:

(i) rolling ! rolling:

cjn * 0 ajn = ajt = ajo = 0; for all j 2 R;

(ii) rolling ! sliding:

cjn * 0 ajn = 0 a2jt + a2jo 6= 0; for all j 2 R; (14)

cjt = \Gamma  _jajtq

a2jt + a2jo

cjn cjo = \Gamma  _j ajoq

a2jt + a2jo

cjn; for all j 2 R; (15)

(iii) sliding ! sliding:

cjn * 0; ajn = 0; for all j 2 S;

(iv) rolling or sliding ! breaking:

ajn * 0 cjn = 0 () cjt = cjo = 0); for j = 1; : : :; nc:

The remaining transitions from sliding to rolling, and from no contact to contact, are indicated when the tangential velocity at a contact point goes to zero, and when the geometric models of two bodies overlap.

Three-Dimensional Multi-Rigid-Body Contact Problem with Coulomb Friction Law: Given q, .q, gobj, hobj , M obj , `, .`, o/ , gman, hman, M man, and the W and J matrices, with .q and

.` satisfying equations (7) and (8); determine "q, "`, cn, ct, co, an, at, ao satisfying equations (4-6),

(9-13), and the complementarity condition:

(cn)T an = 0: (16)

Note that the four contact transitions (i)-(iv) are natural consequences of the defining conditions of the problem.

11

2.4 The friction pyramid law One reason for approximating the friction cone by a pyramid, is that the nonlinearities defining the quadratic cone constraints (13) make it difficult to develop robust, efficient algorithms to compute solutions (when they exist). Figure 2 shows a friction cone and its approximation as a four-sided friction pyramid. The interior and boundary of the pyramid can be expressed as a system of linear inequalities:

_j cjn \Gamma  cjt * 0 _j cjn + cjt * 0

_j cjn \Gamma  cjo * 0 _j cjn + cjo * 0 ) for all j 2 R: (17)

n^j

t^j o^j

Figure 2: Friction Cone and Pyramid When using the friction pyramid, we replace the rolling contact conditions (13) by the following conditions:

_jcjnajt + cjtjajtj = 0

_j cjnajo + cjtjajoj = 0 ) for all j 2 R: (18) Accordingly, the rolling-to-sliding contact interaction constraints (14) and (15) now become

(ii) rolling ! sliding:

cjn * 0 ajn = 0 ajt 6= 0 OR ajo 6= 0; for all j 2 R: (19)

cjt = ( \Gamma 

_jajt

jajtj cjn if ajt 6= 0

2 [\Gamma _jcjn; _jcjn] otherwise

cjo = ( \Gamma 

_jajo

jajoj cjn if ajo 6= 0

2 [\Gamma _jcjn; _jcjn] otherwise

9??????= ??????; for all j 2 R: (20)

Three-Dimensional Multi-Rigid-Body Contact Problem with Friction Pyramid Law: Given q, .q, gobj, hobj , M obj , `, .`, o/ , gman, hman, M man, and the W and J matrices, with .q and

12

.` satisfying equations (7) and (8); determine "q, "`, cn, ct, co, an, at, ao satisfying equations (4-6), (9), (12), and (16-18).

There is one small problem that arises with the friction pyramid law when a rolling contact converts to sliding. If the model predicts that a given rolling contact begins to slide (i.e., the solution of the model satisfies equation (19)), then according to equation (18), the friction force will generally not exactly oppose the direction of sliding. Figure 3 illustrates this problem. Practically speaking, the friction force always points into one of the four corners of the square shown even though the relative acceleration of the contact in the tangent plane does not.

^t o^

Friction

AccelerationRelative

Force

j

j

Figure 3: Projection of Friction Pyramid onto Tangent Plane On the other hand, the friction pyramid law has the important benefit of leading to a computationally tractable model; Coulomb's law (with the quadratic cone) results in a computationally more difficult problem. In addition, the multi-rigid-body contact model can be extended to allow any number of facets in the friction pyramid while retaining the form of an LCP. With this enhancement, our model can be used solve the model with Coulomb's Law by iteratively solving the LCP formulation with increasing numbers of facets in the pyramids until all contact forces at the rolling contacts lie within the quadratic cones they circumscribe.

2.5 Variations of the model There are two specific variations of the multi-rigid-body contact model that are useful to consider and can be studied using the same techniques as will be used below. First, the input, o/ , and output," ` can be interchanged. That is, one can assume that the joint effort vector, o/ , is unknown, while the joint acceleration vector, "` is known. More generally, let o/ ff be the subset of the elements of o/ that are known and o/ _ff denote the unknown elements. The corresponding subsets of elements of" `, "`ff and "` _ff, must be, respectively, unknown and known.

Second, the contact friction model can be extended to what is referred to as a "soft finger" model [25]. This model assumes that in addition to three components of force, a contact can also transmit a moment in the direction of the contact normal. In this extension, the wrench intensity vector, cj

13

is augmented with the signed magnitude of the normal moment, that is cj = [cjn; cjt; cjo; cjo/], and the magnitude of the friction moment about the contact normal, cjo/ is assumed to be bounded by \Gamma _o/ cjn and _o/ cjn when there is no relative rotation at the contact. When there is relative rotation, we have cjo/ = \Gamma sgn( .OE)_o/ cjn where .OE is the normal component of the relative angular velocity and "sgn" is the sign function of a scalar.

There are two interesting variations that change the form of the basic model. The first is the quasistatic, multi-rigid-body contact problem which is formed by setting the terms M obj"q and M man"` to zero. The quasistatic assumption fundamentally alters the model in ways discussed in [22]. The other variation arises when relaxing the rigid body assumption. This case, the unknown contact forces can be replaced by functions of the configuration and velocity of the system, which relate forces to deformations.

3 Complementarity Formulations The nonlinear complementarity problem (NCP) is to find an unknown vector z, satisfying the following complementarity system:

z * 0 f (z) * 0 zT f (z) = 0 (21) where f is a function mapping the n-dimensional Euclidean space into itself. If f (z) = r + M z, where r is a given n-vector and M is a given n \Theta  n matrix, then the NCP is known as a linear complementarity problem, which we denote as LCP (r; M).

While it is clear from the development in the previous section and from previous work on this problem, that the constraints on an and cn are complementary, it is not clear whether or not it is possible to write the complete set of equations and inequalities describing the multi-rigid-body contact problem in the standard form given above. In this section, we will show that it is possible to formulate the multi-rigid-body contact problem as an LCP if one is willing to approximate the quadratic Coulomb friction cone by a four-sided friction pyramid. When the quadratic cones are used, the problem does not quite fit the form of a standard NCP. However, it can be cast as a type of generalized complementarity problem that has been studied previously.

The development of both formulations begins by eliminating the accelerations, "q and "`, as follows:

"q = M \Gamma 1obj hW ncn + W tct + W oco + gobj + hobji (22) "` = M \Gamma 1man ho/ \Gamma  JTn cn \Gamma  JTt ct \Gamma  J To co \Gamma  gman \Gamma  hmani : (23)

Note that because M obj and M man are positive definite and symmetric, their inverses exist.

Substituting these expressions into equation (6), defining an, and the analogous ones defining at and ao yields: 266

4

an

at ao

3775

= A 2664

cn

ct co

3775

+ 2664

bn

bt bo

3775

(24)

14

where

A = 2664

Ann Ant Ano

Atn Att Ato Aon Aot Aoo

377 5 j J

T MJ

with the system inertia and constraint Jacobian matrices, M and J , defined as follows:

M j 24

M \Gamma 1obj 0

0 M \Gamma 1man

35

J j 24

W n W t W o

JTn JTt JTo

35

(25)

and 266

4

bn

bt bo

3775

j .J T 24 .

q

\Gamma  .`

35

+ J T M "

gobj + hobj

gman + hman \Gamma  o/ # : (26)

Notice that the matrix A is positive semidefinite and symmetric.

Next, we can solve the equations (12) for the friction force components, cjt and cjo, to eliminate those corresponding to the sliding contacts. Since there are no additional restrictions on the tangential accelerations, ajt and ajo, for j 2 S, the equations in (24) that define these components can be dropped from the model without affecting its solvability. Let the vectors aSn and aRn be the partitions of an corresponding to the sliding and rolling contacts respectively (i.e., the elements of aSn are ajn for j 2 S, and aRn is defined similarly). Next define aSt, aRt, aSo, aRo, cSn, cRn, cSt, cRt, cSo, and cRo analogously. Eliminating cSt and cSo, and removing the equations defining aSt and aSo result in: 266

66664

aSn aRn

aRt aRo

3777 7775 = ~A

2666 6664

cSn cRn

cRt cRo

3777 7775 +

2666 6664

bSn bRn

bRt bRo

3777 7775

where

~A =

2666 66664

( ~Ann)SS (Ann)SR (Ant)SR (Ano)SR ( ~Ann)RS (Ann)RR (Ant)RR (Ano)RR

( ~Atn)RS (Atn)RR (Att)RR (Ato)RR ( ~Aon)RS (Aon)RR (Aot)RR (Aoo)RR

3777 77775 : (27)

Here 266

66666 4

( ~Ann)SS ( ~Ann)RS

( ~Atn)RS ( ~Aon)RS

3777 77775 j

2666 6664

(Ann)SS (Ann)RS

(Atn)RS (Aon)RS

3777 7775 \Gamma 

2666 6664

(Ant)SS (Ant)RS

(Att)RS (Aot)RS

3777 7775 V St \Gamma 

2666 6664

(Ano)SS (Ano)RS

(Ato)RS (Aoo)RS

3777 7775 V So (28)

15

where V St and V So are, respectively, the diagonal matrices with diagonal entries given by

_j vjtq

v2jt + v2jo and

_j vjoq v2jt + v2jo

; for all j 2 S

and in general, for an N \Theta  N matrix M , if ff and fi are subsets of f1; : : :; N g; Mfffi denotes the submatrix of M consisting of rows and columns indexed by ff and fi respectively.

To manipulate the system of equations closer to the standard complementarity form, for each rolling contact j, we replace the components of the friction forces and tangential contact accelerations, cjt; cjo; ajt, and ajo, with the differences of their nonnegative and nonpositive parts as follows:

cjt j c+jt \Gamma  c\Gamma jt cjo j c+jo \Gamma  c\Gamma jo

ajt j a+jt \Gamma  a\Gamma jt ajo j a+jo \Gamma  a\Gamma jo

9=;

for all j 2 R

where a+jt; a\Gamma jt; a+jo; a\Gamma jo; c+jt; c\Gamma jt; c+jo, and c\Gamma jo are nonnegative. Substituting into equation (24) gives:

2666 6664

aSn aRn

a+Rt a+Ro

3777 7775 = B

2666 66666 666666

66666 64

cSn cRn

c+Rt c\Gamma Rt c+Ro c\Gamma Ro a\Gamma Rt a\Gamma Ro

3777 77777 777777

77777 75

+

2666 6664

bSn bRn

bRt bRo

3777 7775 (29)

where B is the matrix given by26

66666 4

( ~Ann)SS (Ann)SR (Ant)SR \Gamma (Ant)SR (Ano)SR \Gamma (Ano)SR 0 0 ( ~Ann)RS (Ann)RR (Ant)RR \Gamma (Ant)RR (Ano)RR \Gamma (Ano)RR 0 0

( ~Atn)RS (Atn)RR (Att)RR \Gamma (Att)RR (Ato)RR \Gamma (Ato)RR \Gamma I 0 ( ~Aon)RS (Aon)RR (Aot)RR \Gamma (Aot)RR (Aoo)RR \Gamma (Aoo)RR 0 \Gamma I

3777 7775 :

3.1 Formulation with Coulomb's friction law A solution to equation (29) satisfies equations (4-6) and (12), but not necessarily the other relevant constraints listed in the first problem statement, namely, equations (9-11), (13), and (16). These constraints can be expressed in complementarity form as follows:

(an)T cn = 0 (30)

h

_2j c2jn \Gamma  (c+jt \Gamma  c\Gamma jt)2 \Gamma  (c+jo \Gamma  c\Gamma jo)2i h(a+jt \Gamma  a\Gamma jt)2 + (a+jo \Gamma  a\Gamma jo)2i = 0; for all j 2 R (31)

16

i

a+jt \Gamma  a\Gamma jtj ic+jo \Gamma  c\Gamma joj \Gamma  ia+jo \Gamma  a\Gamma joj ic+jt \Gamma  c\Gamma jtj = 0; for all j 2 R (32)

(a+Rt)T c+Rt = (a\Gamma Rt)T c\Gamma Rt = (a+Ro)T c+Ro = (a\Gamma Ro)T c\Gamma Ro = 0 (33)

an; a+Rt; a\Gamma Rt; a+Ro; a\Gamma Ro; cn; c+Rt; c\Gamma Rt; c+Ro; c\Gamma Ro * 0 (34) where the vectors superscripted with + or \Gamma  contain the positive or negative parts of the vectors without the superscripts.

The first constraint, equation (30) applies to all contact points. It enforces the fact that if a contact force is compressive, then the relative acceleration of the bodies in the normal direction at that contact must be zero (i.e., the contact is maintained). Similarly, if the normal component of the relative acceleration at a contact is positive (i.e., the contact is breaking), then the normal component of the contact force must be zero (and since the force must lie within or on the friction cone, the entire contact force must be zero). The remaining constraints apply only to the rolling contacts; those with indices in the set R. Constraint (31) requires that at each rolling contact, either the tangential acceleration components, ajt and ajo are zero or the contact force lies on the boundary of the friction cone. Constraint (32) ensures that the direction of the friction force is parallel or anti-parallel to the sliding direction. Constraint (33) guarantees that the components of the friction forces and the relative sliding directions at the (initially) rolling contacts, have opposite signs. This constraint, when combined with the last one and the previous two restrict the friction forces at the rolling contacts to directly oppose the direction of sliding (if sliding begins). Notice that if the tangential acceleration at a rolling contact is zero, indicating that the contact continues to roll, then equations (31-33) are trivially satisfied. In this case, the only constraint on the contact forces at the rolling contacts is that they remain within their respective friction cones defined by inequality (10).

To summarize the above development, the three-dimensional multi-rigid-body contact problem with Coulomb friction cones is equivalent to equations (29-34). Notice that equations (31) and (32) prevent the formulation from being a standard NCP.

3.2 Formulation with friction pyramid law The key to the formulation of the three-dimensional multi-rigid-body contact problem with the friction pyramid law is the introduction of four nonnegative slack variables, s+jt; s\Gamma jt; s+jo, and s\Gamma jo, at each rolling contact:

s+jt j _j cjn + cjt s\Gamma jt j _jcjn \Gamma  cjt s+jo j _j cjn + cjo s\Gamma jo j _jcjn \Gamma  cjo

9=;

for all j 2 R: (35)

As shown in Figure 4, s+jt represents the amount the friction force component in the ^tj-direction can be decreased before the contact force reaches the boundary of the pyramid. At this point, sliding in the ^tj direction is incipient. Similarly, sliding in the \Gamma ^tj direction is imminent when s\Gamma jt goes to

zero. The other two slack variables, s+jo and s\Gamma jo have analogous interpretations.

17

jt^

jn^ s ForceContact

-sjt +jt

Figure 4: Slack Variables of the Friction Pyramid It is clear that the slack variables have a natural complementarity with their corresponding tangential contact accelerations as follows:

s+jta+jt = s\Gamma jta\Gamma jt = s+joa+jo = s\Gamma joa\Gamma jo = 0: (36) The next step is to substitute equations (35) into equation (29), displacing the elements of the variables c+Rt, c\Gamma Rt c+Ro, and c\Gamma Ro, to yield:2

666666 66666 64

aSn aRn

a+Rt a+Ro

s\Gamma Rt s\Gamma Ro

3777 777777 7775

= M

2666 666666 6664

cSn cRn

s+Rt s+Ro a\Gamma Rt a\Gamma Ro

3777 777777 7775

+

2666 666666 6664

bSn bRn

bRt bRo

0 0

3777 777777 7775

(37)

where

M j

2666 66666 66664

(M nn)SS (M nn)SR (Ant)SR (Ano)SR 0 0 (M nn)RS (M nn)RR (Ant)RR (Ano)RR 0 0

(Mtn)RS (M tn)RR (Att)RR (Ato)RR I 0 (M on)RS (M on)RR (Aot)RR (Aoo)RR 0 I

0 2U R \Gamma I 0 0 0 0 2U R 0 \Gamma I 0 0

3777 77777 77775

18

with U R being the diagonal matrix with diagonal entries _j, for j 2 R, and2

666666 4

(M nn)SS (M nn)SR (M nn)RS (Mnn)RR

(M tn)RS (M tn)RR (M on)RS (M on)RR

3777 7775 j

2666 6664

(Ann)SS (Ann)SR (Ann)RS (Ann)RR

(Atn)RS (Atn)RR (Aon)RS (Aon)RR

3777 7775 \Gamma 

2666 6664

(Ant)SS (Ant)SR (Ant)RS (Ant)RR

(Att)RS (Att)RR (Aot)RS (Aot)RR

3777 7775 " V St 0

0 U R # \Gamma 

2666 6664

(Ano)SS (Ano)SR (Ano)RS (Ano)RR

(Ato)RS (Ato)RR (Aoo)RS (Aoo)RR

3777 7775 " V So 0

0 U R # :

Solutions to equation (37) satisfy equations (4-6), and (12), but not necessarily the other constraints of the model, (9) and (16-18). These constraints have been shown to be equivalent to the nonnegativity of the variables

an; a+Rt; a\Gamma Rt; a+Ro; a\Gamma Ro; cn; s+Rt; s\Gamma Rt; s+Ro; s\Gamma Ro * 0 (38) and the following complementarity relations

(an)T cn = (a+Rt)T (c+Rt) = (a\Gamma Rt)T (c\Gamma Rt) = (a+Ro)T (c+Ro) = (a\Gamma Ro)T (c\Gamma Ro) = 0: (39)

4 Existence and Uniqueness Results In this section, we will present our solution existence and uniqueness results for the contact problems and relate them to previously known results of L"otstedt [14] and Baraff [1].

4.1 Problems without friction In previous work it was established that the frictionless three-dimensional, multi-rigid-body contact problem has a unique solution for the accelerations of the bodies, "q and "`, but that contact force magnitudes, cn, are generally not unique [15]. For example, consider a frictionless right rectangular block oriented so that one of its pairs of parallel faces is perpendicular to the direction of gravitational acceleration. Assume that the other four faces are each contacted near their four corners by fixed active bodies. Clearly, the acceleration of the block will be the gravitational

19

acceleration regardless of the magnitudes of the contact forces. In other words, the vector of relevant contact accelerations is uniquely determined, even though the vector of relevant contact forces is not.

Our two formulations reduce to L"otstedt's LCP formulation when the coefficients of friction are all zero. We have no new results relating to the frictionless case.

4.2 Problems with friction Obtaining existence results for problems with friction is considerably more difficult than the frictionless problems. Indeed, we need to rely on the theory of quasi-variational inequalities.

The all-rolling problem First, consider the restricted problem in which all contacts are rolling at the current time (Note that this restriction does not preclude solutions in which some of the rolling contacts break or convert to sliding.). We refer to this as the "all-rolling problem". Mathematically, it corresponds to the case where S = ;.

Using the theory of quadratic programming, L"otstedt showed that the all-rolling problem in the plane has a solution [14] if the following Jacobian matrix has full column rank:2

4 W n W t

JTn J Tt

35

: (40)

Baraff later showed in [1], that all planar problems with only one contact have a solution (regardless of the rank of this matrix), and he suggested that his proof extends trivially to three-dimensional problems with one contact. However, it is worth noting that Jacobian matrix (40) and its threedimensional counterpart, always have full column rank when there is only one contact. In [21], we used the theory of quasi-variational inequalities to prove that a solution to the three-dimensional allrolling problem must exist for both the cone and pyramid friction models under a mild assumption. This result includes L"otstedt's and Baraff's results as special cases. Our result is repeated below for the reader's convenience.

Suppose S = ;. Let F be the subset of the Euclidean space R3nc consisting of triples (cn; ct; co) satisfying the friction cone constraints (10) and (11): cjn * 0 and _2j c2jn \Gamma  c2jt \Gamma  c2jo * 0, for all j = 1; : : : ; nc. Further, let N be the null space of the system Jacobian matrix (defined in equation (25)). Let FR j F " N .

A solution to the three-dimensional all-rolling problem exists under the following condition [21]:2

4 .q

\Gamma  .`

35T

.J 2664

cn

ct co

377 5 * 0; for all (cn; ct; co) 2 FR: (41)

20

If the quadratic friction cone is replaced by the friction pyramid defined earlier, then by modifying the definition of the cone F accordingly, the satisfaction of the same condition (41) also implies the existence of a solution for the friction pyramid model.

Notice that the condition (41) is trivially satisfied when the columns of the system Jacobian matrix are linearly independent, since then its null space becomes the origin. Thus, our existence condition (when specialized to the planar case) is less restrictive than L"otstedt's and clearly subsumes both his and Baraff's existence results.

A situation in which the Jacobian matrix, J , could be column rank deficient and yet a solution is guaranteed, is when the constant vector, b, defined in equation (26), lies in the column space of J T . This follows because the second summand in the right-hand side of (26) always lies in the column space of J T , hence (41) must hold in this case. Note that in the special case of assembly stability testing, the body velocities (i.e., .q and .`) are zero, so the problem will always have a solution.

The sliding-rolling problem

L"otstedt [14] showed that the planar problem with both rolling and sliding contacts has a unique solution if every matrix in a family of matrices (closely related to ~A defined in equation (27)) belongs to the class of P-matrices; see Appendix A for definition. In effect, this result says that the matrix of every LCP that could arise by fixing the rolling contact force directions anywhere in their respective cones (and substituting to eliminate the tangential components of the contacts forces and accelerations), must be P.

The following result is an extension of L"otstedt's existence and uniqueness result to the threedimensional problem. It also extends our previous all-rolling existence result for three-dimensional systems [21, Theorem 4] to allow sliding at some contacts. The assertion about the uniqueness of solution is new even for the all-rolling problem.

Theorem 1 Suppose that the system Jacobian matrix, J (defined in equation (25)), has full column rank.

(a) There exists a positive scalar __ such that if _j 2 [0; __] for all j 2 S, the three-dimensional

multi-rigid-body contact problem with Coulomb friction cone law has a solution. The same is also true for the problem under the friction pyramid law.

(b) There exists a positive scalar __ such that if _j 2 [0; __] for all j = 1; : : :; nc, the threedimensional multi-rigid-body contact problem with Coulomb friction cone law has a unique solution. The same is also true for the problem under the friction pyramid law.

In the planar case, this theorem, while stated somewhat differently, becomes exactly L"otstedt's result. Our result says that a solution exists if the coefficients of friction are small enough. As we shall see from the proof in Appendix A, the full column rank and small friction coefficient assumptions imply that certain matrices are P. The latter P-property is the key to the validity of Theorem 1. L"otstedt has used the P-matrix assumption directly; we have phrased this assumption

21

in terms of the friction coefficients and the property of the system Jacobian matrix. Thus, both we and L"otstedt are making essentially the same assumptions, but from slightly different points of view.

One difference between the assumptions in (a) and (b) is that in (a) only those friction coefficients at the sliding contacts are assumed small; whereas in (b), all friction coefficients are assumed small. We should point out while it is in general difficult to compute the largest upper bound for __, we will give in Proposition 3 a conservative estimate for such a bound; such an estimate is expected to be much smaller than one would expect to encounter in real systems.

The detailed proof of Theorem 1 is given in Appendix A. The existence proof uses an argument that is similar to the all-rolling case [21]; the uniqueness proof is based on an extension of the argument for the planar case that L"otstedt used in [14].

In general, solutions to the two contact models could be quite distinct. The next result gives a simple condition for a solution of the friction pyramid model to be a solution of the cone model. The extent to which this observation can be put to use will be investigated separately. A proof of the result is also given in Appendix A.

Proposition 1 If a solution to the friction pyramid model satisfies the cone constraints (10), then it must be a solution of the friction cone model.

At first glance, this result seems trivial. Nevertheless a closer look at the two models suggests that they differ not only in the quadratic (10) versus polyhedral (17) friction constraints, but the resulting contact conditions, (13) versus (18), are also different. Thus a careful proof for Proposition 1 is needed; this proof turns out to be quite simple and is motivated by an equivalent formulation of the two models as quasi variational inequalities; see Appendix A. The validity of Proposition 1 requires none of the assumptions in Theorem 1.

5 Solution Methods for Friction Pyramid Model The formulation of the friction pyramid model as the LCP defined by (37)-(39) permits the numerical solution of this model by a host of algorithms as described in [6]. The main concern with the application of these algorithms is whether their convergence criteria are satisfied by the data of this particular LCP. In our study, we have focused on two algorithms: one classical and the other contemporary. The former is Lemke's almost complementarity pivot algorithm; the latter is a feasible interior point algorithm.

We refer the reader to [6] for a comprehensive treatment of Lemke's algorithm and its convergence theory. The following is the main convergence result of this algorithm applied to the LCP in question. A proof of this result is given in Appendix A.

Theorem 2 Suppose that the system Jacobian matrix, J (defined in equation (25)), has full column rank. There exists a positive scalar __ such that if _j 2 [0; __] for all j = 1; : : :; nc, Lemke's algorithm, under a standard nondegeneracy assumption, will in a finite number of pivots successfully compute a solution to the LCP (37)-(39).

22

The family of interior point methods is a recent entry into the field of the LCP. Fueled by their great success for solving linear programs, these methods have received much attention in the mathematical programming literature. For this reason, we have chosen a feasible interior point (FIP) method as a candidate algorithm for solving the LCP arising from the multi-rigid-body problem with friction pyramids. As we shall see in the next section, the interior point method provides a viable alternative to Lemke's method.

Since the engineering community is probably not well acquainted with the interior point method, we give a detailed explanation of the method in Appendix B. Moreover, due to the fact that existing results in the interior point literature are not directly applicable to the particular LCP in question (because of its failing an essential assumption in such results), we also give a detailed convergence analysis of our FIP method. In what follows, we state the main convergence result of this method for solving the friction pyramid model. Proof of this result is given in Appendix B.

Theorem 3 Suppose that the system Jacobian matrix, J (defined in equation (25)), has full column rank. There exists a positive scalar __ such that if _j 2 (0; __] for all j = 1; : : : ; nc, the feasible interior point method, as described in Appendix B and applied to the LCP defined by (37-39), generates a sequence of iterates i

a*n; c*n; a*;\Sigma Rt ; a*;\Sigma Ro ; s*;\Sigma Rt ; s*;\Sigma Ro j ; * = 0; 1; 2; : : : ;

each of which is strictly feasible to this LCP (i.e., equation (37 is satisfied and the inequalities in 38) are strictly satisfied) and

lim*!1(a*n)T c*n = lim*!1(a*;\Sigma Rt )T s*;\Sigma Rt = lim*!1(a*;\Sigma Ro )T s*;\Sigma Ro = 0: Moreover, every accumulation point of the generated sequence solves the LCP; finally, this LCP must have a solution.

We point out that the positivity of the coefficients of friction at the rolling contacts is essential for the validity of the above theorem. This is different from Theorem 2 which permits zero coefficients of friction at the rolling contacts.

6 Numerical Results In this section, we report our computational experience gained through applying Lemke's algorithm and our FIP algorithm (detailed in Appendix B) to some multi-rigid contact problems with friction pyramid law. The algorithms were implemented in MATLAB3, while our data generation code was written in C. Both algorithms were run on the same 328 data sets. These were generated by randomly generating 35 (physically meaningful) problems of varying sizes with varying numbers of passive bodies, contact geometries, and initial conditions. The dimensions of the corresponding M matrices in the LCPs ranged from 2 to 170. Each problem was used to produce about 10 data sets by varying the coefficient of friction, which for convenience, was assumed to be equal at all contact points. While the termination criterion for Lemke's algorithm was standard, our FIP algorithm

3Lemke's algorithm was provided by Michael Ferris.

23

was terminated if the total complementary gap fell below 10\Gamma 6. We set upper limits of 100 and 300 iterations in Lemke's and the FIP algorithms, respectively.

Table 1 summarizes our numerical results. However, due to the large number of data sets, only representative results are included in the table. The column headings in the table are defined as follows: "size" is the dimension of the M matrix appearing in the LCP formulation, "# data sets" is the number of data sets of a given size that were attempted, "frac. solved" is the fraction of the data sets solved by either Lemke's or our FIP algorithm (as indicated by the column headings), "ave. # starts" is the average number of times Lemke's algorithm was started for the data sets which it eventually solved, "ave. # iter." is the average number of iterations required for a solved data set4.

Lemke's algorithm FIP algorithm size # data sets frac. solved ave. # restarts ave. # iter. frac. solved ave. # iter.

2 20 1.00 1.0 2.0 0.50 8.4 10 60 0.92 1.3 8.9 0.35 10.4 16 10 1.00 1.4 20.0 1.00 17.2 30 26 0.62 1.7 44.1 0.42 19.6 50 20 0.35 2.0 107.4 0.30 61.1 90 23 0.00 - - 0.35 30.4 120 10 0.10 1.0 ! 100 0.70 ! 300 150 6 0.00 - - 0.50 46.3 170 10 0.00 - - 0.00 -

Table 1: Summary of numerical results

While some of the interesting trends found in the results are evident in Table 1, not all of them are, so the remainder of this section is devoted to pointing them out. Overall, the two algorithms solved just under 70% (228 out of 328) of the data sets. However, since the data sets were generated randomly, we do not know a priori how many had solutions. It is conceivable that the unsolved data sets did not have solutions. We do not know for sure at this time. Also, since we have not yet developed a practical way to verify the assumptions of the various theorems in the preceding sections, we do not know if these assumptions are satisfied by the data sets whose solutions were found. Our suspicion is that they are not.

Of the 328 problems, 195 were solved by our FIP algorithm, while 136 were solved by Lemke's algorithm using a covering vector with all elements set to 1. On average, the numbers of iterations required by the algorithms for a given data set were comparable. As noted above, whenever Lemke's algorithm failed, it was restarted, up to three times, with a new covering vector chosen at random. With up to three restarts allowed, Lemke's algorithm solved about 10% more problems than the FIP algorithm, but the total number of iterations (including the failed runs) was then correspondingly higher. Even though Lemke's algorithm performed better overall (and could possibly do even better with more restarts and/or iterations allowed), the FIP method was superior when the system Jacobian did not have full column rank.

4The average number of iterations for Lemke's algorithm includes the iterations performed during failed runs on a given data set provided that it was eventually solved by Lemke's algorithm after a later restart.

24

Table 2 shows the results for a problem in which Lemke's algorithm found a solution on its first run only 3 out of 10 times, while our interior point method found solutions to 6 problems. All of these problems were solved by Lemke's algorithm on its second start. The columns headings are defined as follows: "nb" is the number of passive bodies, "nS " is the number of sliding contacts, "nR" is the number of rolling contacts, "_" is the coefficient of friction, "L-starts" is the number of times Lemke's algorithm was run (each time with a different covering vector), "L-iter" is the number of iterations required by Lemke's algorithm to solve the data set, and "FIP-iter" is the number of iterations required by the interior point method. Note that the appearance of "MAX" in either the "L-iter" or "FIP-iter" column means that the algorithms failed to meet our convergence criterion in 300 iterations. In the "L-starts" column, "4F" implies that Lemke's algorithm failed in four starts with different covering vectors. Finally, if "NFP" occurs in the "FIP-iter" column, then no strictly feasible point was found in the first phase of the interior point method, so the algorithm could not be initiated.

size nb nS nR _ L-starts L-iter FIP-iter

26 4 6 4 0.1 1 20 55 26 4 6 4 0.2 1 14 24 26 4 6 4 0.3 2 18 25 26 4 6 4 0.4 2 14 MAX 26 4 6 4 0.5 1 14 MAX 26 4 6 4 0.6 2 14 MAX 26 4 6 4 0.7 2 16 MAX 26 4 6 4 0.8 2 18 21 26 4 6 4 0.9 2 14 25 26 4 6 4 1.0 2 16 20

Table 2: Problem (with 10 data sets) for which Lemke's algorithm

(with up to 1 restart) out-performed the FIP algorithm.

A clear trend developed indicating that the algorithms' success was highly dependent on the geometry of the problems and to a lesser extent, on the value of the coefficient of friction used to generate the particular data sets. This was indicated by failures coming in bunches. That is, for a given problem, one or both algorithms typically failed for all values of the coefficient of friction above some bound. Table 3 shows the results for two problems (but 20 data sets), which exhibit this trend. Notice that in the first problem (10 data sets), Lemke's algorithm solved all ten instances quickly, while the interior point method could not find strictly feasible points. In the second problem, both algorithms worked well until the coefficient of friction reached about 0.5.

25

size nb nS nR _ L-starts L-iter FIP-iter

2 1 2 0 0.1 1 2 NFP 2 1 2 0 0.2 1 2 NFP 2 1 2 0 0.3 1 2 NFP 2 1 2 0 0.4 1 2 NFP 2 1 2 0 0.5 1 2 NFP 2 1 2 0 0.6 1 2 NFP 2 1 2 0 0.7 1 2 NFP 2 1 2 0 0.8 1 2 NFP 2 1 2 0 0.9 1 2 NFP 2 1 2 0 1.0 1 2 NFP 30 5 5 5 0.1 2 31 40 30 5 5 5 0.2 1 25 19 30 5 5 5 0.3 2 23 40 30 5 5 5 0.4 1 27 16 30 5 5 5 0.5 1 24 18 30 5 5 5 0.6 1 29 MAX 30 5 5 5 0.7 4F - MAX 30 5 5 5 0.8 4F - NFP 30 5 5 5 0.9 4F - NFP 30 5 5 5 1.0 4F - NFP

Table 3: Problems showing the dependence of algorithm performance

on problem geometry and the coefficient of friction.

As stated above, while both algorithms tended to have difficulty with the same data sets (particularly larger ones), in some cases, one algorithm performed well while the other one did not. The problem of size 2 in Table 3 shows Lemke's method out-performing the interior point method, while the problem in Table 4 shows the reverse.

size nb nS nR _ L-starts L-iter FIP-iter

55 1 5 10 0.1 2 51 17 55 1 5 10 0.2 2 55 16 55 1 5 10 0.3 2 41 MAX 55 1 5 10 0.4 2 57 15 55 1 5 10 0.5 2 84 23 55 1 5 10 0.6 4F - 43 55 1 5 10 0.7 4F - 17 55 1 5 10 0.8 4F - 15 55 1 5 10 0.9 4F - 14 55 1 5 10 1.0 4F - 40

Table 4: Problem showing the FIP algorithm out-performing Lemke's algorithm.

26

7 Summary and Conclusions In this paper, we have introduced various mathematical formulations for the three-dimensional multi-rigid-body contact problem with Coulomb friction. Existence and uniqueness of a solution to two models of the problem are established under a full column rank assumption on the system Jacobian matrix and a smallness assumption on the coefficients of friction. Two algorithms have been established to compute a solution to the friction pyramid model under the same assumptions. Although the required assumptions for the theoretical results are difficult to verify, the numerical results we have obtained suggest that the algorithms are fairly effective for solving some randomly generated three-dimensional multi-rigid-body contact problems under the friction pyramid law.

Obviously our study is imperfect. Theoretically, the results are not strong enough to handle the case of arbitrary friction coefficients and/or column rank deficiency of the system Jacobian matrix; numerically, we do not know if the two algorithms are capable of solving all friction pyramid problems which possess solutions. In spite of these deficiencies, we believe that we have made an important contribution toward the understanding of the three-dimensional multi-rigid-body contact problem with Coulomb friction. Theoretically, we have provided the most comprehensive results for this problem.

There remains much more to be accomplished. Trying to improve the theoretical results is clearly something worth looking into. The search for alternative algorithms for solving the friction pyramid model would be equally useful. In this regard, the NE/SQP method described in [20] offers a possibility; indeed, this method would also be a candidate for solving the cone model. This method and other approaches for dealing with the latter model are presently under investigation.

Finally, there is practical side to this study. Our motivation was the hope that our results would be useful in the design and development of two classes of systems: simulation systems and planning systems. Both types of systems would be extremely useful: the former helping us to better understand existing mechanical systems, and the latter enabling the automation of a plethora of contact tasks in both hazardous and safe environments.

The primary impediment to the achievement of our practical objectives is the fact that the multi-rigid-body models discussed here do not always have unique solutions, and the sufficient conditions that we have developed are conservative. What is needed are algorithms to always and efficiently compute a solution if one exists, or determine that either no solution or multiple solutions exist. With such algorithms available, one could develop intelligent strategies to deal with the nonuniqueness of the model in both simulation and planning applications. We are currently studying these and related issues.

Acknowledgements The authors would like to express their appreciation to Professor Michael Ferris at the University of Wisconsin at Madison for supplying an implementation of Lemke's algorithm written in MATLAB.

27

References

[1] D. Baraff. Issues in Computing Contact Forces for Non-Penetrating Rigid Bodies. Algorithmica,

10:292-352, 1993.

[2] D. Baraff. Fast contact force computation for nonpenetrating rigid bodies. In Proceedings,

SIGGRAPH '94, pages 23-34, July 1994.

[3] D. L. Brock. Enhancing the dexterity of a robot hand using controlled slip. Master's thesis,

MIT Department of Mechanical Engineering, 1987.

[4] R. C. Brost. Analysis and Planning of Planar Manipulation Tasks. PhD thesis, Carnegie

Mellon University School of Computer Science, January 1991.

[5] D. Chan and J.S. Pang. The generalized quasi-variational inequality problem. Mathematics

of Operations Research 7 (1982) 211-222.

[6] R. W. Cottle, J.-S. Pang, and R. E. Stone. The Linear Complementarity Problem. Academic

Press, 1992.

[7] J.J. Craig. Introduction to Robotics: Mechanics and Control. Addison-Wesley, second edition,

1989.

[8] M. A. Erdmann. On motion planning with uncertainty. Master's thesis, MIT Department of

Electrical Engineering and Computer Science, August 1984.

[9] M. A. Erdmann. On a representation of friction in configuration space. International Journal

of Robotics Research, 13(3):240-271, June 1994.

[10] E. J. Haug, S. C Wu, and S. M. Yang. Dynamic mechanical systems with Coulomb friction, stiction, impact and constraint addition-deletion--I: Theory. Mechanisms and Machine Theory, 21(5):407-416,1986.

[11] C. W. Kilmister and J. E. Reeve. Rational Mechanics. American Elsevier, New York, NY,

1966.

[12] M. Kojima, N. Megiddo, T. Noma, and A. Yoshise. A Unified Approach to Interior Point

Algorithms for Linear Complementarity Problem. Lecture Notes in Computer Science, 538, 1991.

[13] Z. Liu, S.J. Gortler, and M.F. Cohen. Hierarchical spacetime control. In Proceedings, SIGGRAPH '94, pages 35-42, July 1994.

[14] P. L"otstedt. Coulomb friction in two-dimensional rigid-body systems. Zeitschrift f"ur Angewandte Mathematik und Mechanik, 61:605-615, 1981.

[15] P. L"otstedt. Mechanical systems of rigid bodies subject to unilateral constraints. SIAM Journal

of Applied Mathematics, 42(2):281-296, 1982.

[16] K. Lynch. Stable pushing: Mechanics, controllability, and planning. In Proceedings, First

Workshop on Algorithmic Foundations of Robotics. A.K. Peters, Boston, MA, 1994.

28

[17] M. T. Mason and Y. Wang. On the inconsistency of rigid-body frictional planar mechanics.

In Proceedings, IEEE International Conference on Robotics and Automation, pages 524-528, April 1988.

[18] P. Painleve. Sur les lois du frottemont de glissement. C. R. Acad. Sci., 121:112-115, 1895. [19] R. Palmer and J. Cremer. Simlab: Automatically creating physical systems simulators. In

Proc. Automated Modeling Workshop of the 1992 ASME Winter Annual Meeting, November 1992.

[20] J.S. Pang and S.A. Gabriel. NE/SQP: A robust algorithm for the nonlinear complementarity

problem. Mathematical Programming, Series A 60:295-338, 1993.

[21] J.-S. Pang and J. C. Trinkle. Complementarity formulations and existence of solutions of multirigid-body contact problems with coulomb friction. Mathematical Programming. submitted August 1994.

[22] J.-S. Pang, J. C. Trinkle, and G. Lo. A complementarity approach to a quasistatic multi-rigidbody contact problem. Computational Optimization and Applications, 1995, forthcoming.

[23] M. A. Peshkin and A. C. Sanderson. The motion of a pushed, sliding workpiece. IEEE Journal

of Robotics and Automation, 4(6):569-598, December 1988.

[24] V. T. Rajan, R. Burridge, and J. T. Schwartz. Dynamics of a rigid body in frictional contact

with rigid walls. In Proceedings, IEEE International Conference on Robotics and Automation, pages 671-677, March 1987.

[25] J. K. Salisbury. Kinematic and Force Analysis of Articulated Hands. PhD thesis, Stanford

University, Department of Mechanical Engineering, May 1982. Reprinted in, Robot Hands and the Mechanics of Manipulation, MIT Press, Cambridge, Massachusetts (1985).

[26] J. C. Trinkle and J. J. Hunter. A framework for planning dexterous manipulation. In Proceedings, IEEE International Conference on Robotics and Automation, pages 1245-1251, April 1991.

[27] J. C. Trinkle and R. P. Paul. Planning for dexterous manipulation with sliding contacts.

International Journal of Robotics Research, 9(3):24-48, June 1990.

[28] J. C. Trinkle, R. C. Ram, A. O. Farahat, and P. F. Stiller. Dexterous manipulation planning and execution of an enveloped slippery workpiece. In Proceedings, IEEE International Conference on Robotics and Automation, volume 2, pages 442-448, May 1993.

[29] Y.-T. Wang, V. Kumar, and J. Abel. Dynamics of rigid bodies undergoing multiple frictional

contacts. In Proceedings, IEEE International Conference on Robotics and Automation, pages 2764-2769, May 1992.

[30] T. Wang, R.D.C. Monteiro, and J.S. Pang, "An interior point potential reduction algorithm

for constrained equations", manuscript, Department of Mathematical Sciences, The Johns Hopkins University, Baltimore, April 1994.

29

APPENDICES A Mathematical Details and Proofs The goal of this section is to give a rigorous proof of Theorem 1, Proposition 1, and Theorem 2. We begin with Theorem 1. Since the approach to establish the existence claim in this theorem is the same as that in [21] for the all-rolling case, we will simply explain the approach applied to the sliding-rolling case and omit the proofs of the supporting results for part (a) of Theorem 1. We will give the details of the proof for part (b) of this theorem and that of the other two results. For Theorem 1, we will focus on the cone model only; the treatment of the pyramid model is similar and thus omitted.

A.1 The cone model revisited For ease of reference, we restate the essential equations that define the three-dimensional multirigid-body contact problem with Coulomb friction cone as follows:2

66666 64

aSn aRn

aRt aRo

3777 7775 = ~A

2666 66664

cSn cRn

cRt cRo

3777 77775 +

2666 6664

bSn bRn

bRt bRo

3777 7775

(an; cn) * 0 (an)T cn = 0

c2jt + c2jo ^ _2j c2jn _j cjnajt + cjtqa2jt + a2jo = 0 _jcjnajo + cjoqa2jt + a2jo = 0

9????= ????; for all j 2 R

(42)

where the matrix ~A is given by (27), the vector b is given by (26), and R is the index set of the (initially) rolling contacts. Note that in the above formulation, all the tangential force and acceleration variables corresponding to the sliding contacts, cSt; cSo; aSt; aSo, have disappeared, due to the elimination via Coulomb's law for sliding contacts (12) and the unrestricted nature of the sliding accelerations.

Properties of the matrix ~A play a key role in the existence of solution to the model (42). To motivate the assumption that we wish to impose on ~A, consider the frictionless case where all friction coefficients at the sliding contacts, i.e. _j : j 2 S, are zero. In this case, the matrix ~A is equal to 266

66664

(Ann)SS (Ann)SR (Ant)SR (Ano)SR (Ann)RS (Ann)RR (Ant)RR (Ano)RR

(Atn)RS (Atn)RR (Att)RR (Ato)RR (Aon)RS (Aon)RR (Aot)RR (Aoo)RR

3777 7775 ; (43)

30

which is a principal submatrix of A. As we have noted before, A, (and hence each of its principal submatrices), is symmetric positive semidefinite. Moreover, A is positive definite if the system Jacobian matrix, J (defined in equation (25)), has full column rank. Thus under this full rank condition, the matrix ~A is positive definite in the frictionless case.

In the general model where the friction coefficients are positive, the matrix ~A is obtained by modifying A by some correction matrices that involve V St and V So which are diagonal matrices with diagonal entries dependent on the sliding friction coefficients; see (27). By a simple continuity argument, it follows that if these coefficients are (positive and) sufficiently small and if the system Jacobian matrix, J , has full column rank, then ~A will remain positive definite although ~A is no longer symmetric (i.e., the symmetric part of ~A is positive definite in the standard sense). In essence, this is what we will assume about ~A in the frictional case.

A.2 Some known concepts and results We summarize below the key concepts and results about certain matrix classes and the linear complementarity problem that we shall need later. Details of this brief review can be found in the text [6].

Definition. A real square matrix M , not necessarily symmetric, is said to be a P0-matrix if all its principal submatrices have nonnegative determinants; M is called a P-matrix if all its principal submatrices have positive determinants. The matrix M is said to be copositive if xT M x * 0 for all vectors x * 0.

These classes of matrices are extremely important in the theory of the linear complementarity problem. Clearly a P-matrix must be P0; however, the classes of copositive matrices and P0-matrices are different. In general, (asymmetric) positive semidefinite matrices must be both copositive and P0; (asymmetric) positive definite matrices must be P. Moreover, principal submatrices of a matrix belonging to one of these classes are matrices of the same class. Clearly, P-matrices are invariant under small perturbations; that is to say, if M is a P-matrix, then so are all matrices N that are sufficiently close to M .

The following result is key to our analysis; its proof can be found in [6].

Theorem 4 If M is a n \Theta  n P-matrix, then for every n-vector r, the LCP (r; M ) has a unique solution. Moreover, as a function of r, this solution is piecewise linear, and hence Lipschitz continuous.

There is a solvability result of the LCP (r; M ) that pertains to the case where M is a P0- matrix. Subsequently, we shall discuss this result in the context of a feasible interior point method; see Theorem 6.

We will also need a nonlinear generalization of a positive definite matrix.

Definition. A mapping F : Rn ! Rn is said to be strongly monotone if there exists a constant

31

p ? 0 such that for all vectors x; y 2 RN , the following inequality is satisfied:

(x \Gamma  y)T (F (x) \Gamma  F (y)) * pkx \Gamma  yk2:

A.3 Formulation as a quasi-variational inequality Taking into account the definition of the matrix ~A (27) and the above discussion, we state a result that justifies the assumption that we will impose on this matrix.

Proposition 2 Suppose that the system Jacobian matrix, J (defined in equation (25)), has full column rank. There exists a positive scalar __ such that if _j 2 [0; __] for all j 2 S, the matrix ~A is positive definite, and hence P.

In theory, an estimate of the scalar __ can be obtained by computing the smallest eigenvalue of the matrix (43), which is positive definite under the full rank assumption of the theorem, and dividing this eigenvalue by the norm of the correction matrix in (28) to form ~A. Specifically, defining

~B j

2666 6664

(Ann)SS (Ann)SR (Ant)SR (Ano)SR (Ann)RS (Ann)RR (Ant)RR (Ano)RR (Atn)RS (Atn)RR (Att)RR (Ato)RR (Aon)RS (Aon)RR (Aot)RR (Aoo)RR

3777 7775 ; (44)

Et j

2666 6664

(Ant)SS V St 0 0 0 (Ant)RS V St 0 0 0 (Att)RS V St 0 0 0 (Aot)RSV St 0 0 0

3777 7775 ; and Eo j

2666 6664

(Ano)SS V So 0 0 0 (Ano)RS V So 0 0 0 (Ato)RS V So 0 0 0 (Aoo)RSV So 0 0 0

3777 7775 ; (45)

we see that ~

A = ~B \Gamma  Et \Gamma  Eo: (46)

Notice that ~B is a principal submatrix of A and that Et and Eo do not contain the rolling friction coefficients. The following result gives a very conservative bound on the sliding friction coefficients which will guarantee the positive definiteness of ~A.

Proposition 3 Suppose that the matrix ~B is positive definite with smallest eigenvalue aemin ? 0. Let ffit and ffio be the spectral norm of the two matrices2

66666 64

(Ant)SS (Ant)RS (Att)RS (Aot)RS

3777 7775 and

2666 6664

(Ano)SS (Ano)RS (Ato)RS (Aoo)RS

3777 7775

32

respectively. Define

__ j aeminffi

t + ffio :

If _j 2 [0; __) for all j 2 S, then ~A is positive definite.

Proof. We have

xT ~Bx * aeminxT x

for all vectors x. Using this and the Cauchy-Schwartz inequality, we can easily establish that for _j ; j 2 S as given, we must have for x 6= 0,

xT ~Ax ? (aemin \Gamma  __(ffit + ffio))xT x * 0 as desired. Q.E.D.

The estimate obtained in the above result tends to be conservative and not very useful for practical purposes. In spite of the computational difficulty in obtaining the largest possible bound

__ for which the matrix ~A is positive definite, Proposition 2 is useful conceptually since it provides a sound basis for distinguishing the "small friction" case in which the friction coefficients at the sliding contacts are within the theoretically guaranteed bound __. From now on, we shall assume that ~A is positive definite. (Note, however, that ~A is not symmetric.) This implies that the principal submatrix 24

( ~Ann)SS (Ann)SR

( ~Ann)RS (Ann)RR

35

is P. By Theorem 4, for each fixed but arbitrary pair (cRt; cRo), there exists a unique pair (an; cn) satisfying"

aSn aRn # =

24 ( ~

Ann)SS (Ann)SR

( ~Ann)RS (Ann)RR

35 "

cSn

cRn #+"

(Ant)SR (Ano)SR (Ant)RR (Ano)RR # "

cRt cRo #+"

bSn bRn #

(an; cn) * 0 (an)T cn = 0;

(47) moreover, as a function of (cRt; cRo), (an; cn) is Lipschitz continuous. Let an(cRt; cRo) and cn(cRt; cRo) denote the unique solution function of the LCP (47). Also define

F (cRt; cRo) j 24 (

~Atn)RS (Atn)RR (Att)RR (Ato)RR

( ~Aon)RS (Aon)RR (Aot)RR (Aoo)RR

35 2666666

4

cSn(cRt; cRo) cRn(cRt; cRo)

cRt cRo

3777 7775 + " bRt

bRo # :

Let N j 2jRj. The function F maps the Euclidean space RN , which contains the pair of vectors (cRt; cRo), into itself; moreover, F is a piecewise linear map, hence Lipschitz continuous. More importantly, F has the following property.

33

Proposition 4 If the matrix ~A is positive definite, then the mapping F defined above is continuous and strongly monotone.

Proof. See Proposition 1 in [21]. Q.E.D.

We define a set-valued map K : RN ! RN as follows. For a given pair (cRt; cRo) 2 RN , let the set K(cRt; cRo) consist of all vectors (c0Rt; c0Ro) 2 RN such that for all j 2 R,

(c0jt)2 + (c0jo)2 ^ _2j (cjn(cRt; cRo))2: The pair (K; F ) defines a quasi-variational inequality problem which is to find a pair of vectors (cRt; cRo) 2 RN with two properties: (a) (cRt; cRo) 2 K(cRt; cRo); and (b) for all vectors (c0Rt; c0Ro) 2 K(cRt; cRo),

F (cRt; cRo)T " c

0Rt \Gamma  cRt

c0Ro \Gamma  cRo # * 0: We denote this problem by QVI (K; F ).

Theorem 5 Suppose that the matrix ~A is P. The problem (42) is equivalent to the QVI (K; F ) in the sense that if

(an; cn; aRt; aRo; cRt; cRo) (48)

solves (42), then (cRt; cRo) solves the QVI (K; F ); conversely, if (cRt; cRo) solves the QVI (K; F ), then

(_an; _cn; _aRt; _aRo; cRt; cRo)

solves (42), where "

_an _cn # j "

an(cRt; cRo)

cn(cRt; cRo) # and

24 _a

Rt

_aRo 35 j F (cRt; cRo):

Proof. See Theorem 3 in [21]. Q.E.D. A.4 Proof of Theorem 1(a) Theorem 5 allows us to apply an existence result in [5] to establish the existence of a solution to the multi-rigid-body model in the small-friction case. Since the proof is exactly the same as Theorem 4 in [21], we omit the details. Thus part (a) of Theorem 1 that pertains the friction cone model is proved. The proof for the friction pyramid model is similar upon a suitable modification of the definition of the set-valued map K. The details are omitted.

In the cited reference, an extension of Theorem 4 therein was obtained that permitted the system Jacobian matrix, J , to be (column) rank deficient (but still under the all-rolling condition). Presently, we do not know how to generalize Theorem 1 to the rank deficient case. Notice that part (a) of Theorem 1 assumes the "smallness" of the friction coefficients corresponding to the sliding contacts only; there is no restriction on the friction coefficients corresponding to the rolling contacts. This is consistent with the all-rolling case treated in [21] where no assumption on the friction coefficients was made.

34

A.5 Proof of Theorem 1(b) The key idea of the proof is based on the argument given by L"otstedt for the planar case [14]. However, our proof for the three-dimensional problem is somewhat more complicated, especially for the cone model which involves a set of nonlinear frictional constraints. Suppose that

(a*n; c*n; a*Rt; a*Ro; c*Rt; c*Ro) ; * = 1; 2 are two solutions of (42). Then we can rewrite the two frictional constraints (13) equivalently as

c*jt = \Gamma _j c*jn**j cos *j and c*jo = \Gamma _j c*jn**j sin *j where **j is a scalar in the interval [\Gamma 1; 1] satisfying the implication:

r*j j q(a*jt)2 + (a*jo)2 ? 0 ) **j = 1 (49) and where *j is a suitable scalar satisfying

a*jt = r*j cos *j and a*jo = r*j sin *j : For each j 2 R, we have

c1jt \Gamma  c2jt = \Gamma _j c1jn*1j cos 1j + _j c2jn*2j cos 2j = (\Gamma _j *1j cos 1j ) (c1jn \Gamma  c2jn) \Gamma  _jc2jn(*1j cos 1j \Gamma  *2j cos 2j );

similarly,

c1jo \Gamma  c2jo = (\Gamma _j *1j sin 1j ) (c1jn \Gamma  c2jn) \Gamma  _j c2jn(*1j sin 1j \Gamma  *2j sin 2j ):

Introducing the jRj \Theta  jRj diagonal matrices Dt and Do whose diagonal entries are \Gamma _j *1j cos 1j and \Gamma _j *1j sin 1j respectively, we obtain

c1Rt \Gamma  c2Rt = Dt(c1Rn \Gamma  c2Rn) \Gamma  _Rc2Rn[*1R cos 1R \Gamma  *2R cos 2R] c1Ro \Gamma  c2Rt = Do(c1Rn \Gamma  c2Rn) \Gamma  _Rc2Rn[*1R sin 1R \Gamma  *2R sin 2R]

where the notation in the second terms of the right-hand side of these equations has an obvious (componentwise) meaning. Consequently, we have2

66666 64

a1Sn \Gamma  a2Sn a1Rn \Gamma  a2Rn

a1Rt \Gamma  a2Rt a1Ro \Gamma  a2Ro

3777 7775 = ~A

2666 6664

c1Sn \Gamma  c2Sn c1Rn \Gamma  c2Rn

c1Rt \Gamma  c2Rt c1Ro \Gamma  c2Ro

3777 7775

= ~A

2666 6664

I 0 0 0 0 I 0 0 0 Dt I 0 0 Do 0 I

3777 7775

2666 6664

c1Sn \Gamma  c2Sn c1Rn \Gamma  c2Rn _Rc2Rn[*2R cos 2R \Gamma  *1R cos 1R]

_Rc2Rn[*2R sin 2R \Gamma  *1R sin 1R]

3777 7775 :

35

Using the representation (46) for ~A, we see that the matrix

~A

2666 6664

I 0 0 0 0 I 0 0 0 Dt I 0 0 Do 0 I

3777 7775 (50)

is a small perturbation of the matrix ~B in (44) provided that all friction coefficients are sufficiently small. Thus, the matrix (50) is positive definite under the assumptions of Theorem 1(b). We claim26

66666 4

c1Sn \Gamma  c2Sn c1Rn \Gamma  c2Rn _Rc2Rn[*2R cos 2R \Gamma  *1R cos 1R]

_Rc2Rn[*2R sin 2R \Gamma  *1R sin 1R]

3777 7775

T 2666

6664

a1Sn \Gamma  a2Sn a1Rn \Gamma  a2Rn

a1Rt \Gamma  a2Rt a1Ro \Gamma  a2Ro

3777 7775 ^ 0: (51)

By the nonnegativity of (a*n; c*n) and the complementarity (a*n)T c*n = 0, for * = 1; 2, it is easy to deduce "

c1Sn \Gamma  c2Sn

c1Rn \Gamma  c2Rn #

T " a1

Sn \Gamma  a2Sn

a1Rn \Gamma  a2Rn # ^ 0: Consider an index j 2 R. We have

[*2j cos 2j \Gamma  *1j cos 1j ](a1jt \Gamma  a2jt) + [*2j sin 2j \Gamma  *1j sin 1j ](a1jo \Gamma  a2jo)

= [*2j cos 2j \Gamma  *1j cos 1j ][r1j cos 1j \Gamma  r2j cos 2j ] + [*2j sin 2j \Gamma  *1j sin 1j ][r1j sin 1j \Gamma  r2j sin 2j ] = \Gamma r1j *1j \Gamma  r2j *2j + (r1j *2j + r2j *1j ) cos(1j \Gamma  2j ) = \Gamma r1j \Gamma  r2j + (r1j *2j + r2j *1j ) cos(1j \Gamma  2j ) ^ 0 where the last equality follows from the implication (49) and the last inequality follows because j*1;2j j ^ 1. Consequently, the claim (51) is validated. This implies that2

666666

4

c1Sn \Gamma  c2Sn c1Rn \Gamma  c2Rn _Rc2Rn[*2R cos 2R \Gamma  *1R cos 1R]

_Rc2Rn[*2R sin 2R \Gamma  *1R sin 1R]

3777 7775

T

~A

2666 6664

I 0 0 0

0 I 0 0 0 Dt I 0 0 Do 0 I

3777 7775

2666 6664

c1Sn \Gamma  c2Sn c1Rn \Gamma  c2Rn _Rc2Rn[*2R cos 2R \Gamma  *1R cos 1R]

_Rc2Rn[*2R sin 2R \Gamma  *1R sin 1R]

3777 7775

is nonpositive. By the positive definiteness of the matrix (50), it follows that c1n = c2n and

c1Rt = \Gamma _Rc1Rn*1R cos 1R = \Gamma _Rc2Rn*1R cos 1R = \Gamma _Rc2Rn*2R cos 2R = c2Rt; similarly, c1Ro = c2Ro. Consequently, we havei

a1n; c1n; a1Rt; a1Ro; c1Rt; c1Roj = ia2n; c2n; a2Rt; a2Ro; c2Rt; c2Roj ; thereby establishing Theorem 1(b) pertaining to the friction cone model. The proof pertaining to the friction pyramid model is slightly simpler and requires no new argument. The details are omitted.

36

A.6 Proof of Proposition 1 Suppose that (48) is a solution of the friction pyramid model that satisfies the quadratic cone constraints (10). It suffices to verify (13). Consider an arbitrary j 2 R. There is nothing to prove if ajt = ajo = 0. Suppose that both ajt and ajo are nonzero. From (20), we deduce

c2jt = c2jo = _2j c2jn: Since (10) holds, we must have

_j cjn = cjt = cjo = 0:

Thus (13) hold. Suppose that ajt is zero and ajo is nonzero. Then cjt 2 [\Gamma _jcjn; _jcjn] and c2jo = _2j c2jn. Again since (10) holds, it follows that cjt = 0; thus (13) hold. The remaining case where ajt is nonzero and ajo is zero is proved similarly. This completes the proof of the proposition.

A.7 Proof of Theorem 2 We cite the following result which gives a sufficient condition for Lemke's algorithm to successfully terminate with a solution of a general LCP (r; M ).

Proposition 5 Let M be a copositive matrix. If the following implication holds:

x * 0 M x * 0 xT M x = 0

9??= ??; ) r

T x * 0

then under a standard nondegeneracy assumption, Lemke's algorithm will successfully compute a solution to the LCP (r; M ) in a finite number of pivots.

As we see from Section 3, the three-dimensional multi-rigid-body contact problem with friction pyramid law can be formulated as the LCP (r; M) where

r j

2666 66666 66664

bSn bRn

bRt bRo

0 0

3777 77777 77775 and M j

2666 66666 66664

(M nn)SS (M nn)SR (Ant)SR (Ano)SR 0 0 (M nn)RS (M nn)RR (Ant)RR (Ano)RR 0 0

(M tn)RS (M tn)RR (Att)RR (Ato)RR I 0 (M on)RS (M on)RR (Aot)RR (Aoo)RR 0 I

0 2U R \Gamma I 0 0 0 0 2U R 0 \Gamma I 0 0

3777 77777 77775 : (52)

To establish Theorem 2, it suffices to verify the two requirements in Proposition 5. For the copositivity of M , we note that

M = M 1 + M 2 + M 3

37

where

M 1 j

2666 66666 66664

(M nn)SS (M nn)SR (Ant)SR (Ano)SR 0 0 (M nn)RS (Mnn)RR (Ant)RR (Ano)RR 0 0

(M tn)RS (M tn)RR (Att)RR (Ato)RR 0 0 (M on)RS (M on)RR (Aot)RR (Aoo)RR 0 0

0 0 0 0 0 0 0 0 0 0 0 0

3777 77777 77775

M 2 j

2666 66666 66664

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 I 0 0 0 0 0 0 I 0 0 \Gamma I 0 0 0 0 0 0 \Gamma I 0 0

3777 77777 77775 and M 3 j

2666 66666 66664

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2U R 0 0 0 0 0 2U R 0 0 0 0

3777 77777 77775 :

The matrix M 1 is positive semidefinite by the smallness assumption of the friction coefficients; M 2 is skew-symmetric, thus positive semidefinite; and M 3 is nonnegative, thus copositive. Consequently, M is copositive.

Now suppose we have a vector x * 0 such that xT M x = 0. Then we must have xT M 1x = 0. We may write

x = (xSn; xRn; xRt; xRo; uR; vR):

Let us define the principle submatrix ^M of M (and M 1) as follows:

^M j

2666 6664

(Mnn)SS (M nn)SR (Ant)SR (Ano)SR (M nn)RS (M nn)RR (Ant)RR (Ano)RR

(M tn)RS (Mtn)RR (Att)RR (Ato)RR (Mno)RS (M on)RR (Aot)RR (Aoo)RR

3777 7775 : (53)

Then, by the structure of M 1 and the positive definiteness of ^M , it follows that

(xSn; xRn; xRt; xRo) = 0: Hence rT x = 0. Q.E.D.

B Details of The Interior Point Method We present the details of an interior point method for solving the general LCP (r; M) satisfying certain assumptions; in the section on numerical results, we apply this method to the particular

38

LCP that arises from the three-dimensional multi-rigid-body contact problem with Coulomb friction pyramid law.

We refer the reader to the monograph [12] for some background results and detailed historical discussion of the family of interior point methods for solving the LCP.

B.1 Description and convergence analysis Consider an LCP (r; M ) where the vector r 2 Rn and matrix M 2 Rn\Theta n. Let

F j f(w; z) 2 R2n+ : w = r + M zg denote the feasible region of this problem. We write

F+ j F " R2n++ where R2n++ j f(w; z) 2 R2n : (w; z) ? 0g. We postulate the following assumptions on the pair (r; M):

(a) M is a P0-matrix; (b) the LCP (r; M ) is strictly feasible; that is, F+ 6= ;;

(c) for any scalars ff ? fi ? 0, the level set L(ff; fi) j f(w; z) 2 F+ : ffe * w ffi z * fieg is

bounded, where e denotes the n-vector of all ones and u ffi v denotes the Hadamard product of two vectors u and v, that is, u ffi v is the vector whose i-th component is equal to uivi for all i.

It is useful to point out that these assumptions are very common in the family of interior point methods for solving the LCP [12]. Nevertheless assumption (c) is rather distinctive in that it does not imply that the solution set of the LCP (r; M) is necessarily bounded. In the cited reference, there is a similar set of three conditions labeled (i), (ii), and (iii) on page 8 which are parallel to ours. The only difference is our assumption (c) which is weaker than condition (iii) in [12]. (For one thing, the three assumptions therein imply the boundedness of the solution set of the LCP.) It turns out that for the LCP of interest to us, assumption (c) is satisfied whereas condition (iii) fails to hold.

The interior point method for solving the LCP (r; M) satisfying the conditions (a), (b), and (c) starts at a strictly feasible pair (w0; z0) 2 F+. In general, given a pair of vectors (wk; zk) 2 F+, the method consists of two major computational steps. Let fik 2 [0; 1) be a given scalar, called the centering parameter. This parameter plays an extremely important role in the practical success of the interior point method; we refer to [12] for a full discussion of this role. Let

,k j (w

k)T zk

n : We solve the following system of linear equations to compute the direction vector (dwk; dzk):24

\Gamma I M

Zk W k

35 24

dwk

dzk

35

= " 0wk ffi zk \Gamma  fi

k,ke #

(54)

39

where W k = diag(wk) and Zk = diag(zk). By the P0-property of M and the positivity of (wk; zk), it can be shown that the matrix on the left side of (54) is nonsingular; thus (dwk; dzk) is well defined. Next, we need to introduce the merit function for the method, which is defined as

(w; z) j (n + i) log wT z \Gamma 

nX

i=1

log wizi; for (w; z) ? 0;

where i is a positive constant. Again, we refer to [12] for the motivation and derivation of this function. The vector (dwk; dzk) computed above turns out to be a descent direction for the function  at the current iterate (wk; zk). Hence, a line search can be executed on this function starting at this iterate and moving along the generated direction. This search is a standard routine in a nonlinear programming algorithm. It requires two fixed constants, ae; oe 2 (0; 1), where ae controls the step size and oe controls the decrease of .

The following is a detailed description of the interior point method for solving the LCP (r; M) satisfying conditions (a), (b), and (c).

The interior point method Step 0. (Initialization) Let "; oe; ae; j 2 (0; 1) be given scalars. Let i ? 0, (w0; z0) 2 F+, and fi0 2 [0; 1) be arbitrary. Set k = 0.

Step 1. (Computing search direction) Solve the system of linear equations (54) to obtain (dwk;dzk). Step 2. (Computing step size) Calculate

_*k j min (1; (\Gamma  w

ki

dwki : dw

k i ! 0) ; (\Gamma  z

ki

dzki : dz

k i ! 0))

and let *0k j j_*k. Let mk be the smallest nonnegative integer m such that

(wk + aem*0kdwk; zk + aem*0kdzk) \Gamma  (wk; zk) ^ \Gamma oe(1 \Gamma  fik)iaem*0k: (55) Set *k j aemk *0k. Step 3. (Update and termination check) Set *k j aemk *0k and

(wk+1; zk+1) j (wk + *kdwk; zk + *kdzk): If (wk+1)T zk+1 ^ ", terminate; the pair (wk+1; zk+1) is a desired approximate solution of the LCP (r; M ). Otherwise, choose fik+1 2 [0; fik]; replace k by k + 1 and return to Step 1.

We give some further explanation to Step 2. The scalar _*k is the largest step size * 2 (0; 1] for which (wk + *dwk; zk + *dzk) is nonnegative; thus by scaling _*k by the factor j 2 (0; 1), we are ensured that (wk + *dwk; zk + *dzk) is positive for all * 2 [0; *0k]. The integer mk can be determined by starting with m = 0 and increasing m by 1 each time the inequality (55) fails to hold; it can be shown that in a finite number of trials, the desired integer mk can be obtained.

40

The tolerance " is used in practical implementation to check the successful termination of the method. In theory, an infinite sequence of iterates f(wk; zk)g is generated by the method. The following result summarizes the main properties of this sequence and shows that the LCP (r; M) must have a solution. Since this result has not appeared in the LCP literature, we give a detailed proof (except for the omission of a few minor details).

Theorem 6 Under assumptions (a), (b), and (c), the above method generates a well-defined sequence of iterates f(wk; zk)g having the following properties:

(i) f(wk; zk)g ae F+; (ii) f(wk; zk)g is strictly decreasing; (iii) limk!1(wk)T zk = 0;

(iv) every accumulation point of the sequence f(wk; zk)g, if it exists, is a solution of the LCP

(r; M).

Moreover, the LCP (r; M) must have a solution. Proof. The first two parts (i) and (ii) are fairly easy to see. Part (iv) follows easily from (iii). The proof of part (iii) and the final conclusion is quite technical. In what follows, we sketch the principal idea of the proof and leave out some details.

To prove part (iii), we note that for an arbitrary pair (w; z) ? 0, the arithmetic-geometric mean inequality yields

(w; z) = i log wT z + log (P

n i=1 wizi)

nQ

n i=1 wizi * i log

wT z:

Thus the fact that f(wk; zk)g is strictly decreasing implies that f(wk)T zkg is bounded above, say, by the positive scalar ff. Since (wk; zk) is positive, it follows that

wk ffi zk ^ ffe; for all k: Now suppose that for some subsequence f(wk; zk) : k 2 ^g, we have

lim inf k(2^)!1(w

k)T zk ? 0; (56)

which, by the fact that f(wk; zk)g is bounded above, implies that for all i = 1; : : :; n,

lim inf k(2^)!1 w

k i z

k i ? 0:

Hence for some positive scalar fi ! ff, we have

f(wk; zk) : k 2 ^g ae L(ff; fi):

41

By condition (c), the subsequence f(wk; zk) : k 2 ^g is therefore bounded. Without loss of generality, we may assume that f(wk; zk) : k 2 ^g converges, say, to (w\Lambda ; z\Lambda ) which must belong to F+ . The entire sequence f(wk; zk)g must converge. It follows from (55) that

lim k!1 ae

mk *0

k = 0:

The sequence of matrices 8!

:24

\Gamma I M

Zk W k

35

: k 2 ^9=;

converges to the matrix "

\Gamma I M

Z\Lambda  W \Lambda  # (57) where Z\Lambda  j diag(z\Lambda ) and W \Lambda  j diag(w\Lambda ); the limit matrix (57) is nonsingular, by the P0-property of M and the positivity of (w\Lambda ; z\Lambda ). Hence the subsequence of directions f(dwk; dzk) : k 2 ^g is bounded. We can deduce from this conclusion that

lim inf k(2^)!1 _*k ? 0:

Hence we must have

lim k(2^)!1 ae

mk = 0;

or equivalently,

lim k(2^)!1 mk = 1:

Using this limit condition and the definition of mk, we can establish a contradiction; see [30]. Consequently, there exists no subsequence f(wk; zk) : k 2 ^g for which (56) is valid; hence (iii) holds.

We now use the conclusion (iii) to show that the LCP (r; M ) must have a solution. The difficulty with the proof is that we do not know whether the sequence f(wk; zk)g is actually bounded; hence we need to devise an argument which does not make use of this boundedness condition.

There exists a partition of the index set f1; : : :; ng into the disjoint union of two index subsets I and I0 such that for each i 2 I, the sequence fzki g is bounded and for each j 2 I0, the sequence fzkj g is unbounded. Without loss of generality, we may assume, by working with appropriate subsequences, that for each i 2 I,

lim k!1 z

k i = ~zi

for some ~zi * 0, and for each j 2 I0,

lim k!1 z

k j = 1:

By (iii), it follows that for each j 2 I0,

lim k!1 w

k j = 0:

42

We further partition I into the disjoint union of two index subsets:

I+ j fi 2 I : ~zi ? 0g; and I0 j fi 2 I : ~zi = 0g: By (iii) again, we have for each i 2 I+,

lim k!1 w

k i = 0:

Let K be the union of I0 and I+; let _K be the complement of K in f1; : : :; ng. The equation wk = r + M zk can be equivalently written as follows:

r + M I zkI \Gamma  IKwkK = I _Kwk_K \Gamma  M I0zkI0 where AL denotes the columns of a matrix A that are indexed by L. The left-hand vector in the above equation converges to r + M I ~zI . By the closedness of a polyhedron, it follows that there exist nonnegative vectors ~w _K and ~zI0 such that

r + M I ~zI = I _K ~w _K \Gamma  M I0 ~zI0 : The vector ~z is easily seen to be a solution of the LCP (r; M). Q.E.D.

Due to property (i) in the above theorem, the method we have described is called a feasible interior point method. One of its special features is that it requires a pair (w0; z0) which is strictly feasible to the LCP (r; M). As we shall see in the next subsection, such a pair can be obtained without too much difficulty for the multi-rigid-body contact problem with friction pyramid law, by means of a Phase I simplex-type procedure. There is a family of infeasible interior point methods which circumvent the computation of such an initial pair and can be initiated at any positive pair (w0; z0) that is not necessarily feasible to the LCP (r; M); see [30]. At this time, we are not sure whether such infeasible interior point methods will be convergent under the assumptions (a), (b), and (c) in general, and for the contact LCP in particular. Computationally, it is expected that an infeasible method will be at least comparable with the feasible method; however, we have not carried out the experiments for the contact LCP with an infeasible interior point method.

B.2 Application to the friction pyramid model We recall the formulation of the three-dimensional multi-rigid-body contact problem with friction pyramid law as the LCP (r; M) where r and M are given in (52). Recall the definition of ^M given in equation (53). Similar to Proposition 2, we can establish the following result.

Proposition 6 Suppose that the system Jacobian matrix, J (defined in equation (25)), has full column rank. There exists a positive scalar __ such that if _j 2 [0; __] for all j = 1; : : :; nc, the matrix

^M is positive definite.

Our goal in this subsection is to show that under the assumptions of Proposition 6, the conditions (a), (b), and (c) required for the convergence of the feasible interior method are satisfied for the

43

LCP (r; M) with (r; M) given by (52). The demonstration is not particularly difficult but it involves some additional notation and algebraic argument that is not very illuminating.

As a first step to show that the matrix M in (52) is P0, we need to derive an extension of Proposition 6 that deals with matrices similar to ^M . For this purpose, we recall the matrix ~B defined in (44) and introduce certain matrices similar to Et and Eo in (45). Specifically, let \Sigma  denote a diagonal matrix of order jRj with diagonal entries equal to 0 or \Sigma 1. Define

F t j

2666 6664

0 (Ant)SR\Sigma U R 0 0 0 (Ant)RR\Sigma U R 0 0 0 (Att)RR\Sigma U R 0 0 0 (Aot)RR\Sigma U R 0 0

3777 7775 and F o j

2666 6664

0 (Ano)SR\Sigma U R 0 0 0 (Ano)RR\Sigma U R 0 0 0 (Ato)RR\Sigma U R 0 0 0 (Aoo)RR\Sigma U R 0 0

3777 7775 :

Proposition 7 Suppose that the system Jacobian matrix, J , (defined in equation (25)), has full column rank. There exists a positive scalar __ such that if _j 2 [0; __] for all j = 1; : : :; nc, the matrix

~B \Gamma  Et \Gamma  F t \Gamma  Eo \Gamma  F o (58)

is positive definite, for any diagonal matrix \Sigma  with diagonal entries equal to 0 or \Sigma 1.

Proof. This follows easily from a continuity argument and the fact that there are finitely many diagonal matrices \Sigma  of the type specified. Q.E.D.

We note that with \Sigma  being the identity matrix, the matrix in (58) reduces to the matrix ^M .

Proposition 8 Suppose that the matrix, J , (defined in equation (25)), has full column rank. There exists a positive scalar __ such that if _j 2 [0; __] for all j = 1; : : :; nc, the matrix M defined in (52) is P0.

Proof. We need to show that if each _j * 0 is sufficiently small, then all principal minors of M are nonnegative. The proof consists of examining the principal submatrices of M and showing that these submatrices either are positive semidefinite or, by means of a determinantal argument, have nonnegative determinants. Throughout the proof, the matrix ^M is assumed to be positive definite. We label the columns of M according to the variables with which they are associated; that is:

(cSn; cRn; s+Rt; s+Ro; a\Gamma Rt; a\Gamma Ro): The principal submatrices of M are of several types:

(I) a principal submatrix of ^M ; (II) a principal submatrix of ^M together with some columns (and same rows) corresponding to

either some of the a\Gamma Rt or a\Gamma Ro variables but not both;

(III) a principal submatrix of ^M together with some columns (and same rows) corresponding to

at least one a\Gamma Rt-variable and at least one a\Gamma Ro-variable;

44

(IV) a zero matrix. Clearly, we need to be concerned only with the principal submatrices of the middle two types. We consider a principal submatrix of M of type (II) of the following form:2

66666 64

(M nn)S0S0 (M nn)S0R1 (Ant)S0R2 0 (Mnn)R1S0 (M nn)R1R1 (Ant)R1R2 0

(M tn)R2S0 (M tn)R2R1 (Att)R2R2 IR2R3

0 2(U R)R3R1 \Gamma IR3R2 0

3777 7775 (59)

where S0 ` S and Ri for i = 1; 2; 3 are (not necessarily disjoint) subsets of R with R3 6= ;. Clearly, if R2 6= R3, then the matrix (59) would contain a column of zeros; thus its determinant is equal to zero. If R2 = R3, then since"

(Att)R2R2 IR2

\Gamma IR2 0 #

\Gamma 1

= " 0 \Gamma I

R2

IR2 (Att)R2R2 #

we can compute the Schur complement of the submatrix"

(Att)R2R2 IR2

\Gamma IR2 0 #

in the matrix (59); this Schur complement turns out to be"

(M nn)S0S0 (M nn)S0R1 + 2(Ant)S0R2(U R)R2R1 (M nn)R1S0 (Mnn)R1S0 + 2(Ant)R1R2 (U R)R2R1 #

which is a principal submatrix of a matrix of the type (58). Thus this Schur complement is positive definite, provided that the friction coefficients are sufficiently small. Since

det " (A

tt)R2R2 IR2

\Gamma IR2 0 # = 1

by the Schur determinantal formula, the matrix (59) has a positive determinant.

By slightly modifying the above argument, we can establish that all principal submatrices of M of types II and III have nonnegative determinants. Consequently, M is a P0-matrix. Q.E.D.

The proof of the strict feasibility of the LCP (r; M) turns out to be much easier, provided that all friction coefficients are positive (as opposed to just nonnegative) and sufficiently small.

Proposition 9 Suppose that the system Jacobian matrix, J (defined in equation (25)), has full column rank. There exists a positive scalar __ such that if _j 2 (0; __] for all j = 1; : : :; nc, then there is a vector x ? 0 such that M x ? 0, where M is defined in (52). Consequently, for any vector r, the LCP (r; M) has a strictly feasible solution.

45

Proof. By scaling, it suffices to establish the first assertion. Let eR be the vector of all ones of dimension jRj. Since ^M is positive definite under the assumptions, it follows that there exists (xSn; xRn) ? 0 such that

(M nn)SS xSn + (M nn)SRxRn + (Ant)SReR + (Ano)SReR ? 0 (M nn)RS xSn + (Mnn)RRxRn + (Ant)RReR + (Ano)RReR ? 0 2U RxRn \Gamma  eR ? 0:

To complete the proof, simply choose (uR; vR) ? 0 large enough so that

(M tn)RS xSn + (M tn)RRxRn + (Att)RReR + (Ato)RReR + uR ? 0 (M on)RSxSn + (M on)RRxRn + (Aot)RReR + (Aoo)RReR + vR ? 0:

The desired vector x is then

(xSn; xRn; eR; eR; uR; vR):

This completes the proof. Q.E.D.

Finally, we verify that condition (c) required in the feasible interior point method is satisfied for the pair (r; M) given by (52).

Proposition 10 Suppose that the system Jacobian matrix, J , (defined in equation (25)), has full column rank. There exists a positive scalar __ such that if _j 2 (0; __] for all j = 1; : : : ; nc, condition (c) holds for the pair (r; M ) given by (52).

Proof. Let f(wk; zk)g be an arbitrary sequence in the level set L(ff; fi) for some ff ? fi ? 0. Consistent with the notation for the rigid-body model, we write

wk j

2666 666666 66664

akSn akRn ak;+Rt ak;+Ro

sk;\Gamma Rt sk;\Gamma Ro

3777 777777 77775

and zk j

2666 666666 66664

ckSn ckRn sk;+Rt sk;+Ro ak;\Gamma Rt ak;\Gamma Ro

3777 777777 77775

:

Since fwk ffi zkg is bounded, it follows that for all k,2

66666 64

ckSn ckRn sk;+Rt sk;+Ro

3777 7775

T

^M

2666 6664

ckSn ckRn sk;+Rt sk;+Ro

3777 7775 + (sk;+Rt )T ak;\Gamma Rt + (sk;+Ro )T ak;\Gamma Ro ^ constant (60)

46

which implies, by the positive definiteness of ^M and the nonnegativity of the variables, that f(ckSn; ckRn; sk;+Rt ; sk;+Ro )g is bounded. Thus fakRng is also bounded.

It remains to show that f(ak;\Gamma Rt ; ak;\Gamma Ro )g is bounded. Suppose that for some j 2 R, we have

lim k!1 a

k;\Gamma  jt = 1:

Then (60) implies

lim k!1 s

k;+ jt = 0:

Since

ff * ak;\Gamma jt i2_jckjn \Gamma  sk;+jt j * fi;

and _j ? 0, it follows that

lim k!1 c

k jn = 0:

By the boundedness of fakjng, we deduce

lim k!1 a

k jnc

k jn = 0:

But this is impossible because for all k,

akjnckjn * fi ? 0: Consequently, the sequence fzkg, and thus fwkg, is bounded. Q.E.D.

By combining Propositions 8, 9, and 10, we conclude that Theorem 3 holds.

47