UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Efficient and transparent haptic rendering of rigid body motion with constraints Constantinescu, Daniela 2005

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2005-994554.pdf [ 14.4MB ]
Metadata
JSON: 831-1.0064992.json
JSON-LD: 831-1.0064992-ld.json
RDF/XML (Pretty): 831-1.0064992-rdf.xml
RDF/JSON: 831-1.0064992-rdf.json
Turtle: 831-1.0064992-turtle.txt
N-Triples: 831-1.0064992-rdf-ntriples.txt
Original Record: 831-1.0064992-source.json
Full Text
831-1.0064992-fulltext.txt
Citation
831-1.0064992.ris

Full Text

Efficient and Transparent Haptic Rendering of Rigid Body Motion with Constraints by Daniela Constantinescu B.Sc.  (Mechanical Engineering), University of Brasov, 1986  M.A.Sc. (Mechanical Engineering), University of British Columbia, 1998  A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F THE REQUIREMENTS D O C T O R  FOR THE DEGREE O F  OF  PHILOSOPHY  in THE FACULTY OF GRADUATE  STUDIES  DEPARTMENT OF ELECTRICAL A N D COMPUTER  T H E UNIVERSITY OF BRITISH  ENGINEERING  COLUMBIA  December 2004 © Daniela Constantinescu, 2004  Abstract This work presents general purpose simulation and control techniques for efficient and transparent haptic rendering of rigid body motion with constraints.  Transparent interaction is achieved by enabling users to  feel collisions and to manipulate both virtual objects and virtual linkages. Efficient rendering is accomplished through fast approximations of rigid body interaction implemented in a local model with haptic performance. Hapic rendering of impacts is based on a new representation of rigid body contact. In this representation, contacts are infinitely stiff when they arise and have limited stiffness thereafter. Multiple impacts are resolved simultaneously, in a manner consistent with conservation of energy principles and with the force capabilities of the haptic device. Haptic rendering of impacts is beneficial in training simulators for dental procedures and bone surgeries, as well as C A D and virtual prototyping systems with force feedback. Realistic linkage manipulation is enabled by permitting users to operate linkages from any link and through singularities while restricting their motion according to the virtual environment geometry and the linkage topology. Linkage topology is imposed on users through penalizing users' departure from the configuration manifold of the virtual linkage. Operation of links with insufficient degrees of freedom is important in applications like training for laparoscopy, where the scope limits the tool motion at the entry point. Efficient rendering of rigid body motion with constraints is enabled by interfacing the device to a simulation through a local model of interaction. The model comprises constraints imposed on the virtual tool by virtual objects within an e-active neighborhood of the virtual tool and a dynamic proxy of the virtual tool. This model is the first that can be used to constrain both the translation and the rotation of the device and to add realistic forces to virtual environments generated using any commercial simulation package with interactive performance. The model is beneficial in consumer-grade haptic applications, because it allows easy development of haptic applications by users without detailed haptic knowledge. It can also be used to enable cooperative haptic manipulations in applications that involve two-handed operations and/or multiple users.  11  Table of Contents Abstract  ii  Table of Contents  iii  List of Tables  vi  List of Figures  vii  Glossary  xii  Acknowledgements 1  2  x  v  Introduction  1  1.1  Objective  4  1.2  Contributions  4  1.3  Organization of the thesis  7  Physically-based rigid multibody simulations  10  2.1  Interactive simulations  10  2.1.1  Coordinates for rigid multibody dynamics  13  2.1.2  Dynamic models of rigid multibody systems  14  2.1.2.1  Penalty-based simulations  15  2.1.2.2  Constraint-based simulations  16  2.1.2.3  Impulse-based simulations  18  2.2  Impact simulation  19  2.3  Haptic interaction within rigid multibody virtual environments  22  2.3.1  Haptic controllers  23  2.3.2  Local models of interaction  24  2.3.3  Haptic rendering of impact  25  2.3.4  Haptic rendering of friction  26  2.3.5  Summary  26 iii  CONTENTS  iv  3  Experimental setup  28  4  Haptic rendering of rigid contact  32  4.1  Synopsis of the impulse-augmented penalty simulation and of the haptic controller  33  4.2  The contact model  35  4.3  Resting contact  37  4.3.1  39  4.4  Colliding contact  46  4.4.1  Passive collision resolution  48  4.4.1.1  Independent constraints  49  4.4.1.2  Overdetermined constraints  52  4.4.2 4.5  4.6 5  6  Friction modeling  Rendering collisions to users  54  The performance of the impulse-augmented penalty simulation  54  4.5.1  Best performance  56  4.5.2  Accounting for device limitations  57  4.5.3  Experimental validation  65  Discussion  69  Manipulation of serial linkages  70  5.1  Control architecture for realistic linkage manipulation  71  5.2  Linkage contacts at user's hand  72  5.3  Haptic rendering of linkage inertia  74  5.4  Haptic rendering of linkage topology  76  5.4.1  Simulations  78  5.4.2  Experiments  81  Efficient haptic rendering of rigid body motion with constraints  85  6.1  Synopsis of the local model  86  6.2  Local geometry  87  6.2.1  Active geometry  87  6.2.2  Local proxy deformation  91  6.2.3  e-active local geometry  92  6.2.4  Contact clustering and constraint coherence  95  6.3  6.4  Local dynamics  98  6.3.1  Local contact interactions applied to users  98  6.3.2  Local control  99  Experiments  100  6.4.1  Passive collision dynamics in the impulse-augmented penalty local model  100  6.4.2  Transparency in the local model of interaction  102  CONTENTS  7  v  6.4.3  The effect of the teleoperation controller on the achievable contact stiffness  108  6.4.4  Dynamic contact  109  Conclusion  112  7.1  Summary of Contributions  112  7.2  Future Work  114  Bibliography  117  A Proof of implicit elimination of constraint overdeterminancy  126  B Weighted generalized inverses for linkage manipulation B . l Weighted generalized inverses of matrices  128 128  B.2 Weighted generalized inverses applied to robotics  129  C Equality of null spaces of A^" and  132  D The dynamics and control of the planar haptic interface D . l Interface actuation and dynamics D.2 Interface control  134 134 136  1  List of Tables 4.1  5.1  Parameters of the teleoperation controller and the unilateral coupler connecting the virtual environment to the haptic device  56  Parameters of the three links planar virtual linkage operated by the user in the simulations. . .  78  vi  List of Figures  1.1  The key components of a haptic interaction system and their conflicting requirements for stability and realism  1.2  2  Alleviation of the effect of the computational delay and the geometry-dependent environment stiffness via control  3  1.3  Elimination of the computational delay of the simulation via new collision detection algorithms.  3  1.4  Mitigation of the effect of the computational delay of the simulation via a local model of interaction.  4  1.5  The contact model the enables users to feel collisions: contacts have infinite stiffness upon contact and finite stiffness during contact. Users feel impulsive forces upon contact and penalty and friction forces during contact  1.6  5  Users perceive topological constraints through penalties applied to the device when their hand leaves the configuration manifold of the linkage (vh is the position and orientation of the user's hand on the virtual linkage, while h is their real position and orientation)  1.7  6  The local model of rigid body interaction includes a dynamic proxy of the virtual tool and active and predicted motion constraints imposed on the virtual tool by nearby objects  7  2.1  The simulation loop and the scope of the survey: prior work in dynamic response  11  2.2  Coordinates used to describe the system state in rigid multibody simulations  13  2.3  The contact states in a penalty-based rigid multibody simulation  15  2.4  A rectangular object interacting with a virtual wall. Each surface contact is represented through two point contacts  16  2.5  The contact states in a constraint-based rigid multibody simulation  17  2.6  The corner law  17  2.7  The contact states in an impulse-based rigid multibody simulation  18  2.8  The virtual coupler  23  3.1  The planar haptic simulation system used for experiments throughout this thesis  28  3.2  Communication between the virtual environment simulation and the force control loop in the  3.3  planar haptic interaction system used in this thesis  30  Testbed virtual environments used in the experiments  31  vii  LIST OF FIGURES  viii  4.1  Haptic manipulation in an impulse-augmented penalty virtual world  33  4.2  The impulse-augmented penalty-based rigid body contact model: contact rigidity is exactly enforced upon contact and only approximately enforced during contact  35  4.3  Contact information received from collision detection  35  4.4  Example transition from free motion to resting contact during a fixed step haptic simulation, cu is the angular velocity of the body, v is the velocity of its center of mass,  is the velocity of  point A on the body, and t is the time step of the simulation 4.5  36  Example contact forces and and hand wrenches arising during the haptic manipulation of a contact group with dynamics computed in configuration space  38  4.6  Simulated peg-in-hole manipulation  40  4.7  Peg-in-hole simulation. Friction is simulated using the reset-integrator model [61] with various damping coefficients  41  4.8  Classical approximation of Coulomb friction  43  4.9  Peg-in-hole simulation. Friction is simulated using the critically damped reset-integrator model [61], the model proposed by Nahvi et al. [113], the modified Dahl model [67], and the classical dry friction approximation  44  4.10 Forces felt by users and C O M trajectories of a virtual peg sliding with friction in a tight hole due to a sinusoidal force applied to its C O M . Various friction models are used to simulate dry friction. Note that all models render slip and stick. Perceivable oscillations occur during stiction when friction is simulated using the modified Dahl model [67] (Figure 4.10(b)) and the model based on human finger-pad characteristics [113] (Figure 4.10(c)). Some slip to stick transitions are missed by the model based on the human finger-pad characteristics [113] (Figure 4.10(c)). .  45  4.11 Two link planar manipulator whose loss of kinetic energy for various values of the coefficient of restitution is shown in Figure 4.12 for two contact geometries. The simplifying assumption in Equation (4.24) imposes the particular solution shown in Figure 4.12(b)  51  4.12 Loss of kinetic energy of the planar two link manipulator in Figure 4.11 during one frictionless collision for two contact geometries  51  4.13 Planar virtual world used in simulations and experiments. Starting from rest and the position shown, the rectangular object is pushed into the lower right corner by a controlled constant force. 55 4.14 1 D O F mechanical equivalents of the controllers connecting the haptic device and the virtual environment. Controller parameters are given in Table 4.1  56  4.15 Simulated hand trajectories obtained when constraint-based ("ideal"), impulse-augmented penaltybased ("IAPB"), and penalty-based ("PB") interactions are applied to users by the four channel controller, and when constraint-based interactions are applied by the unilateral coupler ( " C B " ) . UC  The device applies the simulated impulses to users in one step  58  L I S T OF  ix  FIGURES  4.16 Simulated forces along the x-axis and torques applied to users by the four channel controller when the virtual world is generated using the constraint-based ("ideal"), the impulse-augmented penalty-based ("IAPB"), and the penalty-based ("PB") methods, and applied by the unilateral coupler when the virtual environment is generated using the constraint-based method ( " C B " ) . UC  The device applies the simulated impulses to users in one step  59  4.17 Kinetic energy of the users' hand when constraint-based ("ideal"), impulse-augmented penaltybased ("IAPB"), and penalty-based ("PB") interactions are appllied to users by the four channel controller, and when constraint-based interactions are applied by the unilateral coupler ( " C B " ) . UC  The device applies the simulated impulses to users in one step  60  4.18 Simulated hand trajectories during user interaction with the impulse-augmented penalty-based virtual world when the device limitations are ignored and when they are taken into account ("ideal" - full constraint-based collision impulses are applied to users; "p/ (i" - full collision u  impulses are applied to users; "Pscaied" -  "p ,turated" sa  -  collision impulses are saturated on the device;  collision impulses are scaled in the virtual environment;  "p read" sp  -  collision impulses  are spread over several simulation steps when necessary)  62  4.19 Simulated hand trajectories obtained when constraint-based ("ideal"), impulse-augmented penaltybased ("IAPB L M " ) , and penalty-based ("PB") interactions are transmitted to users by a four channel controller, and when constraint-based interactions are transmitted to users by a unilateral coupler ( " C B " ) . Collision impulses computed using the proposed method are spread over UC  several simulation steps when necessary  63  4.20 Simulated forces and torques applied to users by the four channel controller when the virtual world is generated using the constraint-based ("ideal"), the impulse-augmented penalty-based ("IAPB L M " ) and the penalty-based ("PB") methods, and applied by the unilateral coupler when the virtual environment is generated using the constraint-based method ( " C B ) " . ColliUC  sion impulses computed using the proposed technique are spread over several simulation steps when necessary  64  4.21 Kinetic energy of the users' hand when constraint-based ("ideal"), impulse-augmented penaltybased ("IAPB L M " ) , and penalty-based ("PB") interactions are applied to users by the four channel controller, and when constraint-based interactions are applied by the unilateral coupler ("CB "). UC  Collision impulses computed using the proposed method are spread over several  simulation steps when necessary  65  4.22 Device trajectories when the peg is pushed into the corner by a controller force. The virtual environment is generated using the impulse-augmented penalty-based and the penalty approaches ("IAPB" - impulse-augmented penalty-based world, full collision impulses applied to the device; " I A P B L M " - impulse-augmented penalty-based world, limited impulses applied to the device; " P B " - penalty world). The "ideal" trajectory is generated using a constraint-based simulation.  67  L I S T OF  FIGURES  x  4.23 Environment wrenches applied to the device when the peg is pushed into the corner by a controlled force. The virtual environment is generated using the impulse-augmented penalty-based and the penalty approaches ("IAPB" - impulse-augmented penalty-based world, full collision impulses applied to the device; "IAPB L M " - impulse-augmented penalty-based world, limited impulses applied to the device; " P B " - penalty world). The "ideal" wrenches are generated using a constraint-based simulation  68  4.24 Kinetic energy of the device while it is pushed towards the corner by a controlled force. The virtual environment is generated using the impulse-augmented penalty-based and the penalty approaches ("IAPB" - impulse-augmented penalty-based world, full collision impulses applied to the device; " I A P B L M " - impulse-augmented penalty-based world, limited impulses applied to the device; " P B " - penalty world). The "ideal" wrenches are generated using a constraint-based simulation  69  5.1  Realistic haptic manipulation of virtual linkages in an impulse-augmented penalty virtual world.  71  5.2  1 D O F mechanical equivalent of the control architecture  72  5.3  Example manipulations restricted by the topology of the virtual linkage  76  5.4  Penalty wrench constraining users to the configuration manifold of the virtual linkage that they manipulate  5.5  77  Simulated manipulations of a planar virtual linkage. Initial linkage position is shown in black. Linkage positions during manipulation are shown in grey  78  5.6  Simulink diagram of the haptic manipulation of a planar virtual linkage  79  5.7  Simulated planar manipulation of a three links virtual tool held from the C O M of the distal link. 80  5.8  Simulated planar manipulation of a three links virtual tool held from the C O M of the second link. 80  5.9  Testbed virtual environment for controlled linkage manipulation from any user-selected link. . .  82  5.10 Haptic manipulation from the distal link of the virtual linkage shown in Figure 5.9. Users apply a constant wrench F  h  = (0.4N  ON 0 N m )  r  82  5.11 Haptic manipulation from the second link of the virtual linkage shown in Figure 5.9. Users apply a constant wrench F  h  6.1  = (0.4N  ON 0 N m )  T  Realistic haptic manipulation of virtual tools using the local model of interaction proposed in this thesis  6.2  6.4  85  Communication between the virtual environment (VE), the local model of interaction, and the haptic device  6.3  83  86  Example haptic interaction during which the computational delay of the virtual environment results in delayed updating of the virtual tool contacts in the local model of interaction  88  The influence of the geometry of the virtual environment on the interactions of the virtual tool.  88  LIST OF FIGURES  6.5  xi  Geometric information provided by typical simulation packages for each type of virtual tool contacts (Figures 6.5(a), 6.5(b), and 6.5(c)) and its representation in the local contact geometry (Figure 6.5(d))  90  6.6  Partial virtual environment geometry available in the proposed local model of interaction. . . .  91  6.7  Local proxy deformation due to violation of new constraints  92  6.8  Neighboring constraints included in the local model through the use of an e-active virtual tool, obtained by sweeping a sphere of radius e over the volume of the virtual tool  6.9  93  The e-active virtual tool eliminates the singularity in the constraint normal computation at a vertex-vertex contact between the virtual tool and the virtual environment by selecting a unique constraint normal direction  93  6.10 One dimensional peg-in-hole task used to illustrate the influence of local geometry on user's perception of locally cluttered virtual environments  94  6.11 Hand position for the simulated peg-in-hole task depicted in Figure 6.10(a)  95  6.12 Contact clustering in the virtual environment  96  6.13 Typical commercial simulation packages use polygonal representations (in a planar virtual world) for both rounded objects and polygonal objects  97  6.14 Experimental interaction used to investigate the passivity of the collision resolution approach proposed in this thesis  101  6.15 Passivity of proxy dynamics in an impulse-augmented penalty local model of interaction. The user operates the virtual linkage in Figure 3.3(b) from its distal link  102  6.16 Environment wrenches applied to users during user manipulation of the distal link of the articulated object in Figure 6.14(a) in an impulse-augmented penalty local model of interaction with perfectly plastic contacts, e = 0  103  6.17 Environment wrenches applied to users during user manipulation of the distal link of the articulated object in Figure 6.14(a) in an impulse-augmented penalty local model of interaction with perfectly elastic contacts, e = 1 6.18 Experiment used to investigate the haptic rendering of static contact  104 105  6.19 Device trajectories obtained when the virtual tool shown in Figure 6.18(a) is pushed with a constant force towards the bottom-most horizontal wall of the virtual environment, as depicted in Figure 6.18(b).  The intermediate representation ("IR"), the local model including active  geometry ( " L M " ) , and the local model including e-active geometry within 5mm from the virtual 0  tool ( " L M 5 " ) interface the haptic device to the virtual environment simulation  105  6.20 Experiment used to demonstrate the role of e-active geometry in improving user's perception of tight virtual clearances  107  L I S T OF  xii  FIGURES  6.21 Device trajectories obtained when the virtual tool shown in Figure 6.20(a) is rotated by 90°, inserted into the tight-fitting hole at the bottom of the virtual world, and shaken horizontally with a sinusoidally varying force, as depicted in Figure 6.20(b). The intermediate representation ("IR"), the local model including active geometry ("LMo"), and the local model including eactive geometry within 5mm from the virtual tool ( " L M 5 " ) interface the haptic device to the virtual environment simulation  107  6.22 Device trajectories obtained when the virtual tool is inserted in the tight-fitting hole at the bottom of the virtual world shown in Figure 6.20(a) and shaken horizontally with a sinusoidally varying force, as depicted in Figure 6.20(b). The local model including active geometry ("LMo") and the local model including predicted geometry within 5mm from the virtual tool ( " L M " ) 5  interface the haptic device to the virtual environment simulation 6.23 Experiment used to investigate the realism of the haptic rendering of dynamic contact  109 110  6.24 Device trajectories obtained when the virtual tool pushes the last link of the virtual linkage in Figure 6.23(a) with a constant force, as shown in Figure 6.23(b). The intermediate representation ("IR"), the local model including predicted geometry and static constraint ( " L M s c " ) ' and the local model including prdicted geometry and moving (kinematic) constraints ( " L M x c " ) interface  D.l  the haptic device to the virtual environment simulation  110  The twin pantograph planar haptic interface  134  D.2 The parameters of one pantograph linkage. The base joints are shown in black. The pantograph endpoint P is the joint between the two not grounded links e  135  D.3 The four channel architecture for 1DOF haptic interaction within virtual environments. Simulated blocks are shaded  137  Glossary P  damping coefficient  J € TZ  contact Jacobian  A6K  operational space inertia matrix  dxc  c  6 x 6  f € TZ  Cartesian space force  3  F G TZ  Cartesian space wrench (force and torque)  6  J * G TZ  weighted generalized inverse of J  J h G TZ  Jacobian matrix of a linkage computed at the user's hand  cxd  6xd  M Q €TZ ,M  €lZ  dxd  metrics on the linear vector spaces Q, V  6x6  v  n G TZ , t G TZ  Cartesian space directions (of the contact, of the relative sliding velocity)  ri( S 7Z  configuration space direction  p G 7Z  vector of magnitudes of frictionless collision impulses  3  3  d  C  C  q G 7Z , q G 7Z , q G TZ  configuration space position, velocity, accceleration  R G TZ  rotation from world to proxy coordinates  T G 7Z  configuration space torque  u G TZ  device control signal  d  d  d  3x3  d  6  x-h € TZ , x/j G 7Z , x/j G TZ  body position, velocity, acceleration at the user's hand  cu G TZ  angular velocity of a rigid body  T G TZ  Cartesian space torque  6  6  3  3  J  G TZ  cxd  c  fi, e  6  dynamically consistent inverse of the contact Jacobian J  c  coefficient of friction, coefficient of restitution  Glossary  xiv  £ € TZ , v € 7£ , a € TZ  Cartesian position, velocity, acceleration of a point  £j  i-th singular value of a matrix  3  3  D, G , n ,M y c  3  value of D , G , n , M y in the virtual environment at the moment of the c  simulation update a  stiction gradient  c  number of contacts of a contact group  Ct>, C  rigid body constraints (bilateral and unilateral)  d  dimension of the configuration space of a contact group  F, T, Q, V  linear vector spaces  KE  kinetic energy  s  separation distance between two rigid bodies  Z  one dimensional impedance  z  strain of the bond between two contacting bodies  zo  breakaway distance (maximum strain of the bond between two contacting bodies)  u  1Z ( A ) , N (A)  range space, null space of A  Acknowledgements I would like to acknowledge my thesis advisors, Professors T i m Salcudean and Elizabeth Croft, for their guidance and support throughout this work. They have allowed me the freedom to roam and have asked all the right questions. I feel very fortunate to have worked with them. If one day somebody will be as grateful for having worked with me as I am today for having worked with T i m and Elizabeth, then I have learned something from them and their efforts have not been wasted. I would also like to thank my colleagues in the Robotics and Control Laboratory for the times we shared and for the friendships we built. I hope our friendship will endure through years to come. M y special thanks go to Shahin, for his patience to discuss system dynamics with me. I am also greatful to my colleagues in the Industrial Automation Laboratory for our instructive and enjoyable laboratory meetings. Above all, I thank my family for supporting all I do. Everything would be worthless without you.  xv  Chapter 1  Introduction This thesis is concerned with providing realistic force, or haptic, feedback to users interacting with rigid multibody virtual environments. Lifelike forces give users a sense of "presence" in the computer-generated virtual environment; that is, they make users feel as if they are directly manipulating the rigid virtual objects. Realistic force interactions are important for physical skill training and for gaining an understanding of the virtual world that cannot be obtained using other human-computer communication paradigms. For example, astronauts may practice the operation of specialized tools under conditions similar to those during the space mission before the tools are manufactured; engineers may recognize limitations of their designs by manipulating virtual prototypes; students may become familiar with various physical phenomena by interacting with virtual objects. Despite significant potential benefits, commercial applications enabling force interactions with rigid multibody virtual worlds do not exist yet. Such interactions continue to be a focus of on-going research, including of this thesis. It is challenging to apply lifelike forces to users manipulating rigid virtual objects because the key components of a haptic interaction system impose conflicting requirements on the interaction (see Figure 1.1). Users require virtual forces at rates of the order of 500Hz or higher during kinesthetic interaction with rigid objects [28, 34,135] so that they do not perceive artifacts such as vibrations. Additionally, they exchange energy with the simulation. The device controller needs similar high "haptic rates" in order to create an adequate perception of rigidity. It also requires fixed force update rates and passive virtual environments to stably transfer energy between the user and the simulation [3]. Moreover, the device can display only limited stiffnesses to the user [40]. However, typical interactive physically-based simulations of rigid multibody virtual 1  2  environments either have rendering frequencies only of the order of tens of Hz [140,149], and can slow down further during complex interactions, or have geometry-dependent stiffness [42] that cannot be guaranteed for arbitrary interactions [59]. In other words, while the device requires realistic forces updated at guaranteed high rates and arising from a limited environment stiffness, existing simulations provide these forces with uneven computational delay or cannot guarantee the stiffness that they represent.  In prior haptics research, three  approaches have been developed to address the de-stabilizing effect of the simulation. User + device controller  Virtual environment simulation  - fixed force update rates (hundreds of Hz) - limited virtual stiffness  - variable update rates (tens of Hz on average) - geometry-dependent stiffness  Figure 1.1: The key components of a haptic interaction system and their conflicting requirements for stability and realism.  A first approach used control techniques to guarantee stable user interaction within virtual environments with arbitrary stiffness and/or computational delay (Figure 1.2). The most common method separates the simulation from the device through a "virtual coupler" [42]. The virtual coupler is a controller that limits the stiffness displayed to the user, eliminating the potentially destabilizing effect of the geometry-dependent environment stiffness.  To date, it has been used for rendering both point [128, 140] and rigid body [22]  interaction within rigid multibody virtual environments. The virtual coupler guarantees stable manipulation of passive nondelayed virtual worlds [3]. However, the restriction to passive nondelayed environments is quite limiting, since existing nondelayed rigid multibody simulations use penalty-based dynamic algorithms and fixed step forward integrators that are nonpassive [31]. Furthermore, the virtual coupler introduces perceptual artifacts (such as contacts that pull on the user) and does not allow users to perceive physical phenomena that rely on fast force transitions (such as collisions and stick-slip friction). Recently, a control technique was designed in [65] that monitors the energy produced by one or more simulation blocks and dissipates this energy when necessary. The technique guarantees stable user interaction within any virtual environment that produces  3  less energy than the device can dissipate, regardless of the computational delay of the simulation. However, the energy dissipation mechanism introduces perceptual artifacts (e.g., it allows contacts to pull on the user) that compromise the realism of the interaction.  Human user •  Haptic device high Virtual „ fixed or fixed & variable^ environment haptic rates coiuiol .ilgorithms rates simulation  «  Figure 1.2: Alleviation of the effect of the computational delay and the geometry-dependent environment stiffness via control.  A second approach improves the real time performance of the simulation via new collision detection methods that eliminate the collision detection bottleneck and significantly decrease the computational delay of the virtual environment. The new methods take advantage of hierarchical spatial bounding volumes [59,75,85], incremental distance computations [76,83], and multiresolution bounding hierarchies [117] to detect contacts and to approximate penetration distances at the high haptic speeds. Combined with penalty-based dynamic algorithms and fixed step forward integration routines that require little additional computation time, the new collision detection techniques allow the simulation to run synchronously with the haptic controller (see Figure 1.3). However, instability may arise if the user exerts too much force which causes a large amount of penetration [59], or in non-convex contact configurations [85].  Human user  r . . > Haptic device , high Virtual environment simulation & fixed fixed ( tt Dynamics haptic rates control algorithms haptic rates jfpenalty-based' dicci-riim  Figure 1.3: Elimination of the computational delay of the simulation via new collision detection algorithms.  A third approach alleviates the haptic rate demand on the virtual environment by decoupling the force control loop from the simulation (see Figure 1.4). Therefore, virtual environments of increased complexity, based on passive algorithms, can be rendered. In this approach, users feel realistic forces only if the setpoints of the control loop over one simulation step suitably approximate the forces acting on the virtual object manipulated by the user (the virtual tool). Hence, a local model of the interaction (i.e., a reduced simulation that runs at the fixed haptic rate) must be available in the force control loop. A first local model for point interaction with rigid virtual objects was proposed in [2]. It incorporates the geometry of the constraint  1.1 Objective  4  restricting the user's motion. This local model shifts the computational delay of the simulation from delay in computing the interaction forces to delay in updating the local geometry. The model was refined in [100,149] to ensure geometry continuity at the model updates. However, its extension to haptic interaction with a rigid multibody virtual world through a virtual tool led to instability during tightly constrained motions of the virtual tool, such as during peg-in-hole insertion [21]. Therefore, the local geometry was used to constrain only the translation of the device and a virtual coupler was used to constrain its rotation. Human user  Haptic device high fixed ^ & fixed haptic rates control algorithms haptic rates  l.ocnl modi/I ut' inrcHcrtuii  > Virtual environment ^ ' w simulation update^ rates .variable  >  Figure 1.4: Mitigation of the effect of the computational delay of the simulation via a local model of interaction.  None of the approaches described above allows realistic and guaranteed stable haptic manipulation of rigid virtual objects. Such interaction continues to be a challenging proposition. Both new control and new simulation methods need to be developed in order to increase the physical accuracy of the haptic rigid object manipulation while maintaining its stability.  1.1  Objective  The goal of this thesis is to increase the realism of the force interaction within rigid multibody virtual environments through modeling and simulation methods that enable efficient and transparent haptic rendering of rigid body motion with constraints. In particular, general-purpose simulation methods are proposed that model interactions convincingly, yet suitably for being applied to the user's hand by the device controller. Such simulation methods have not been available in prior research. A further objective of this work is to develop local approximations of the proposed simulation approaches that adequately represent the interaction and run at the high haptic speeds regardless of the complexity of the virtual environment.  1.2  Contributions  This work contributes to: (i) physically-based modeling of rigid multibody virtual environments; (ii) haptic manipulation of virtual linkages; and (iii) local models of interaction for rigid multibody virtual worlds.  5  1.2 Contributions  • Physically-based modeling of rigid multibody virtual environments Collisions are ubiquitous during physical manipulations of rigid objects. In haptics, collisions felt by the user enhance the perceived rigidity of the virtual contacts [90,125]. To increase the realism of the haptic interaction within rigid multibody virtual environments, this thesis develops a simulation approach that enables users to perceive collisions. The approach uses a contact model compatible with fixed step integrators and closed-form dynamics algorithms. Specifically, contacts are infinitely stiff upon contact and have limited stiffness during contact (see Figure 1.5). Furthermore, a new simultaneous collision resolution algorithm is employed upon contact. The algorithm is proven to never increase the kinetic energy of the virtual environment. Users feel collisions through impulsive forces. They feel contacts through penalty and friction forces.  (a)  (b)  (c)  (d)  (e)  Figure 1.5: The contact model the enables users to feel collisions: contacts have infinite stiffness upon contact and finite stiffness during contact. Users feel impulsive forces upon contact and penalty and friction forces during contact.  Although developed for haptics, the simulation approach proposed in this thesis can be used in any application that requires rigid multibody virtual environments with interactive performance. The approach provides a compromise between the efficiency of penalty-based simulations [62, 111] and the accuracy of constraint-based simulations [7,10,13,35,63,66,70,110,127,141.143]. • Haptic manipulation of virtual linkages In contrast to prior work that allows operation of spatial (planar) virtual linkages from links with at least 6 DOFs (3 DOFs) [114,128,140], this thesis develops methods that enable realistic rigid body manipulation of linkages from any user-selected point. In particular, users are enabled to manipulate  1.2 Contributions  6  links with motion limited by the linkage topology and to convincingly perceive motion restrictions due to other virtual objects and to the linkage itself. This is important in applications like training for laparoscopic procedures, for example, when the scope restricts the motion of the virtual tool at the point of entry. Motion constraints due to other virtual objects are imposed through contact (impulsive, penalty, and friction) forces transformed to the user's hand in a coordinate-invariant manner. The linkage topology is imposed through penalizing users when they violate the configuration manifold of the linkage (see Figure 1.6). Direction of topological constraint  Figure 1.6: Users perceive topological constraints through penalties applied to the device when their hand leaves the configuration manifold of the linkage (vh is the position and orientation of the user's hand on the virtual linkage, while h is their real position and orientation).  • Local models of interaction within rigid multibody virtual worlds Regardless of the efficiency of the algorithms employed, the computational performance of rigid multibody virtual environments depends at least linearly on the complexity of the virtual scene to be simulated (i.e., its contact geometry [12,59,76,83,84] and the dimension of its configuration space [48,91,127]). To enable haptic interaction within virtual worlds of arbitrary complexity, this thesis develops a local model of rigid body interaction that is simple enough to run at the haptic rates and general enough to interface a device to a virtual environment regardless of the algorithms used to simulate the virtual world. This model is the first local model that includes a dynamic proxy of the virtual tool manipulated by the user and active and predicted motion constraints imposed on the virtual tool by nearby objects (see Figure 1.7). Geometry  1.3 Organization of the thesis  7  prediction enables realistic manipulations of the virtual tool through tight clearances. The dynamic proxy encapsulates motion constraints imposed on users by articulated structures, allowing them to manipulate linkages. Furthemore, the haptic controller can coordinate both forces and motions between the haptic device and the proxy. Coordination of forces increases transparency and enables users to feel collisions, thus enhancing the perceived local contact rigidity. Position coordination allows stiffer objects to be manipulated.  Figure 1.7: The local model of rigid body interaction includes a dynamic proxy of the virtual tool and active and predicted motion constraints imposed on the virtual tool by nearby objects.  1.3  Organization of the thesis  The remainder of the thesis is arranged as follows: Chapter 2:  The state of the art in physically-based rigid multibody simulation and impact modeling is presented, focusing on methods with interactive performance. Simulation and control techniques employed for haptic interaction within rigid multibody virtual environments are also over viewed.  Chapter 3:  This chapter describes the planar haptic interaction system used to experimentally validate the simulation and the control methods proposed in the present thesis.  1.3 Organization of the thesis  Chapter 4:  8  A new approach to haptic rendering of rigid body contact is introduced that enables users to feel collisions and, thus, enhances the perceived rigidity of the v i r t u a l objects.  The  chapter describes the underlying contact model, as well as the simulation of impulsive and contact interactions. A new, passive, multiple collision resolution technique is proposed and modified to account for the practical limitations of haptic devices.  Simulations  and controlled experiments are used to validate the proposed rigid contact rendering approach. A friction model developed i n computational dynamics is, for the first time, used to haptically render dry friction.  T h e efficacy and the advantages of this m o d e l  are evaluated by comparing it to models priorly used i n haptics through simulations and experiments. C h a p t e r 5:  In this chapter, simulation and control methods are devised that enable realistic m a n i p u lation of v i r t u a l linkages. Linkage contacts are modeled as described i n C h a p t e r 4. T h e y are applied to the user's hand using a coordinate-invariant representation.  The inertia  and the topology of the v i r t u a l linkage are rendered to users through device control. T h e limitations of such control are overcome through a physically-motivated technique that convincingly restricts the m o t i o n of the device as required by the topology of the v i r t u a l linkage. Realistic and unrestricted operation of v i r t u a l linkages is demonstrated through simulated a n d experimental manipulations of a planar v i r t u a l linkage from various links. Chapter 6:  A local model of rigid body interaction is developed and used for haptic m a n i p u l a t i o n of rigid v i r t u a l objects.  T h e communication between the local model and the device,  and the local model and the v i r t u a l environment is described. Geometry and dynamics techniques ensuring model continuity at simulation updates are detailed. A numerically efficient implementation of the methods i n Chapters 4 and 5 is proposed.  T h e novel  local model of rigid body interaction is validated experimentally by using it to interface a planar haptic device to a virtual environment generated using a commercial s i m u l a t i o n package.  1.3 Organization of the thesis  Chapter 7:  This chapter summarizes the contributions and limitations of the work described in this thesis and discusses directions for future work.  9  Chapter 2  Physically-based rigid multibody simulations In haptic applications, both the computational speed and the physical accuracy of the simulation are essential for the user's sense of presence in the virtual environment. Physically-based rigid multibody simulations achieve interactive performance by balancing efficiency and accuracy in various ways. This chapter starts by surveying interactive simulations and the modeling methods they employ. The survey emphasizes the potential of existing techniques to achieve the requirements of haptic interaction, namely physical realism and guaranteed real time performance. The chapter ends by overviewing haptic manipulation of rigid multibody virtual environments, including simulation and interaction control techniques, local models of interaction, and impact and friction rendering.  2.1  Interactive simulations  Simulation algorithms implement a loop consisting of three main steps (see Figure 2.1): • collision detection: the body contacts are computed; • dynamic response: the equations of motion (EOMs) are solved in accordance with the new set of contacts; • time integration of the system state: the body positions and velocities at the next time step are computed. A l l three steps contribute both to the physical accuracy and to the efficiency of the simulation. For example, recent collision detection algorithms [75, 76, 83, 85,108] have significantly alleviated the collision detection 10  2.1 Interactive simulations  11  bottleneck and correspondingly improved the simulation performance.  However, only work addressing the  dynamic response is relevant to the scope of this thesis and, therefore, surveyed in this section. The focus is on the ability of existing dynamic response techniques to realistically represent rigid multibody interaction while achieving the real-time performance required by the haptic controller.  Collision detection  Time integration  1 )\ I1.1I111C-  rcsponsc  Figure 2.1: The simulation loop and the scope of the survey: prior work in dynamic response.  Dynamic response includes the formulation of the E O M s of rigid multibody systems starting from the fundamental principles of the Newtonian mechanics (principle of virtual displacements, Gauss' least action principle, etc.), and their solution. The E O M s of a multibody system are derived as a set of ordinary differential equations (ODEs) using either positions and velocities or positions and momenta as state variables. However, general purpose simulations use positions and velocities as state variables exclusively. In their most general form, the E O M s of a system of b rigid bodies with d DOFs with positions and velocities as state variables are a set of n > d second order ODEs [146]:  M(q)q + Q (q,q,i) = Q(q,q,t) c  In Equation (2.1), M G TZ  nxn  => M (q) q +  (q, q, t) + Q  c ni  (q, q, t) = Q (q, q, t) .  (2.1)  is the system mass matrix; q G TZ , q G TZ , q G TZ" are the system acceleration, n  n  velocity, and position vectors, respectively; Q G TZ is a vector that includes external forces as well as Coriolis, n  centrifugal, and gravitational effects; Q = Q? + c  are  G TZ is the vector of constraint forces, whose components n  G TZ , the vector of ideal constraint forces (they do no work and are along the constraint directions), n  and Q£j G TZ", the vector of nonideal constraint forces (they do work and are in the constraint plane). Note that Equation (2.1) does not represent the E O M s as they are directly derived from the principles of mechanics; rather, it represents the E O M s rearranged for computational purposes. Unknown forces and accelerations are grouped on the left hand side of this equation, while known forces are on its right hand side. If the b bodies have c contacts, the n-dimensional Equation (2.1) has n + 2c unknowns: the system acceleration q, and the constraint forces  and Q^. Additional equations are needed to find the dynamic response of the system,  12  2.1 Interactive simulations  i.e., to compute the system acceleration q to pass it to the time integration routine. The additional equations are given by the constraints and by the laws of dry friction. The constraints enforce the rigidity assumption by requiring the bodies not to overlap. Bilateral constraints represent permanent connections between bodies, or joints, whereas unilateral constraints represent temporary connections, or contacts. If the system has j joints and c contacts, then there are j bilateral constraints:  G\(q,q,t) = 0  Vt = l , - - - j ,  (2.2)  and c unilateral ones: C (q,£)>0 U i  Vi = l , - - - , c .  (2.3)  Additionally, the Coulomb law of friction requires the nonideal and the ideal constraint forces to be related through: \\Qnu\\<Vi\\Q u\\ c  where the index i selects the i-th vector component and  Vi =  l,---J+C,  (2.4)  is the coefficient of dry friction at the i-th contact.  The constraints in Equation (2.2) and the E O M s form a set of second order differential algebraic equations (DAEs). Integration of DAEs is unsuited for real-time performance. Therefore, in interactive simulations, these D A E s are transformed into ODEs by appending the constraints to the dynamics at the acceleration level (i.e., by differentiating Equation (2.2) twice with respect to time):  Ajxn(q,q,i)q = b j i ( q , q , t ) . X  (2.5)  In Equation (2.5), A is the Jacobian of the constraints (the constraint matrix) and b groups all other terms obtained by differentiation. Dynamic response requires the systems acceleration q to be computed. However, Equations (2.1) and (2.5) form a system of n + j equations with n + 2j + 2c unknowns, Qf, and Q ^ . The additional equations necessary for computing the dynamic response depend on the employed model of rigid body contact, as described in Section 2.1.2. The numerical efficiency and the physical accuracy of the dynamic response depend both on the coordinates used to describe the system state and on the model of rigid body contact, as described in Sections 2.1.1 and 2.1.2.  13  2.1 Interactive simulations  2.1.1  Coordinates for rigid multibody dynamics  Rigid multibody simulations describe the system state using Cartesian, generalized (dependent), or independent coordinates (see Figure 2.2). The coordinate system used depends on whether bilateral constraints are distingushed from unilateral ones.  (a) C a r t e s i a n c o o r d i n a t e s .  (b) G e n e r a l i z e d c o o r d i n a t e s .  (c) I n d e p e n d e n t c o o r d i n a t e s .  Figure 2.2: Coordinates used to describe the system state in rigid multibody simulations.  Cartesian coordinates (Figure 2.2(a)) do not differentiate bilateral constraints from unilateral ones. They incorporate no topological information into the EOMs. The dynamic model of each body is independent of all others and constraint equations are added or deleted from the constraint matrix A as needed [56]. Cartesian coordinates are best suited for general purpose simulators. They are used in many implementations, developed in computational dynamics [56,57,66,116,153,154]), in computer graphics [5,12,16,32,37,62,104,109,110, 122,134], and in haptics [34,85,117]. However, models in Cartesian coordinates result in the largest number of dynamic equations (three per body for planar interaction and six for spatial interaction) and are prone to constraint violation. The numerical drift is addressed by constraint stabilization [5,19, 37] and by projection methods [82,98], which integrate only d independent coordinates and use the constraint equations to compute the remaining n — d position variables. Generalized coordinates (Figure 2.2(b)), e.g., relative coordinates [48] and natural coordinates [51], incorporate the bilateral (joint) constraints into the EOMs and explicitly represent only the unilateral (contact) constraints. Hence, the number of E O M s to be solved is smaller, equal to the number of DOFs of the system. Moreover, the system evolves on its configuration manifold and the bilateral constraints are automatically satisfied. Efficient serial algorithms are proposed in [25,48,73,74,91,107] for loop-free topologies and in [72] for topologies with internal loops, while efficient parallel algorithms for linkages with or without loops are presented  2.1 Interactive simulations  14  in [49,50]. Collision resolution and contact force determination methods for tree-like systems are described in [127], where the increased numerical efficiency of the coordinate representation allows the simulation of a 26 DOFs tree-like robot at the haptic rates. The significant computational advantage of generalized coordinates has led to the development of several simulators with interactive performance [13,51,128,140,141], some of which have also been used for haptic interaction [128,140]. Independent coordinates (Figure 2.2(c)) are the minimum number of coordinates required for describing a multibody system, i.e. n = d. Minimal representations result in a minimum number of E O M s and no numerical drift, since all constraints are incorporated in the choice of independent coordinates. However, the set of independent coordinates changes as the constraint configuration changes. This set cannot be derived numerically at present. Hence, in interactive applications, the efficiency gained by embedding the constraints in the E O M s may be offset by the symbolic calculations needed to derive a new system model at configuration changes. A n application-specific simulator using independent coordinates has been demonstrated in [52,54] for haptic rendering of piano playing. In this simulator, all constraint configurations are known prior to run time. Hence, dynamic models corresponding to these configurations are derived off-line and then combined on-line using a finite state machine. Recently, the design of a general purpose haptic simulator using independent coordinates has been outlined in [53]. This design proposes to achieve interactivity by decoupling the symbolic computation of the E O M s from their numerical solution through running the symbolic model formulation in a thread outside the main simulation loop. Further work is required to demonstrate the interactive performance of this design. In addition to the choice of coordinates, the enforcement of the rigidity and of the dry friction assumptions also influence the computational efficiency and the physical accuracy of rigid multibody simulations. This influence is discussed in the following section by overviewing dynamic models of rigid multibody systems available in prior work.  2.1.2  Dynamic models of rigid multibody systems  Rigid body dynamics have been a topic of active investigation for more than three hundred years. This long-lived interest has arisen from the mathematical difficulties accompanying two widely used modeling simplifications: the rigid body and the Coulomb dry friction assumptions. These simplifications bring about issues related: (i) to  15  2.1 Interactive simulations  the existence and the uniqueness of solutions of the E O M s of rigid multibody systems with unilateral constraints and dry friction [7,47,92,94,102,118]; and (ii) to the algorithms approximating these solutions [6,7,143,144]. Physically-based rigid multibody simulations differ in how they enforce these simplifications and, therefore, in how they use the dynamics given by Equation (2.1). Under the rigid body assumption, Newtonian mechanics distinguishes three contact states: free motion, collision, and contact. Under the Coulomb friction assumption, contact can further be divided into sliding and sticking contact. Collisions (or impacts) are phenomena with much faster time constants then contacts. They are often modeled as instantaneous events which result in velocity discontinuities. Contacts are modeled as finite duration events which result in acceleration discontinuities. Based on the contact states they model, rigid multibody simulation methods can be classified in: (i) penalty-based; (ii) constraint-based; and (iii) impulse-based approaches. 2.1.2.1  Penalty-based simulations  Penalty-based simulations model only free motion and contact (Figure 2.3). They enforce rigidity only approximately, by introducing restoring forces when objects interpenetrate.  Hence, they incorporate the penalized  constraints into the EOMs. Since they only loosely enforce constraints, penalty methods sacrifice accuracy for robustness (singular configurations pose no problem) and a reduced set of equations (reaction forces are eliminated from Equation (2.1), as they are computed from local contact geometric information). Thus, penalty methods are compatible with the fixed time step integrators required in haptics because they do not need to exactly determine when new contacts arise.  Figure 2.3: The contact states in a penalty-based rigid multibody simulation.  A main limitation of penalty simulations is that large penalties are necessary for acceptable accuracy. This induces numerical ill-conditioning during integration and requires the penalties to be adjusted for different simulation conditions. However, highly realistic computer animations of rigid bodies were demonstrated in [109, 121] and a general purpose simulation package was developed in [57,58]. A more important limitation for haptics is the fact that the environment impedance depends on the contact  16  2.1 Interactive simulations  geometry and cannot be bounded for interactions within arbitrary virtual worlds (i.e., for arbitrary contact geometries).  This is illustrated by the example in Figure 2.4. In this example, the contact between the  rectangular object and the environment is represented through c = 8 point contacts and the environment stiffness is represented as: (2.6) where k „t ct co  a  is the contact stiffness.  Regardless of the chosen contact stiffness, contact geometries may  exist that induce instability. Instability may arise because the environment stiffness is outside the region of conditional stability corresponding to the given time step of the forward integration routine used in the haptic simulation. As discussed in [40], instability may also arise due to an environment stiffness larger than the Z-width of the device (the impedance range that the device can render stably) [29].  Figure 2.4: A rectangular object interacting with a virtual wall. Each surface contact is represented through two point contacts.  Despite the difficulties ensuing from the conditional stability of penalty methods, their low computational requirements and ease of implementation are very attractive for haptics. Several implementations of such methods exist [59,75,76,83,85,103,117]. 2.1.2.2  Constraint-based simulations  Constraint-based simulations model all contact states and enforce exact rigidity by computing either collision impulses or contact forces to prevent body interpenetration (Figure 2.5). Interactive constraint-based simulations are more difficult to implement because they need to detect contact state changes at run time. This means that, given the set of geometric contacts reported by the collision detection algorithm, constraint-based simulations must determine the set of "active" contacts, i.e., the set of contacts at which objects interact. Such detection is possible if the system dynamics is derived based on the complementarity rule (Figure 2.6), also known as the corner law or Signorini's conditions [120]. In the direction normal to the contact, the rule  17  2.1 Interactive simulations  Figure 2.5: The contact states in a constraint-based rigid multibody simulation.  states that either the relative acceleration a  is zero and the contact force / „ is nonnegative , or the relative 1  r > n  acceleration is nonnegative and the contact force is zero. In the contact plane, it states that either the relative tangential velocity v  T>n  is zero and the friction force / / is in the friction cone (i.e., / i /  n  — / / > 0, where / i is  the coefficient of dry friction), or the relative tangential velocity is nonnegative and the friction force is on the boundary of the friction cone (i.e., jif — ff = 0). n  a..  (a) The normal contact direction.  (b) The contact plane.  Figure 2.6: The corner law.  Two types of complementarity-based dynamics formulations exist: acceleration-force methods [93,94] and velocity-impulse methods [112]. Correspondingly, event-based and time-stepping rigid multibody simulations were developed [27]. Event-based simulations [11,120,122,127] integrate the system dynamics between collision events, resolve occurring collisions, and reset the integrator before continuing to advance the system state in time. Such algorithms are not directly useful in haptics because they require variable step size integrators; i.e., there is no guaranteed completion time. Time-stepping simulations [6-8,37,104,142,143] solve a time discretization of the system dynamics that includes the impact rules and the perfectly rigid constraints enforced at the velocity level. These algorithms are compatible with the fixed step haptic integrators, but have no guaranteed completion time. This is because T h e contact force is positive if it pushes the bodies apart. The relative acceleration and the relative velocity are positive if they are such that the bodies separate. 1  18  2.1 Interactive simulations  they can avoid constraint drift only through accurate collision detection and the reformulation of the system dynamics at each contact configuration. Since there exists no a priori bound on the number of collisions that can occur during one simulation step, the computation time is unpredictable. Time-stepping algorithms with guaranteed completion time include suitable constraint stabilization techniques [6,37]. In [37], constraint stabilization is formulated as a complementarity problem solved at the end of the simulation step. Besides the computational price paid by solving two linear complementarity problems (LCP) per step, the method is not shown to ensure the passivity of the simulated system.  In [6], constraint stabilization is achieved  by augmenting the complementarity-based system dynamics with a term that depends on the geometrical constraint infeasibility. Hence, no additional computational costs are incurred. However, the method is proven to ensure the passivity of the simulated system only for perfectly plastic collisions between smooth convex objects.  Further work is needed to extend the method to elastic and partially plastic collisions between  polyhedral objects, which are only locally smooth and may not be convex. 2.1.2.3  Impulse-based simulations  Impulse-based simulations model only free motion and collisions (Figure 2.7). They enforce rigidity exactly, by representing any interaction as staggered collisions and computing impulses that prevent body interpenetration. Contact is represented as a series of microcollisions. Friction is dealt with during microcollisions. The impulsebased technique was proposed in [62] and was developed into a general purpose simulation package in [110]. Recently, the physical accuracy of the method was traded off for efficiency by adding heuristics that enable interactive performance for simulations as large as 1000 cubes in an "hourglass" [134]. These heuristics are visually acceptable, since the number of virtual objects is too large for humans to be able to predict the natural motion of individual bodies. Impulse-based simulations are suitable for simulating systems with frequently occurring impacts (such as a part feeder, for example), but inefficient at simulating systems with joints and persisting contacts.  Figure 2.7: The contact states in an impulse-based rigid multibody simulation.  2.2 Impact simulation  19  A haptic implementation of the impulse-based method was proposed in [34], but the results were not encouraging. In particular, the perception of contact was not convincing and users could not feel friction, regardless of the fact that it was built into the collision model [34]. There are two possible causes for this. First, the representation of contact through series of microcollisions could be satisfactory for visual feedback, but inappropriate for kinesthetic feedback. Second, it might be challenging to reconcile fixed step integration with exact enforcement of rigidity in a physically meaningful manner. For example, invalid impact states can result during a fixed step impulse-based simulation. These are states where two objects are in contact, but they are separating from each other and collision impulses cannot be computed. Because only collisions are modeled, contact can be maintained only by assigning a "computational" contact state to the objects. Assignments of a "computational" contact state other than those used in [34] might be needed to enable users to feel friction. This overview of physically-based simulation techniques with interactive performance reveals that no type of coordinates or simulation method is ideally suited for generating the rigid multibody virtual environment in haptic applications. For example, further work is required for guaranteed stable interaction with penalty-based simulations and constraint-based rigid multibody haptic simulations have yet to be demonstrated. Additional difficulties in rendering rigid contact arise from the fact that the haptic controller has limited capabilities (limited bandwidth and motor power, for example). Therefore, it cannot prevent user penetration into the constraints.  Nevertheless, psychophysical experiments [90,125] have demonstrated improved perception of  contact rigidity if large forces are applied to users upon contact. The present research associates such large forces with collisions and develops a collision resolution method suitable for haptic interaction. To put this method in perspective, relevant work in impact simulation is overviewed in the following section. As before, the emphasis is on the ability of prior models to realistically model collisions while achieving real-time performance.  2.2  Impact simulation  Impact simulation is important in a large number of applications, including analysis of vibrations in machines [120], simulation of robotic pushing [101] and orienting [55] operations, virtual prototyping [1], and interactive computer graphics [134]. At the same time, it is a challenging endeavor, mostly because the many physical phenomena occuring during collisions (restitution, tangential compliance, jamb, vibrations) are difficult to capture into simple models that apply to both single and multiple collisions. Therefore, many collision  20  2.2 Impact simulation  (or impact or restitution) laws (or rules) were proposed that aim to realistically predict the collision outcome. They fall into three categories: compliant, differential, and algebraic collision laws [27]. Compliant [57, 86, 87,106] and differential [15,18, 23, 24, 78, 99,115,145] impact laws treat collisions as phenomena with a finite duration, derive the corresponding equations of motion (under the assumption of constant body positions), and integrate them in order to compute post-impact body velocities. Such rules are adequate for resolving simultaneous collisions with friction [86,99], i.e., they predict collision outcomes that agree with basic constraints like nonnegative energy dissipation and the Coulomb dry friction inequality at the impulse level. However, compliant and differential restitution laws require much smaller integration steps during collisions. Hence, they are not suitable for simulations with interactive performance. Algebraic collision laws [4,7,26,35,63,120,137,142,143,151] ignore the mechanics of the interaction and predict the collision outcome using various generalizations of Newton's and Poisson's rules. Newton's law relates the pre-collision (w ,no) and post-collision (v ) relative normal contact velocities through the coefficient of 2  r  Ttn  restitution e: (2.7) Poisson's rule decomposes the collision into compression and restitution and relates the normal impulses during restitution p ^ and compression p n  r  n c  through the coefficient of restitution e:  Pn,r —  €Pn,c-  (2.8)  The coefficient of restitution e € [0,1] describes the nature of the collision, with e = 1 corresponding to perfectly elastic collisions (zero energy loss), e = 0 corresponding to perfectly plastic collisions (maximum energy loss), and intermediate values corresponding to plastic collisions (some energy loss). Friction is incorporated into collisions through Coulomb's dry friction law, by requiring that:  Pt < HP,  (2-9)  In this work, "normal" indicates the component of a vector along the direction normal to the contact plane, while "tangential" indicates the component of the vector in the contact plane. The normal direction and the tangent plane are well defined almost everywhere. 2  21  2.2 Impact simulation  if Newton's restitution hypothesis is employed, or:  Pt,r Pt,c  (2.10) <  A»Pn,c  if Poisson's rule describes the collision. In Equations (2.9) and (2.10), p is the coefficient of friction, and the subscripts n and t indicate vector components along the contact normal and in the contact (tangent) plane, respectively. Multi parameter algebraic collision rules were also proposed in prior work [35,120]. These laws model phenomena like tangential compliance in addition to restitution and friction. They are employed primarily when the need for accuracy outweights the need for performance. Algebraic collision laws are suitable for simulations with interactive performance because they do not require integration during collisions. However, their direct application to multiple collisions may result in energetically inconsistent results (i.e., an increase in the kinetic energy of the colliding bodies or dissipation of energy for perfectly elastic collisions), as pointed out in [26,27,78,145]. Some event-driven constraint-based simulations avoid such difficulties by sequencing collisions according to various heuristics [35,134]. However, collision sequencing is not suitable for haptics because it is not guaranteed to terminate in a finite number of steps [35]. Techniques for resolving multiple collisions simultaneously are needed. Such techniques were developed by formulating the dynamic response based on the complementarity rule and combining it with the time integration in time-stepping simulations. Time-stepping algorithms were developed based both on Newton's [112] and on Poisson's [6-8,37,104,142,143] restitution hypotheses. However, most existing time-stepping algorithms resolve the simultaneous collisions by iteration [7,8,104,112,142,143], thus having no guaranteed completion time. Time-stepping methods with guaranteed completion time [6,37] have been proven not to increase the kinetic energy of the system only for perfectly plastic collisions between smooth convex objects [6]. The extension of these algorithms to arbitrary coefficients of restitution and polyhedral objects is an open research question. This brief survey of restitution laws employed for impulse simulation indicates that only algebraic rules are suitable for algorithms with real time performance. Moreover, these rules must carefully be extended to simultaneous collisions to ensure that the kinetic energy of the system does not increase during collisions.  2.3 Haptic interaction within rigid multibody virtual environments  2.3  22  Haptic interaction within rigid multibody virtual environments  To enable haptic interaction with rigid multibody virtual worlds, prior haptics research either developed specialpurpose dynamic response algorithms or employed methods already proposed in computational dynamics, robotics, and graphics. Specific force computation techniques were developed: (i) for point interaction with static collections of polyhedra [129,150,155]; (ii) for point interaction with static objects described by implicit functions [97,133]; (iii) for line interaction with static polyhedral environments [17]; and (iv) for rigid body interaction with static environments represented as voxel maps [103,123]. Realistic rigid body interactions within dynamic multibody environments were computed using the general-purpose simulation methods described in Section 2.1.2. Most special-purpose techniques developed in haptics are penalty-based methods that achieve improved performance by taking advantage of the specific representation of the virtual environment [97,103,123,133] or the virtual tool [17] to speed up collision detection. The "god-object" [150,155] and the "virtual proxy" [129] techniques are constraint-based algorithms designed to overcome difficulties inherent in fixed step penalty-based simulations. They represent the user's finger through a point (sphere) that follows the user's motion as close as allowed by the objects in the virtual environment. These methods eliminate disadvantages inherent in fixed step penalty-based techniques. Namely, they prevent users from moving through thin objects in one simulation step and uniquely define a penetration distance at environment corners.  In turn, users perceive constant  environment stiffness throughout the interaction (while they perceive higher stiffness at concave corners and lower stiffness at convex ones in penalty-based simulations). However, the extension of the god-object and the virtual proxy methods to rigid body interaction is not straightforward [126]. Moreover, the perceived rigidity of virtual worlds simulated using these methods is similar to that of penalty-based environments. This is because god-object and virtual proxy simulations do not compute constraint-based forces and impulses and users feel only penalty-like forces applied to them by the interaction controller. Haptic rendering of rigid body interaction with dynamic multibody virtual environments was enabled in prior work by employing all general-purpose simulation methods described in Section 2.1.2. The impulsebased technique did not adequately render contact and friction in an initial implementation [34] and was not pursued further.  Penalty-based simulations are currently used in conjunction with advanced collision  detection algorithms. These algorithms achieve real time performance through combining local and global  2.3 Haptic interaction within rigid multibody virtual environments  23  search techniques [59,75,76,83,85] and eliminate artifacts due to discrete object representation through mesh decimation and mesh filtering methods [117]. Event-driven [22,128,149] and time-stepping [139] constraintbased simulations were used in haptics by bridging the gap between the virtual environment performance and the required haptic rate either through control or through a local model of interaction.  2.3.1  Haptic controllers  Control approaches aim to eliminate the instability caused by the (nonuniform) computational delay of the virtual environment through added virtual damping. They curb the energy flow from the virtual environment by using a "virtual coupler" [3,42] (Figure 2.8). The virtual coupler is a stiffness and damping connection between the haptic interface and the virtual environment that limits the maximum impedance displayed by the device. It guarantees stable interaction with passive, nondelayed virtual environments [3]. Although the passivity and the nondelayed assumptions are quite restrictive, the virtual coupler is often used for enabling interaction with simulations with variable delay [22,128,138,139]. Nonlinear damping was proposed in [65] to passivate the discrete implementation of virtual contacts and design guidelines were derived in [105] for adding virtual damping that guarantees stable interaction with both delayed and nondelayed nonlinear mass-springdamper (i.e., penalty-based) environments.  Human user  high fixed haptic rates  (  Haptic device &  Virtual coupler ^ fixed or f high \ fixed Virtual variable ^ P —i haptic rates environment rates 1 simulation  <  -  Figure 2.8: The virtual coupler.  One problem with adding virtual damping is that it degrades the transparency of the interaction by filtering out phenomena that rely on fast force transitions, such as collisions and stick-slip friction. In addition, the virtual coupler is a conservative design. The damping must be chosen such that stable interaction is guaranteed for the worst case interaction. A less conservative design was proposed in [95]. Specifically, for virtual environments with dynamics represented through several nonlinear functions, a virtual coupler was used for each function and stable switching between the various virtual couplers was devised based on passivity of hybrid systems arguments. However, stability is achieved by delaying the switches compared to the user commands, which, again, degrades transparency.  2.3 Haptic interaction within rigid multibody virtual environments  2.3.2  24  Local models of interaction  The local model of interaction approach decouples the force computation loop from the virtual world simulation. The decoupling alleviates the haptic rate demand on the simulation. Therefore, virtual environments of increased complexity can be rendered and multiple networked users may be allowed to manipulate the same virtual environment. However, realistic haptic feedback can be applied to users only if the setpoints of the control loop are suitable approximations of the forces acting on the virtual tool (i.e., the virtual object manipulated by the user). The techniques proposed for approximating the virtual interactions within rigid multibody virtual environments are presented in this section. The decoupling of the force computation from the virtual environment was first proposed in [79], in a haptic extension of the X-Windows graphical user interface (GUI). A haptic model of the GUI, running on a microcontroller, was used to compute forces at a fixed control rate of 1kHz. The model included constraints corresponding to window borders, pull-down menus, cursor forces towards icons, etc., and was updated asynchronously by the X-Windows host according to the status of the GUI. A local model was not needed due to the relative simplicity of the virtual environment (few constraints, typically horizontal and vertical). The earliest local model of force computation used in the control loop was proposed in [2] for point interaction within virtual environments. It included the position of the active constraint and its outward normal. Hence, it was a geometrical local model of the virtual environment. The local geometry shifted the computational delay of the simulation to delay in updating the local geometry. Therefore, larger contact stiffnesses were achievable, but force discontinuities arose at model updates due to geometry discontinuities. The force discontinuities introduced perceptual artifacts or destabilized the interaction. Improvements to this geometric local model aim to preserve stability and eliminate the perceptual artifacts by applying continuous forces to users. In [149], kinematic constraints were proposed for diminishing force discontinuities due to point contact with dynamic objects. The kinematic constraints are geometric constraints that move locally with constant speed, equal to their speed in the virtual environment at the moment of the update. Hence, the local model uses a first order predictor to compute the constraint position between simulation updates.  Rather than augmenting the local model, a force smoothing scheme based on linear  interpolation between old and new constraints was designed in [100]. The most advanced local model of point interaction was presented in [119]. This model includes all virtual environment geometry that might be reached  2.3 Haptic interaction within rigid multibody virtual environments  25  by the user until the next simulation update. Predicted geometry is selected based on the velocity of the user. This model was shown to ensure the continuity of the local geometry for user speeds of 0.2m/s in an environment with an average edge length of 7.5mm, typically running at 20Hz. One disadvantage of this model is that it is computationally expensive, at times barely finishing in one simulation step [119]. A local model of rigid body interaction, called the intermediate representation, was proposed in [21]. The intermediate representation includes all active constraints, enabling the user to feel crisp, stiff contacts. However, instability arose during tightly constrained motions of the virtual tool, such as during peg-in-hole insertion. Therefore, the intermediate representation was used to constrain only the translation of the haptic interface and a virtual coupler was employed to constrain its rotation. Local models were also proposed for haptic interaction with deformable objects [9,14,96]. They aim to accurately render variations in model stiffness to users. Hence, these models are not suitable for rigid body interaction, which requires users to feel constant environment stiffness, but varying inertia and restrictions imposed on their motion by other virtual objects. Increased transparency can be achieved within local models of interaction if continuity is ensured at updates. This is because users can feel stiffer contacts, impacts, and dry friction. A synopsis of the techniques employed for rendering impacts and friction to users is presented in the following sections.  2.3.3  Haptic rendering of impact  A key factor in enhancing the realism of the haptic manipulation of rigid multibody virtual worlds is the perceived contact rigidity. As demonstrated by psychophysical studies [90,125], the perceived rigidity of virtual contacts can be improved through applying large forces to users when new contacts arise. However, little haptics work exists that enhances the perception of rigidity through applying large forces to users when new contacts arise. The earliest large forces rendered to users are the "braking pulses" [131]. The braking pulses arise from a virtual wall model with high initial contact damping. They are designed to dissipate the entire kinetic energy of the haptic device during point interaction within virtual environments. Hence, the braking pulses can also be interpreted as impulsive forces arising due to a perfectly plastic collision. They are shown to improve both the stability and the perceived rigidity of the contact.  2.3 Haptic interaction within rigid multibody virtual environments  26  Large forces applied on users upon contact during rigid body manipulation of virtual worlds were physically motivated either as "pre-contact braking forces" [103] or as impulsive forces [43,44]. The pre-contact braking forces arise from the interaction with a viscous layer surrounding the virtual objects and dissipate the translational kinetic energy of the virtual tool along the direction of approaching contact. The impulsive forces arise from a planar model of rough impact employing on Poisson's restitution rule. Simultaneous collisions satisfy this model in a least squares sense and their passivity is not proven. Recently, open loop force pulses [71] were used to increase the perceived rigidity of surfaces. In addition to canceling user's momentum upon point contact with a rigid environment, the open loop pulses generate high frequency vibrotactile content that approximates the frequency content of contact with real surfaces.  2.3.4  Haptic rendering of friction  Since friction is ever-present in the physical world, realistic haptic manipulation of virtual environments requires it to be suitably simulated and applied to users. Several friction rendering techniques exist in haptics, all of them approximating Coulomb (dry) friction. Of course, Coulomb friction is only a macroscopic approximation of microscopic phenomena. Nonetheless, user studies [124] have shown that haptic rendering of dry friction affects task performance in a manner similar to that of real friction. Haptic rendering of friction is based on the simulation models proposed by Dalh [45], Karnopp [77], and Haessig and Friedland [61]. Various modifications of the Karnopp model were employed in [113,124,131,155]. For example, viscous friction replaced kinetic friction in [131], the sticking region was designed based on the characteristics of the human finger pad in [113], and both kinetic and viscous friction were used in [124]. The bristle model proposed by Haessig and Friedland was used to simulate both friction and adhesion in [36]. Lastly, a modification of the Dahl model was proposed in [67] that eliminates the spurious position drift of the original model. Additionally, the modified model depends only on position, which makes it robust to noise and suitable for implementation in event-based simulations with non-uniform sampling. A n implementation of the modified Dahl model for haptic interaction with a virtual environment through a virtual tool was reported in [85].  2.3.5  Summary  Rigid multibody modeling and simulation techniques with interactive performance have been developed in computational mechanics, robotics, and graphics research.  The potential of these techniques to meet the  2.3 Haptic interaction within rigid multibody virtual environments  27  specific demands of haptic manipulation of virtual environments (guaranteed real time performance for stable and transparent interaction, and physical realism for meaningful force feedback) has been considered. The investigation has shown that constraint-based and impulse-based techniques represent the physical interaction more accurately than penalty-based techniques, because they enforce exact contact rigidity. However, perfectly rigid contacts are compatible with the fixed step haptic integrators only if multiple collisions are resolved simultaneously and a constraint stabilization mechanism is available. In existing research, only constraint-based simulations of virtual environments with perfectly plastic contacts have been developed that are compatible with fixed step integrators. Haptic multibody simulation techniques, approaches to guaranteeing the force computation rates in arbitrarily complex virtual environments, and methods for increasing the realism of the interaction through impact and friction rendering have also been overviewed. Regardless of how the virtual world is generated, haptic force update rates are achieved in prior research by connecting the simulation to the device through a virtual coupler, i.e., by applying penalty-like forces to users through control. Decoupling of the simulation from the control loop through a local model of interaction has been demonstrated only for point interaction, while impact-like forces have been applied only to users interacting with planar virtual environments. A novel haptic simulation approach that allows users to feel impacts during their rigid body interaction within spatial virtual environments is presented in the following chapter. Its implementation in a local model of rigid body interaction is detailed in Chapter 6.  Chapter 3  E x p e r i m e n t a l setup This chapter presents the haptic interaction system that will be used throughout this thesis to validate newly proposed simulation and control techniques that enable efficient and transparent haptic rendering of rigid body motion with constraints. All experiments described in later chapters use the planar haptic simulation system depicted in Figure 3.1. This sytem comprises: a haptic device together with its input/output (I/O) interfaces and power amplifiers [136]: a haptics (real time) processor; and a graphics (host) processor with a graphical display. The haptic device is a 3 D O F interface that has a workspace larger than 100mm x 100mm and allows unlimited rotation. Hence, it enables planar rigid body haptic manipulations. The haptics processor is a 700MHz Pentium III personal computer running V x W o r k s ™ , while the graphics processor is a 2.4GHz Pentium IV personal computer running Windows2000™. The two processors communicate via a U D P socket. As depicted in Figure 3.2, the communication between the virtual environment simulation and the force control loop is implemented using two software architectures.  Current  DAC Haptics processor Decoder Ethernet  Haptic device  Graphical display  Figure 3.1: The planar haptic simulation system used for experiments throughout this thesis.  28  29  Synchronous communication is implemented by running both the virtual environment and the device control on the haptics processor (see Figure 3.2(a)), at a frequency of 512Hz. In this architecture, the graphics processor only displays the virtual environment, at a frequency of approximately 30Hz. The graphics processor acts as the server and polls the socket for new data asynchronously, each time it has finished displaying the old data. The haptics processor serves as the client and updates the socket at the end of each control step. Asynchronous communication is implemented by running the force control loop and the local model of interaction proposed in Chapter 6 on the haptics processor at 512Hz (see Figure 3.2(b)). In this architecture, the graphics processor both generates and displays the virtual environment, at frequencies varying between 30Hz and 60Hz. The haptics processor acts as the server. Hence, the local model polls the socket for new data at the beginning of each control step. When new data are available, the local model updates its state and acknowledges the receipt of the packet by sending back the proxy state. The virtual environment sends packets asynchronously, each time it has completed a simulation step. V o r t e x ™ , a physics based engine developed by CMLabs Simulations Inc. , is used to generate the virtual environment in the asynchronous communication 1  architecture. The force control loop comprises two components: the haptic interaction controller and the device controller. The haptic interaction controller [136] implements a four channel teleoperation architecture [89]. In the synchronous communication paradigm, the four channel controller transmits wrenches (i.e., forces and torques) between the device and the virtual object manipulated by the user (hereafter called virtual tool), and coordinates positions between them. In the asynchronous communication paradigm, the four channel controller transmits wrenches between the device and a proxy of the virtual tool in the local model of interaction, and coordinates positions between the haptic interface and the proxy. The impedance device controller [136] shapes the dynamics of the haptic interface to match those of the virtual tool or of the proxy when the synchronous or the asynchronous communication architectures are used, respectively. Both controllers are presented in further detail in Appendix D. A l l experiments later described in this thesis use the planar haptic simulation system described in this chapter and the synchronous and the asynchronous testbed virtual environments depicted in Figures 3.3(a) and 3.3(b), respectively. In both virtual environments, the moving objects have mass m = 2kg and moment of 1  www.criticalmasslabs.com  30  Haptics processor Virtual (Controller environment simulation  Graphics processor Ethernet  Haptic device  Graphical display  (a) Synchronous communication.  Haptii :s processor  o  Controller  Local model of interaction  *9'  Ethernet  Haptic device  Graphical display  (b) Asynchronous communication.  Figure 3.2: Communication between the virtual environment simulation and the force control loop in the planar haptic interaction system used in this thesis.  inertia / = 0.005kg-m . The links of the virtual linkage in Figure 3.3(b) are numbered from the proximal link 2  to the distal one. They have masses m\ — 3kg, m-i = 3kg, 7713 = 1kg, moments of inertia I\ = 0.015kg-m , 2  I2 = 0.015kg-m , ^3 = 0.005kg-m , and lengths h = 4.2cm, I2 = 4.2cm, I3 = 3cm, respectively. 2  2  Several experiments presented in later chapters require performance comparisons between the simulation and control methods developed in this thesis and prior approaches. Since users are not able to apply the exact same interaction wrenches during successive trials, these experiments cannot involve manipulations performed by humans. To enable comparison, these experiments are controlled experiments. The same initial conditions and the same "user" are ensured during one experiment through replacing the user's hand by controlled wrenches. Since the haptic device is an impedance-type interface, elimination of the adaptive damping associated with the user manipulation of the device represents a worst-case scenario for stability [64].  31  (a) V i r t u a l environment synchronous with the haptic loop.  (b) V i r t u a l environment asynchronous with the haptic loop.  Figure 3.3: Testbed virtual environments used in the experiments.  Chapter 4  H a p t i c rendering of rigid contact A n impulse-augmented penalty simulation approach is proposed in this chapter that enables users to feel collisions during their haptic manipulation of rigid multibody virtual environments. The chapter starts with a synopsis of the proposed approach and of the haptic controller that applies the simulated interactions to users. Then, the rigid body contact model underlying the impulse-augmented penalty approach is discussed. The dynamics employed during contact are subsequently presented. These dynamics include a friction model developed in computational dynamics [61] that is for the first time employed for haptic rendering of dry friction, and is compared through simulations and experiments to existing models. The dynamics used upon contact follow. Next, the performance of the impulse-augmented penalty simulation approach is compared to the performance of existing approaches through simulations. Limitations imposed by the haptic device are discussed and a solution is devised to address them. Lastly, experiments are presented to validate the simulation results presented in earlier section. Coordinate invariant representation of contact at the user's hand is deferred to the following chapter. The performance of the impulse-augmented penalty simulation approach is contrasted to the performance of prior approaches through investigating the realism and the stability of the haptic interaction. The realism of the interaction is assessed through evaluating the transient response of the closed loop system that comprises the user, the device, and the virtual environment simulation. A closed loop system with smaller overshoot and shorter settling time represents more realistic interaction. This is because smaller overshoot means that users violate the virtual constraints less, and shorter settling time means that users stop faster upon contact with a static virtual environment. In other words, the trajectory imposed on users by the virtual constraints is closer to 32  4.1 Synopsis of the impulse-augmented penalty simulation and of the haptic controller  33  the trajectory that a physical rigid environment would impose on them. The stability of the haptic interaction is defined as the coupled stability (in the sense of limited velocities and forces [39]) of the closed loop system that comprises the user, the device, and the virtual environment simulation. Since a feedback interconnection of passive systems is necessarily stable [148] and the haptic device and the user are passive [69], the passivity of the simulation guarantees the coupled stability of the haptic interaction. Furthermore, enhanced passivity of the virtual environment enhances the passivity of the closed loop system, thereby improving the stability of the haptic interaction. Therefore, interaction stability is assessed in this thesis through examining the passivity of the virtual contacts by monitoring the kinetic energy of the user's hand during the interaction. The schematic of haptic interaction within rigid multibody virtual environments generated using the methods introduced in this chapter is shown in Figure 4.1. In this figure, and torque), F  env  is the user-applied wrench (force  is the environment wrench represented at the user's hand, and fi and  are contact forces  between the virtual tool and other virtual objects. Furthermore, x/j and x/j are the body position (position and orientation) and the body velocity (linear and angular velocity) of the user's hand, and x j / i and x im,h g  m )  S  are simulated body position and body velocity of the user's hand. Haptic device  Haptic controller  Coordinate invariant contact representation  F  Impulse-augmented penalty simulation I*  4 —~.—, ^  X  ,—"  j y .  w  —"-*\.f  sin,.h  -/AVFigure 4.1: Haptic manipulation in an impulse-augmented penalty virtual world.  4.1  Synopsis of the impulse-augmented penalty simulation and of the haptic controller  The impulse-augmented penalty simulation approach is motivated by user studies [90,125] that demonstrate that the perceived rigidity of virtual contacts can be improved by applying large forces to users upon contact. The approach computes impulsive forces upon contact and then penalty and friction forces during contact. The impulsive forces are derived using a new collision resolution method that never increases the kinetic energy of the simulated system. When new contacts arise, the impulsive forces generate large hand accelerations without  4.1 Synopsis of the impulse-augmented penalty simulation and of the haptic controller  34  requiring increased contact stiffness and damping. The approach generalizes earlier work in [131] in two ways: (i) it allows the energy dissipated upon contact to be adjusted through the coefficient of restitution; and (ii) it is suitable for rigid body manipulation, as opposed to point interaction. Compared to earlier work in [43] and [44], it uses a new contact model and Newton's restitution hypothesis for collision resolution. The contact model allows transitions to collision from all contact states. Therefore, rigid body contact is represented more accurately (see Section 4.2) and the simulation can account for device limitations (see Section 4.5.2). Newton's restitution law allows multiple collisions to be resolved such that the kinetic energy of the system does not increase during collisions. Manipulation of both virtual objects and linkages is enabled through simulating the dynamics of the virtual environment in configuration space. In this space, virtual objects and linkages are represented as points and no distinction between them is necessary. Ultimately, transparent haptic interaction is achieved by applying the simulated contact forces to users through a suitable haptic controller. The controller used for interaction with impulse-augmented penalty virtual worlds is a four channel teleoperation controller [136]. It coordinates both positions and forces between the virtual tool and the haptic device, as shown in Figure 4.1. The two position coordination channels implement a generalized (translational and rotational) spring-damper connection (i.e., a proportional derivative - P D controller) between the virtual tool and the device, maintaining kinematic correspondence and eliminating drift. The two force coordination channels apply the hand wrench (force and torque)  to the virtual tool in  the simulation, and the environment wrench to the user. Unlike haptic controllers that employ only position coordination and apply just penalty-like wrenches to the user's hand [3,42], the four channel teleoperation controller directly applies the simulated impulsive interactions to users through its force channels. Hence, it enables users to feel collisions and improves the stability and the perceived rigidity of the virtual contacts. Figure 4.2 illustrates the behavior of the impulse-augmented penalty simulation: upon contact (Figures 4.2(b) and 4.2(d)), the rigid body contact model exactly enforces contact rigidity, while during contact (Figures 4.2(c) and 4.2(e)), contact rigidity is only approximated.  35  4.2 The contact model  (a)  (b)  (c)  (d)  (e)  Figure 4.2: The impulse-augmented penalty-based rigid body contact model: contact rigidity is exactly enforced upon contact and only approximately enforced during contact.  4.2  The contact model  The impulse-augmented, penalty-based, rigid-body-contact model is a dynamic model rather than a geometric one. It is used for computing interaction forces and impulses, based on the geometric information provided by a collision detection algorithm that is implemented in the virtual environment. Typically, the collision detection algorithm decomposes each rigid object into a collection of convex polyhedra and computes contacts between pairs of these polyhedra [85]. In three dimensional (3D) virtual environments, these contacts represent either vertex-face or edge-edge contacts , Figure 4.3. For each contact, it provides a contact point c , a penetration 1  depth (equal and opposite to the separation distance s between the bodies), and a contact normal direction n.  (a) Vertex-face contact.  (b) Edge-edge contact.  Figure 4.3: Contact information received from collision detection.  1  A l l other types of contacts, for example face-face contacts, are represented as a finite number of vertex-face a n d / o r edge-edge  contacts [152].  4.2 The contact model  36  In the impulse-augmented penalty simulation approach, a contact is defined by this geometric information plus the contact velocity v. The contact velocity is the relative velocity between the contacting polyhedra at the contact point. It is defined such that the normal contact velocity v  n  (i.e., the component of the contact  velocity along the contact normal v = n • v) is negative if the polyhedra move into each other. A penetrating T  n  contact is a contact with negative normal contact velocity v  n  < 0. A separating contact is a contact with  nonnegative normal contact velocity v„ > 0. Given a rigid multibody virtual environment, two rigid bodies are said to be in the same contact group if there exists a chain of contacting moving rigid bodies between them [122]. The rigid body contact model has 2  three states: free motion, colliding contact, and resting contact. A rigid body is said to be in free motion if it has no contacts. A rigid body is said to be in colliding contact if at least one new penetrating contact exists within its contact group. Finally, a rigid body is said to be in resting contact if it is neither in free motion nor in colliding contact. Note that a body can have non-zero acceleration and velocity during "resting contact". This terminology is used in the present work to maintain consistency with prior literature on the subject [152]. Note that, due to the fixed time step of the haptic simulation, bodies may transition from free motion to resting contact. A n example is depicted in Figure 4.4, where contact A does not exist at time t and is already a separating contact at time t + 1. Hence, collisions might be missed in an impulse-augmented penalty virtual environment. Note also that bodies may remain in colliding contact for multiple consecutive simulation steps because new penetrating contacts may appear at every step.  Figure 4.4: Example transition from free motion to resting contact during a fixed step haptic simulation, w is the angular velocity of the body, v is the velocity of its center of mass, v^ is the velocity of point A on the body, and t is the time step of the simulation.  By including the colliding contact state, the impulse-augmented penalty-based rigid body contact model approximates rigidity better than the penalty-based model. The tradeoff is that, by allowing body interpene2  a s opposed to static bodies, i.e., rigid walls.  4.3 Resting contact  37  tration during resting contact, the model provides a poorer approximation of rigidity than the constraint-based contact model. The combination of approximate and exact constraint enforcement requires penalty-based resting contact dynamics and constraint-based collision resolution. The numerical methods employed to simulate resting contact are presented in the following section. The simultaneous resolution of multiple collisions is discussed in Section 4.4.  4.3  Resting contact  This work considers the case in which users manipulate both virtual objects and virtual linkages (i.e., the virtual tool can be a single object or a chain). Realistic forces during these types of interaction can be computed by representing the virtual world dynamics either in Cartesian space or in configuration space. In Cartesian space, constraint equations must be added to maintain the bilateral constraints and the simulation must integrate a computationally expensive differential algebraic system of equations for which constraint satisfaction may be problematic. In configuration space, the bilateral constraints are embedded in the coordinate representation. In this case, only a reduced number of coordinates must be integrated and bilateral constraint satisfaction is guaranteed. Therefore, configuration space dynamics are used in the proposed simulation to compute the interactions between the virtual tool and the virtual world. Since bilateral constraints are incorporated in the coordinate representation, only contact and user applied forces must be included in the dynamics equations. Consider a contact group with d degrees of freedom (DOFs) and c resting contacts. In configuration space, its dynamics are: c  (4.1)  In Equation (4.1),  D (q) € 7Z  dxd  is the configuration space inertia matrix of the contact group,  represent Coriolis and centripetal effects, G (q) G lZ are the gravitational terms, J; (q) € TZ d  computed at the i-th contact,  3xd  B (q, q) S lZ  d  is the Jacobian  (q) € 7Z  fj € TZ is the Cartesian space contact force at the i-th contact,  6xd  3  is the Jacobian computed at the user's hand, f/j € TZ being the user-applied force and r 3  = (fj  T  £ TZ is the wrench applied by the user (with 6  £ TZ the user-applied torque), and q € lZ , q G lZ , and q € lZ 3  h  T£)  d  d  are the configuration space positions, velocities, and accelerations, respectively (see Figure 4.5). If the contact group consists of both virtual objects and virtual linkages, the matrices and vectors in  d  4.3 Resting contact  38  J ,(q)f,+J ('M+JTH(q)Fj3 T  T  2  U (q)f+J (q)f+J (q)FKl2 T  1  1  T  2  2  T  k  Pi  U (q)f +J (q)f +JTH(q)Fj T  1  1  T  2  2  Figure 4.5: Example contact forces and and hand wrenches arising during the haptic manipulation of a contact group with dynamics computed in configuration space.  Equation (4.1) are obtained by concatenating the matrices and vectors corresponding to each object and linkage. For example:  and: D i (qi) D(q) =  (4.3) 0  •••  D (q )J m  ra  In Equation (4.3), m is the total number of virtual objects and virtual linkages and the configuration space dynamics of a virtual rigid body are the same as its Cartesian space dynamics. The contact forces in Equation (4.1) have a component contact rigidity, and a component  along the contact normal direction rij, modeling  along the direction t* (orthogonal to the contact normal), modeling dry  friction: U = /n,»nj + / / , i t j .  (4.4)  Resting contact is enforced using penalties. Hence, the normal component of the contact force at the z-th contact is computed by: fn,i — —KcontactSi (q) — B tactV ,i con  In Equation (4.5), K  c o n t a c t  and B  c o n t a c t  n  (q) •  (4-5)  are the contact stiffness and damping, Sj is the separation between  bodies at contact (because bodies overlap, Sj is negative and equal to the penetration depth of the contact), and i> ,i is the normal contact velocity. Dry friction can be simulated using any friction model that employs n  4.3 Resting contact  39  only local contact state information, such as [67] or [113]. However, the reset-integrator model proposed by Haessig and Friedland [61] is implemented in this work, for reasons discussed in the following section.  4.3.1  Friction modeling  The reset-integrator model was first used in haptics in [43]. According to this model, the friction force at each contact is computed by:  r  (4.6) (1 + a)K z + Bvt  if stick  Kz  if slip  r  T  and the input to the integrator is: 0  Vt >  i= 0  if  {  z >  ZQ  (slip)  or vt < 0  z = vt  and  and  (4.7) z < —zo (stick)  otherwise.  In Equation (4.6), z is the strain of the bond between the contacting bodies, zo is the breakaway distance (i.e., the maximum strain of the bond), K is a spring stiffness, a is the "stiction gradient", v is the sliding velocity T  t  (i.e., the projection of the contact velocity on the contact plane vt = t v), T  and B is a damping coefficient.  In a haptics implementation, the model may exhibit chattering because of noisy velocity measurements. This difficulty is avoided by implementing a dead band of width 2v i m  n  around zero sliding velocity. Specifically,  Equation (4.7) is replaced by:  V  >  t  and  -Vr,  z >  ZQ  ii {  i= 0  (slip) (4< v < v  -Vmin  z =  V  t  t  otherwise.  min  and  z <  -z  0  (stick)  In the reset-integrator model, the contacting points on the two bodies are connected by a spring with stiffness (1 + a)K  r  during sticking and K  r  during sliding. A stiction gradient a = 0 describes contacts which  have a coefficient of static friction fi equal to the coefficient of kinetic friction fif s  A contact transitions  4.3 Resting contact  40  from stick to slip when the strain of the bond exceeds the maximum strain zo and tends to increase further. Otherwise, the contact transitions to slip. The oscillations occurring when entering stiction are damped by the j3vt term. The role of the damping factor in controlling the transition from slippage to stiction is illustrated through simulating the peg-in-hole manipulation depicted in Figure 4.6. In the simulation, the peg has mass m = 2kg and barely fits the hole, the normal contact force at each contact is / „ — 0.175N, and all contacts have the same coefficients of static /J. = 0.25 and kinetic /xj, = 0.2 friction. The breakaway distance is ZQ = 0.01mm and S  the velocity dead band is v i m  n  = O.lmm/s. To match the experiments, the simulation time step is chosen gj^s.  The simulated position of the C O M of the peg, the sliding velocity at each contact (equal to the velocity of the C O M of the peg), and the resultant friction force on the peg are shown in Figure 4.7 for three values of the damping coefficient: (i) no damping, j3 = 0; (ii) critical damping, j3 = 118; and, (iii) overdamping, j3 = 1180. For clarity, the figures depict the sliding velocity and the friction force only during the first 1.5 seconds of the first stiction period. |/= -0A-0Asm(t) [N]  Centre of Mass (COM)  Figure 4.6: Simulated peg-in-hole manipulation.  The results in Figure 4.7 illustrate that the reset-integrator model successfully renders the slip-stick phenomenon regardless of damping. In addition, underdamping results in oscillatory sliding velocity and friction force during stiction. Critical damping and overdamping eliminate these oscillations when transitioning from slippage to stiction. In haptics, the friction force oscillations may be used for changing the feel of the contact. Psychophysical studies are required to clarify the relationship between the model damping and the user's perception of the contact characteristics. By representing stiction through compliance in the contact plane, the reset-integrator model is similar to other haptic dry friction simulation techniques reviewed in Section 2.3.4, such as the model based on the human finger pad characteristics proposed by Nahvi et al. [113] and the modified Dahl method proposed by Hayward  4.3 Resting contact  41  COM position along y-axis  4.5 £ o o  3.5 3  2.5 2.  no damping critical damping overdamping 10 time [sec] (a) C O M  COM velocity along y-axis  15  20  position.  Friction force  0.01 |0.005 o  0  -0.005 -0.01 -0.015 -0.02.  no damping critical clamping overdamping 0.5 1 time [sec] (b) S l i d i n g velocity.  0.5  time [sec]  1  (c) R e s u l t a n t f r i c t i o n f o r c e .  Figure 4.7: Peg-in-hole simulation. Friction is simulated using the resetintegrator model [61] with various damping coefficients.  and Armstrong [67]. Like these methods and unlike the classical approximation of Coulomb friction (Figure 4.8), the model does not exhibit drift. In addition, it can either render oscillations at the transitions from slippage to stiction similar to the models proposed by Nahvi et al. [113] and by Hayward and Armstrong [67], or it can eliminate these oscillations similar to the classical Coulomb model. This is illustrated through simulating the manipulation shown in Figure 4.6 using the critically damped reset-integrator technique, the models proposed  4.3 Resting contact  42  by Nahvi et al. and Hayward and Armstrong, and the classical approximation for implementing friction. In the model based on the human finger pad characteristics [113], the slipping friction force f/ is computed by: f)=Md||n||^,  (4.9)  where p<j is the coefficient of dynamic frcition, ||n|| is the magnitude of the normal force, and v is the tangential t  velocity of the haptic interface. Stiction begind when ||v || becomes smaller than a threshold velocity t  During stiction, the haptic interface is trapped in a circle of radius  v i„. m  centered at c by a spring of stiffness kf.  When the interace transitions from slippage to stiction, the stiction centre c is computed by:  c =a - —-—-,  (4.10)  kf I N I where a is the position of the haptic interface. The spring has a rupture limit of p ||n|| and slipping begins s  when the haptic interface is  away from the stiction centre.  In the modified Dahl model [67], the dry friction force f/ is computed by:  f y ^ l H l ^ - i - , ||v || z t  (4.11)  max  where z is an internal position variable representing the distance from the sticking centre to the haptic interface and  z ax m  is the maximum value of this distance. Stiction occurs when z < z  , z is updated using forward  max  Euler integration, and the time rate of change of z is computed by:  vt  if stiction, i.e., ||z|| < z  max  z=< v (lt  The  J^t: 9 ( t)} si  n  v  (4.12)  if slippage, i.e., ||z|| > z„  simulation results are depicted in Figure 4.9. For clarity, the sliding velocity (equal to the peg  COM  velocity) and the friction force are plotted only during the first 1.5 seconds of the first stiction period. The  damping included in the reset-integrator model allows an extra degree of freedom in designing the  contact characteristics: the higher the damping, the faster the oscillations occuring at the transition between slip and stick are suppressed. Therefore, this model is considered most suited for force feedback and proposed for use in haptics applications.  4.3 Resting contact  43  M/. v, >  Figure 4.8: Classical approximation of Coulomb friction.  To experimentally validate the reset-integrator dry friction model, the same controlled interaction is implemented on the haptic simulation system presented in Chapter 3, using the synchronous communication paradigm. In the experiment, users insert the peg into the hole, release the device handle, and their hand is replaced by the unidirectional sinusoidal force f = (—1 — sin(t))N aligned with the hole axis. A l l investigated y  friction models are used for simulating contacts with friction. Their parameters are adjusted such that the contacts have coefficients of static and kinetic friction [i = 0.25 and Hk = 0.2, respectively, a breakaway distance s  zo = 0.01mm, and a velocity dead band v i m  n  = O.lmm/s. The reset-integrator model is critically damped,  (/? = 118). The trajectory of the peg C O M and the forces applied to users along the hole axis are plotted in Figure 4.10. The experimental results confirm that all friction models exhibit slip and stick. However, some transitions from stick to slip are missed when dry friction is simulated using the model based on the human finger pad characteristics. The experiments also demonstrate that the driftless Dahl model and the model based on the human finger pad characteristics introduce perceivable compliance during stiction. This might not always be desirable for rendering frictional contact between rigid bodies. The critically damped reset-integrator model successfully eliminates such compliance. At the same time, the model can easily be adjusted to allow users to perceive oscillations during stiction by decreasing its damping coefficient 8. Hence, the reset-integrator friction simulation technique can render a wider range of behaviors during stiction. The advantage of the reset-integrator model over the classical friction model is that it does not exhibit drift [61]. Presently, the reset-integrator model of dry friction is implemented by critically damping the oscillations during stiction. A physically meaningful choice of a damping coefficient requires future user studies.  The friction and contact models in Equations (4.5) and (4.6) show that only state information is used  4.3 Resting contact  44  COM position along y-axis  4.5 ^ E  4  o .5 o 3  reset-integrator Nahvi et al. modified Dahl classical  2.5 2  10 time [sec]  15  20  (a) C O M position.  Friction force  COM velocity along y-axis  §0.005  0.5  time [sec]  (b) Sliding velocity.  1  0.5  time [sec]  1  (c) Resultant friction force.  Figure 4.9: Peg-in-hole simulation. Friction is simulated using the critically damped reset-integrator model [61], the model proposed by Nahvi et al. [113], the modified Dahl model [67], and the classical dry friction approximation.  in the impulse-augmented penalty simulation to compute the contact forces in Equation (4.5). Therefore, Equation (4.1) can be directly solved for the configuration space acceleration:  (4.13)  4.3 Resting contact  45  COM position along the y-axis User^s hand is replaced by controlled force  COM position along the y-axis  Oscillations during sticking  Missed stick->slip transition 10  15  20  25  30  Force along the y-axis  10  (c) M o d e l b a s e d o n h u m a n  finger-pad  c h a r a c t e r i s t i c s [113].  15  time [s]  20  (d) R e s e t - i n t e g r a t o r m o d e l [61].  Figure 4.10: Forces felt by users and C O M trajectories of a virtual peg sliding with friction in a tight hole due to a sinusoidal force applied to its C O M . Various friction models are used to simulate dry friction. Note that all models render slip and stick. Perceivable oscillations occur during stiction when friction is simulated using the modified Dahl model [67] (Figure 4.10(b)) and the model based on human finger-pad characteristics [113] (Figure 4.10(c)). Some slip to stick transitions are missed by the model based on the human finger-pad characteristics [113] (Figure 4.10(c)).  4.4 Colliding contact  46  In Equation (4.13), the dependence on the instantaneous state of all terms on the right hand side of the equation is implied. The configuration space acceleration is then integrated using a fixed step size integrator compatible with the requirements of the haptic control loop. The dynamics of resting contact discussed in this section are typical for penalty-based simulations. They are used to approximately enforce rigidity during contact. Exact enforcement of rigidity upon contact is accomplished through implementing constraint-based colliding contact dynamics. These dynamics are presented in the following section.  4.4  Colliding contact  The colliding contact state is introduced in the impulse-augmented penalty approach to improve the penaltybased approximation of unyielding contacts. In conjunction with the fixed step size of the haptic simulation, this state requires multiple collisions to be resolved simultaneously. Rather than incorporating contact rigidity and dry friction into a complementarity formulation, the proposed approach uses three simplifying assumptions: • that impulses develop at all contacts where bodies move into each other during a collision; • that velocities of all colliding contact points obey Newton's restitution law; • and, that collisions are frictionless. Unlike more accurate multiple collision models based on complementarity formulations [6,63,112,120], this new technique is non-iterative and requires no further assumptions on the value of the coefficient of restitution or on the shape of colliding objects. The colliding contact dynamics are obtained through time integration of Equation (4.1):  Dq = D q + £ ] f j f f ^ D q o 0  »=1  J  t  °  +  XXpi.  In Equation (4.14), D q and D q are the pre- and post-collision configuration space momenta and 0  ( ) 414  i=l  = / ' fjdi t  is the impulse at the i-th contact. Since collisions are modeled as instantaneous events, i.e., t -> to, the hand wrench and the gravitational forces do not contribute impulses to the impulse and momentum balance of the system. Furthermore, apart from the collision impulses, no other external impulses are applied to the contact group.  4.4 Colliding contact  47  In addition, collisions are assumed frictionless:  (4.15)  Pi^PiUi.  In Equation (4.15), pi is the magnitude of the i-th contact impulse. Then, the configuration space dynamics of colliding contacts become: c  Dq^Dqo + ^ j f n ^ i =Dq  0  + jjp,  (4.16)  »=l  where:  (4-17)  P=(vi...v<...v) ...p ...p, t  is the vector of contact impulses, and:  (4.18)  Jc =  Jfri! . . . Jf n, .. . J j n  c  is the contact Jacobian. For a contact group with d DOFs and c colliding contacts, Equation (4.16) represents a set of d equations with d + c unknowns, the post-collision configuration space velocity q and the contact impulses p. The additional assumptions needed to solve this system are provided by the various collision laws proposed in the literature [112], [63], or [6,120]. In this work, Newton's restitution rule is used, because it allows the development of a non-iterative solution that imposes no restrictions on the coefficient of restitution or the shape of the colliding objects. Moreover, the kinetic energy of the contact group does not increase during collisions. For one colliding contact, Newton's restitution hypothesis relates the pre-collision (v o) and post-collision n  (v ) normal contact velocities through the coefficient of restitution e: n  v = -ev . n  (4-19)  n0  The coefficient of restitution e € [0,1] describes the nature of the collision, with e = 1 corresponding to a perfectly elastic collision (no energy loss), and e = 0 corresponding to a perfectly plastic collision. In configuration space, Equation (4.19) becomes:  n J(,,q = -en J .qo T  T  6l  (4.20)  4.4 Colliding contact  48  at a collision between body 6, of the contact group and a static environment, and:  n  T  (J  bi  - J .) q = - e n  T  b  (J  - J) q  bi  bj  0  (4.21)  at a self-collision, i.e., a collision between bodies bi and bj of the contact group. A more restrictive condition is imposed at a self-collision in the proposed approach. Namely, the second simplifying assumption is imposed on the contact group by requiring it to obey:  n J q=-en J  q  0  (4.22)  n J q = -en J q  0  (4.23)  T  r  6 i  f t i  and: T  T  6 i  6 i  simultaneously. Equations (4.22) and (4.23) ensure both that Newton's restitution law is observed and that the proposed collision resolution technique maintains system passivity, as shown in subsequent derivations. Using Equations (4.20), (4.22), and (4.23), Newton's restitution law is restated as:  JA = -eJAo-  (4.24)  Equation (4.24) represents a set of c equations, where self-collisions are counted once on each colliding body. Note also that the equal sign in Equation (4.24) embeds the first simplifying assumption used for resolving collisions. In the impulse-augmented penalty simulation approach, Equations (4.16) and (4.24) describe the dynamics of colliding contact of a contact group with d DOFs and c simultaneous collisions. Their resolution, theorems for system passivity for both independent and overdetermined constraints, and the computation of the configuration space impulsive interactions applied to users are presented in the following subsection.  4.4.1  Passive collision resolution  This section starts by showing that a contact group with one frictionless contact is passive if its colliding contact dynamics are resolved using Newton's hypothesis. The result is then extended to a contact group with multiple independent and overdetermined contacts for which Equation (4.24) is used to ensure that Newton's  4.4 Colliding contact  49  collision law is obeyed. 4.4.1.1  Independent constraints  Passivity of a contact group with a single colliding contact is shown by proving that: Theorem 1 If a contact group described by the momentum equation: Dq = D q + J p T  0  c  has one frictionless colliding contact and the post-collision normal contact velocity is given by: v = -ev , n  n0  where e € [0,1] is the coefficient of restitution, then the post-collision kinetic energy of the contact group is: KE = KE - (1 - e ) q j jJ J T> J J q 2  T  C  0  where J = ~D~ J^ ( ^ D l  c  Proof  - 1  JJ)  c  c  0  < KE  (4.27)  0  is the dynamically consistent inverse of J [81].  1  c  The proof starts by computing the contact impulse by substitution of Equation (4.24) into Equa-  tion (4.16):  -eJ qo = Jcqo e  + JcD~ JcP l  =>  P = " (1 + e) (JcD"  1  jJY  1  JAo-  (4.28)  Then, the post-collision configuration space velocity results after substitution from Equation (4.28) in Equation (4.16):  q = qo - (1 + e) D " JJ (JcD^jJY 1  JAo = (I - (1 + e)J Jc) qo  1  e  (4.29)  where I is the d x d identity matrix. Next it is shown that: DJcJe - jJj VJ J T  c  c  c  = 0.  (4.30)  4.4 Colliding contact  50  Indeed, using the definition of the dynamically consistent inverse of J and the symmetry of the inertia matrix: c  ~D~J Jc - JJJ VJCJC  =  T  c  C  =  D D  =  jj  (JcD^jJ)  =  jj  {jcD^jjy  1  jj  je - JI (p~ Jl  {j -D- jjy l  v  _1  l  (JeD" jj) ~ J T>~ jJ  j  {jcD- jjy  c  - jj  In Equation (4.31), (JcD" jj)~ 1  {DJCJC  l  J - jj c  1  {jjy- jjy ) vT>-  l  c  T  1  T  c  l  v  j  l  c  c  =  J =  1  c  1  r  j  1  c  (4.31)  = o.  1  = JJJ ~D - jJj VJ J  T  T  {Jji^jjy  [j T>- jjy  and D ~ = D " since they are symmetric. Furthermore:  = {JcD^jJ)'  - JJJ VJCJC)  jj  x  (jcD- j j y  jcD-ij?  l  T  T  = 0.  T  C  C  c  c  c  (4.32)  Then, the post-collision kinetic energy of the system is be computed as follows:  = iq Dq=i  KE  Since jj J DJ Jc c  c  " ( l - ( l +e)J jJ)D(l-(l+e)J Jc)qo = T  T  q o  =  i q o Dqo - (1 + e) hfi (nj J  =  KE - (1 + e) h$ (2jJj VJ J  =  KE - (1 - e ) i&jj J VJ JA .  c  + jJj T> T  c  c  T  0  c  2  0  c  c  c  c  - (1 + e) jJj VJ J ) c  c  c  - (1 + e) jJj TVJ J ) c  c  C  0  q =  c  0  (4.33)  T  C  q =  T  0  is symmetric, it is positive semi-definite and:  KE < KE  0  Ve £ [0,1]  Q.E.D.  (4.34)  The contact impulse due to one collision given in Equation (4.28) is equal to the contact impulse computed in prior complementarity formulations [112]. Hence, the simplifying assumptions embedded in (4.24) involve no further approximation for the case of a single frictionless colliding contact. Moreover, the proposed collision resolution method uses the dynamically consistent inverse of the collision Jacobian, i.e., it is coordinate invariant. Kinetic energy is conserved during a perfectly elastic collision (e = 1). The loss of kinetic energy during a plastic collision (e < 1) depends both on contact properties, as given by the coefficient of restitution e, and on contact geometry and the contact group topology and geometry, embedded in J . Various contact group c  4,4 Colliding contact  51  topologies and geometries and various contact geometries result in either a total or a partial loss of kinetic energy during a plastic collision. This is illustrated in Figure 4.12 for two different contact geometries of the two link planar manipulator depicted in Figure 4 . I I . Note that the manipulator loses all kinetic energy during 3  a perfectly plastic collision between the distal link and a static environment, and it loses only a part of its kinetic energy during a perfectly plastic collision between the base link and a static environment.  (a) Base link colliding with the environment.  (b)  Distal  link colliding with the  environ-  ment.  Figure 4.11: Two link planar manipulator whose loss of kinetic energy for various values of the coefficient of restitution is shown in Figure 4.12 for two contact geometries. The simplifying assumption in Equation (4.24) imposes the particular solution shown in Figure 4.12(b).  Pre- and post-collision kinetic energy  (a)  Base link colliding with  environment (L\  Pre- and post-collision kinetic energy  the  (b) Distal link colliding with the  — 0.7m).  environment (L2 = 1.5m).  Figure 4.12: Loss of kinetic energy of the planar two link manipulator in Figure 4.11 during one frictionless collision for two contact geometries.  3  T h e manipulator has link lengths h — h = l m , link masses m i =  and pre-collision configuration space velocity qn = (1 and l m < L2 < 2m. respectively.  1)  = 1kg, configuration space position q = (5  0)  r a d / s . T h e environment constraints at A and B are such that L\  rad, <  lm  4.4 Colliding contact  52  The collision resolution method employed for one colliding contact can be directly applied to resolve multiple collisions simultaneously if the contact constraints are independent, i.e., J is full row rank, and Equation (4.24) c  is imposed to ensure that the post-collision configuration space velocity obeys Newton's hypothesis. Then, J is c  full row rank and the matrix JcD~ jj  is invertible. This can be easily shown by considering that its singular  l  and J jJ, Smin ( D ) and S  values of bmD~  x  _1  c  q^cD- J q > S 1  where E  (D- ) q J J  T  c  1  m i n  T  c  T c  q > S  qq > 0 Vq # 0,  (4.35)  . Equation (4.35) proves that J ~D~ jJ  is positive  (D" ) S 1  r a i n  ( D ) is the minimum singular value of D _1  m i n  respectively, obey:  {JcJj),  m i n  _ 1  m i n  {JcjJ)  T  l  c  definite, hence invertible. 4.4.1.2  Overdetermined constraints  If the contact constraints are overdetermined, J is rank deficient and the matrix J T>~ J^F is not invertible. 1  c  Nevertheless, its pseudo-inverse ( 7 D l  - 1  c  C  7 ) * can be used to compute the contact impulses p, and the postT  t  c  collision configuration space velocity q according to:  p  =  - ( 1 + e) ( j D - J  q  =  qo-a  1  c  T c  ) J q  (4.36)  t  c  0  + ^ D - ^ ^ D - ^ V c q o -  (4.37)  Passivity of a contact group with overdetermined colliding constraints resolved simultaneously according to Equation (4.37) results from the following theorem: Theorem 2 If a contact group described by the momentum equation: Dq = Dq + J p,  (4.38)  T  0  c  has c overdetermined frictionless colliding contacts and its post-collision configuration space velocity is given by:  JA = -eJAo,  (4.39)  where e € [0,1] is the coefficient of restitution, then the post-collision kinetic energy of the contact group is: KE = KE In Equation  (4-40), J  n  =  D~ J J 1  r  - (1 - e ) q ^ J 2  0  ( J n D jl) - 1  1  T n  J ^ D J „ J „ q o < KE .  (4.40)  0  is the dynamically consistent inverse of  J , and jj n  =  4.4 Colliding contact  53  T  3n  , with J full row rank.  3  n  r  P r o o f The proof follows the same reasoning as the proof of Theorem 1, where (j T> jj) 1  1  c  (J T>- J^.  is replaced by  The proof holds since:  1  C  (JcD" J )t 1  - i jTf  T  c  =  (  J  c  D  J  c  D  - i jT ( ^ D - i j ^ t  (  4 4  1  )  Furthermore, it is shown in Appendix A that:  Jj ( J c D - J 1  where jj  =  KE  JJY  n  d 3n  full  1S  r  o  w  ) J = Jl ( J n D - ^ ) t  1  =  KEo - (1 - e )  (D- J  =  KE -(l-e )c£jJ  (J T>- jj)  =  K £ - (1 - e )  =  K £ - (1 - e ) o j j f jlBJ J q .  D  2  2  J  2  T C  I  ( ^ D "  q o  "j  1  1  J q  ]  c  0  T C  ( J C D "  (j' D- J ) J q  l  1  C  c  ( J c D " Jj)  T c  ^ ) ^ D (D-" J  1  J T>~ jJ  ]  C  2  c  JAO =  C  1  0  (4.42)  n  n  0  J,  rank, i.e., rank{J ) = rank{J ) = n. Then:  KE - (1 - e ) tg jj  0  - 1  c  =  c  r  t  c  0  1  J?)*) JAo =  =  = (4.43)  2  0  Since J% J DJ Jn n  a  T c  n  n  n  0  is positive semi-definite and e g [0,1], it follows that the kinetic energy of the contact group  does not increase during simultaneous collision resolution when the constraints are overdetermined regardless of the value of the coefficient of restitution, i.e. KE < KEo Ve € [0,1] Q.E.D.  Similar to the case of independent constraints, the loss of kinetic energy depends both on contact properties, through the coefficient of restitution e, and on the topology and geometry of the contact group and the geometry of contact, through J . Kinetic energy is conserved during perfectly elastic collisions (e = 1). n  From Equation (4.42), it also follows that:  q  =  qo-(l +e)D- J  c  =  q - ( l +e)D- J  n  1  T  1  0  (j D- J  c  (j D- J  n  1  c  T  1  n  T  ) Jcqo =  T  )~  +  1  J„q . 0  (4.44)  4.5 The performance of the impulse-augmented penalty simulation  54  In other words, the post-collision state of the contact group is the same regardless of whether the collisions are resolved using the pseudo-inverse technique or the constraint overdeterminancy is eliminated before collision resolution. Hence, the pseudo-inverse method is equivalent to selecting a set of independent constraints and simultaneously resolving the collisions at these contacts as described in Section 4.4.1.1.  4.4.2  Rendering collisions to users  Collision impulses are rendered to users as impulsive forces through the haptic controller shown in Figure 4.1 at the beginning of this chapter. The configuration space impulsive torques to be applied to users' hand T „ e n  are computed such that, when integrated over one time step of the haptic simulation, they induce the same change in the configuration space momentum of the contact group as the simulated collision impulses:  T,env  —  (4.45)  At '  In Equation (4.45), At is the time step of the haptic simulation and p are the contact impulses, computed according to Equations (4.28) and (4.36).  The key feature of the haptic interaction within the impulse-augmented penalty virtual world described above is that it enables users to feel collisions when new contacts arise. This is in contrast to existing haptic interaction paradigms which apply only penalty-like forces to users. To demonstrate the advantage of the proposed approach, the next section compares it to existing methods, as well as to the "ideal" case (i.e., physical interaction) through simulations and experiments.  4.5  The performance of the impulse-augmented penalty simulation  In this section, the performance of the impulse-augmented simulation approach is evaluated against the performance of the penalty-based and the constraint-based approaches through simulated and experimental user interactions within a planar virtual environment . The evaluation is initially carried out assuming that the 4  Note that the dynamics in Sections 4 . 3 a n d 4 . 4 are suitable for rigid b o d y interaction within spatial v i r t u a l environments. However, the approach is validated only for rigid b o d y interaction within a planar virtual world in this thesis. T h i s is due to the fact that only planar (fx, f , T ) a n d point (f , f , f ) haptic interaction systems are available in the Robotics a n d C o n t r o l Laboratory. Furthermore, full rigid body force and torque feedback is required (fx, fy, T during planar interaction and f , f , fz, T , T , T during spatial interaction) to guarantee that the kinetic energy of the user's hand will not increase during collisions as a result of the impulsive feedback. Unstable interaction m a y arise due to the lack of torque feedback. 4  y  Z  x  y  z  Z  X  Y  Z  x  y  4.5 The performance of the impulse-augmented penalty simulation  55  haptic interface can fully apply the simulated impulses to users. Device limitations are then taken into account and solutions are designed to address them. Figure 4.13 depicts the virtual environment employed in the simulations. Users manipulate a rectangular peg by applying forces and torques at its center of mass. The virtual world is connected to the haptic interface through a unilateral coupler [155], [127] when generated using a constraint-based approach, and through a four-channel teleoperation controller [136] when generated using the impulse-augmented penalty and the penalty-based approaches. The unilateral coupler extends to rigid body interaction the controller used in the god-object [155] and the virtual proxy [127] point interaction techniques . It implements two spring-damper 5  connections (one translational and one rotational) between the virtual tool and the device during constrained motion and is inactive during free motion. The teleoperation controller coordinates both forces and positions between the virtual tool and the haptic interface, as explained in Section 4.1. Schematics of the mechanical equivalents of both controllers for 1 D O F interaction are shown in Figure 4.14.  Figure 4.13: Planar virtual world used in simulations and experiments. Starting from rest and the position shown, the rectangular object is pushed into the lower right corner by a controlled constant force.  In the simulations, the user applies a wrench F ^ = (0.32N  — 0.4N 0 N m ) on the peg, i.e., they push T  the peg towards the lower right corner of the virtual enclosure. The user-applied wrench is such that the peg hits the vertical wall first. It then slides along this wall until it hits the horizontal wall and stops. This simple interaction is chosen because it illustrates the performance of the approach for overdetermined constraints. The peg has dimensions li = 0.042m and lo — 0.021m, mass m = 2kg, and moment of inertia I = O.OOSkgm . 2  The virtual walls have stiffness K n wa  = 4000N/m and damping B u = 30N/(m/s). Collisions are considered wa  perfectly plastic (e = 0). The stiffness and damping of the position coordination channels of the teleoperation controller connecting the proposed and the penalty-based virtual environments and of the unilateral coupler However, since it is not clear how the god-object and the virtual proxy simulations can be extended to rigid body interaction, the virtual world is evolved using forward dynamics algorithms similar to those developed in graphics [11]. 5  56  4.5 The performance of the impulse-augmented penalty simulation  Haptic device  Virtual environment  Unilateral coupler  Haptic device  Virtual environment  Four channel controller  VsAAA(b) The four channel controller.  (a) The unilateral coupler during constrained motion.  Figure 4.14: 1 D O F mechanical equivalents of the controllers connecting the haptic device and the virtual environment. Controller parameters are given in Table 4.1.  connecting the constraint-based virtual world to the device are given in Table 4.1. The parameters of the teleoperation controller are optimized for transparency [136], while those of the unilateral coupler are chosen to match the impedance of the virtual contacts. To match the planar interface used for experiments, the haptic device is modeled as an impedance device. Furthermore, the time step of the simulation and of the controllers connecting the haptic interface and the virtual environment is chosen equal to  gjjS.  In addition, the haptic  device is considered to have purely inertial dynamics and to be kinematically equivalent to the proxy, i.e., the local device controller is not modeled. The user's hand is modeled as a pure force source, which is a worst-case scenario for stability when an impedance device is used [64]. Table 4.1: Parameters of the teleoperation controller and the unilateral coupler connecting the virtual environment to the haptic device. Four channel Controller K B  c o o r d  COO  = (lOON/m  r d = (70N/(m/s)  lOON/m  Unilateral coupler 0.5N/rad)  70N/(m/s)  J  0.375N/(rad/s))  K , = (lOOON/m c p  B , = (50N/(m/s) cp  lOOON/m 50N/(m/s)  2.5N/rad)  i  0.125N/(rad/s))  The performance of the haptic interaction paradigm proposed in this thesis is compared to that of existing techniques in the following two sections.  4.5.1  Best performance  In this section, the planar interaction described above is simulated by assuming that the haptic interface can fully apply the virtual impulses to the user. Hence, the best performance potentially achievable during haptic manipulation of impulse-augmented penalty virtual worlds is investigated. Interface limitations are considered and addressed in the following section. Figures 4.15, 4.16, and 4.17 depict the users' hand trajectories, the forces felt by users, and the kinetic energy of the users' hand, respectively, while users move towards the virtual rigid corner, contact the vertical wall, and slide down along it until they contact the horizontal wall. In these figures, the various interaction paradigms are identified as follows. "Ideal", " I A P B " , and " P B " represent interaction with a constraint-based, an impulse-  4.5 The performance of the impulse-augmented penalty simulation  57  augmented penalty-based, and a penalty-based virtual environment, respectively, all rendered to users through the four channel controller shown in Figure 4.14(b). " C B " represents interaction with a constraint-based U C  virtual world rendered to users through the unilateral coupler shown in Figure 4.14(a). Note that only the P B interaction has been demonstrated in existing haptics research. Note also that the "ideal" trajectory represents the peg interaction with a perfectly rigid (i.e., real) environment. Therefore, this trajectory separates the free motion of the peg from its constrained motion. Positions of the C O M of the peg to the left and above the lowest point on this trajectory represent free motion of the peg (i.e., peg not in contact with the environment). Positions of the C O M of the peg to the right and below the lowest point on this trajectory correspond to peg penetration into the vertical and the horizontal wall, respectively (i.e., peg in contact with the vertical and the horizontal wall, respectively). Figure 4.15 demonstrates that users' penetration into the virtual walls is smaller when users interact with a world generated using the impulse-augmented penalty approach than when they interact with a penaltybased or with a constraint-based virtual environment. This means that the trajectory of the user's hand is closer to the ideal trajectory, i.e., the interaction is more realistic, when the virtual world is simulated using the proposed approach. The user-perceived forces are closer to the ideal forces, too. Users feel large wrenches when new contacts arise (see Figure 4.16), and they feel wrenches that only balance the small hand wrench F/j = (0.32N  — 0.4N 0 N m ) during contact. Figure 4.17 demonstrates that users lose more kinetic T  energy upon contact when they interact with an impulse-augmented virtual world than when they interact with a penalty-based or with a constraint-based virtual environment. This means that the impulse-augmented penalty contacts improve the passivity of the virtual world compared to the other simulation paradigms, thereby improving the stability of the haptic interaction. Note that, while the large impulsive forces improve the perceived rigidity of the virtual contacts, they may exceed the force capabilities of the device. Constraints imposed by the haptic device on the virtual manipulation of an impulse-augmented penalty simulation are considered in the following section.  4.5.2  Accounting for device limitations  Several techniques can be used to account for device limitations during haptic interaction with impulseaugmented penalty virtual environments: • Collision impulses can be saturated on the device. When this strategy is used, the simulation computes collision impulses according to Equation (4.36) and sends them to the four channel controller according to Equation (4.45). The controller then saturates the impulses to the maximum value that the actuators can apply to the device. Hence, full collision impulses are applied to the virtual tool in the simulation and saturated impulses are applied to the device. As a result, different amounts of kinetic energy are extracted from the virtual tool and the device during collision and the kinematic correspondence between the two can be significantly changed, depending on the system dynamics. Post-collision kinematic correspondence is re-established through the position coordination channel of the four channel controller. This channel is more compliant than the contact, and thus larger user violations of constraints can result.  4.5 The performance of the impulse-augmented penalty simulation  58  Hand trajectory  -163 -164 ^-165 E E  "^-166  o o  >»-167 -168 -16 1  x  C0M[  m m  ]  2.5  2  (a)  2  Hand rotation  x 10"  1-  -  0 =o-1  -  -  CO i _  CD-2  -3 -4 -5  /  -  -  J  -1 *  -  —  | / 1  ZJ  /  .  |  I  1  l 1  0.3  ideal IAPB PB CB uc  0.4  0.5 0.6 time [seel  0.7  0.8  (b)  Figure 4.15: Simulated hand trajectories obtained when constraint-based ("ideal"), impulse-augmented penalty-based ("IAPB"), and penalty-based ("PB") interactions are applied to users by the four channel controller, and when constraint-based interactions are applied by the unilateral coupler ( " C B " ) . The device applies the simulated impulses to users in one step. UC  59  4.5 The performance of the impulse-augmented penalty simulation  Figure 4.16: Simulated forces along the x-axis and torques applied to users by the four channel controller when the virtual world is generated using the constraint-based ("ideal"), the impulse-augmented penaltybased ("IAPB"), and the penalty-based ("PB") methods, and applied by the unilateral coupler when the virtual environment is generated using the constraint-based method ( " C B " ) . The device applies the simulated impulses to users in one step. UC  60  4.5 The performance of the impulse-augmented penalty simulation  Kinetic energy of the user's hand  time [s] Figure 4.17: Kinetic energy of the users' hand when constraint-based ("ideal"), impulse-augmented penalty-based ("IAPB"), and penalty-based ("PB") interactions are appllied to users by the four channel controller, and when constraint-based interactions are applied by the unilateral coupler ( " C B " ) . The device applies the simulated impulses to users in one step. UC  • Collision impulses can be scaled in the simulation. When this strategy is used, the simulation computes collision impulses according to Equation (4.36) and scales them to the maximum value that the device actuators can apply on users. The scaled impulses are applied to the virtual tool and sent to the four channel controller. As a result, the same amount of kinetic energy is extracted from the virtual tool and the device during collision and their kinematic correspondence is maintained. However, the dissipated kinetic energy is less than that prescribed by the coefficient of restitution and the contact geometry. The simulation is altered and larger user penetrations into constraints occur. • Collision impulses can be spread over several steps of the simulation. When this strategy is used, the simulation scales the collision impulses as explained above. However, the colliding contact group does not transition to resting contact if scaling is necessary. Rather, it transitions back to colliding contact and new collision impulses are computed at each step of the simulation until the force levels return to the range of the haptic device. As a result, the amount of kinetic energy extracted from the virtual tool and the device at each simulation step is the maximum allowable by the actuators, while the amount extracted over several steps is equal to that prescribed by coefficient of restitution and the contact geometry in Equation (4.40). Furthermore, the kinematic correspondence between the virtual tool and the device is preserved. Simulated hand trajectories obtained by using these techniques (assuming that the maximum wrench capability of the device is Fn it m  = (15N  15N l N m ) ) are presented in Figure 4.18. In this figure, "p/ ;/" T  u  4.5 The performance of the impulse-augmented penalty simulation  61  is the trajectory obtained when the device can fully apply the simulated collision impulses, obtained when collision impulses are saturated on the device,  "p a/ed" S C  "^saturated"  is  is obtained when collision impulses are  scaled in the virtual environment, and "pspread" is obtained when collision impulses are spread over several steps of the simulation. Moreover, "ideal" is the trajectory obtained when the device fully applies the impulses generated by a constraint-based simulation. The hand trajectory degrades as the wrench capabilities of the device decrease. The loss of performance is highest if the interaction forces are saturated on the device. In this case, full collision impulses are applied to the virtual tool which stops abruptly and only limited wrenches are applied to the user's hand, which continues to move. After collision, the user's hand is coordinated with the virtual tool through the position coordination channels of the haptic controller, whose stiffness damping  Hcoord  are much lower than those of the virtual walls, K u wa  and B u. wa  K  c  o  o  r  d  and  Hence, constraint violation  is largest and settling time increases (transient response is poorest). The loss of performance due to limited wrench capabilities of the haptic interface is diminished most by spreading the collision impulses over several steps of the simulation. Therefore, this technique is adopted to overcome device force limitations during haptic manipulation of impulse-augmented penalty virtual worlds. Both scaling and spreading of collision impulses are equivalent to adapting the coefficient of restitution to the device capabilities. As a result of this adaptation, the effective coefficient of restitution may be negative and the post-collision normal contact velocities may be negative (i.e., bodies may move into each other after collision resolution). Nevertheless, the energetic passivity of the proposed collision resolution approach is not affected by an adaptive coefficient of restitution . Hence, the coupled stability of the haptic interaction system 6  is unaffected by the adaptation of the coefficient of restitution. As the force capabilities of the haptic interface decrease, the haptic and visual performance of the proposed approach diminishes. For e = —1, the impulse-augmented penalty-based simulation reduces to a penaltybased simulation. The reduction in haptic performance can be seen by comparing Figures 4.15 and 4.19. In Figure 4.15, it is assumed that the device can fully apply the collision impulses to the user's hand. In Figure 4.19, it is assumed that the device can apply at most Fu i  m t  = (15N  15N l N m ) and collision impulses are spread T  over more time steps when necessary. In this figure, "ideal", " I A P B L M " , and " P B " represent interaction with a constraint-based, an impulse-augmented penalty-based, and a penalty-based environment, respectively, all connected to users through the four channel controller. " C B " represents interaction with a constraint-based U C  virtual world connected to users through the unilateral coupler. The hand trajectories representing interactions with the impulse-augmented penalty-based and penaltiesonly virtual environments are closer to each other in Figure 4.19 than they are in Figure 4.15. The visual performance diminishes correspondingly, because the virtual tool penetrates the constraints deeper when spreading is required than when the device can fully apply the collision impulses. Nevertheless, constraint penetration is smaller in the proposed simulation than in the penalty-based one regardless of the device limitations. Both the haptic and the visual performance of the impulse-augmented penalty-based virtual world is better than that E q u a t i o n (4.40) shows that the post-collision kinetic energy of the contact group is at most equal to its pre-collision kinetic energy for any e £ [—1,1]. 6  62  4.5 The performance of the impulse-augmented penalty simulation  Hand trajectory  -163 -164 _-165 E —166 o ^167 -168 -16!  X  COM t  m m  l  1.5  (a)  Hand rotation  .x 10  03  04  tinfeW'  16  07  08  (b)  Figure 4.18: Simulated hand trajectories during user interaction with the impulse-augmented penalty-based virtual world when the device limitations are ignored and when they are taken into account ("ideal" - full constraint-based collision impulses are applied to users; "p/ ;/" - full collision impulses are applied to users; "p saturated" - collision impulses are saturated on the device; "Pscaied" - collision impulses are scaled in the virtual environment; " P s p r e o d " - collision impulses are spread over several simulation steps when necessary). u  4.5 The performance of the impulse-augmented penalty simulation  63  (b)  Figure 4.19: Simulated hand trajectories obtained when constraint-based ("ideal"), impulse-augmented penalty-based ("IAPB L M " ) , and penalty-based ("PB") interactions are transmitted to users by a four channel controller, and when constraint-based interactions are transmitted to users by a unilateral coupler ( " C B " ) . Collision impulses computed using the proposed method are spread over several simulation steps when necessary. UC  4.5 The performance of the impulse-augmented penalty simulation  Figure 4.20: Simulated forces and torques applied to users by the four channel controller when the virtual world is generated using the constraintbased ("ideal"), the impulse-augmented penalty-based ("IAPB L M " ) and the penalty-based ("PB") methods, and applied by the unilateral coupler when the virtual environment is generated using the constraint-based method ( " C B ) " . Collision impulses computed using the proposed technique are spread over several simulation steps when necessary. UC  64  4.5 The performance of the impulse-augmented penalty simulation  65  Kinetic energy of the user's hand  Figure 4.21: Kinetic energy of the users' hand when constraint-based ("ideal"), impulse-augmented penalty-based ("IAPB L M " ) , and penaltybased ("PB") interactions are applied to users by the four channel controller, and when constraint-based interactions are applied by the unilateral coupler ( " C B " ) . Collision impulses computed using the proposed method are spread over several simulation steps when necessary. UC  of the penalty-based world. In addition, the perceptual advantage obtained by applying abrupt forces to the user's hand upon contact and the increased stability of the interaction are maintained. Note that the wrenches felt by users are larger when users interact with the proposed simulation than when users interact with the penalty-based virtual environment (Figure 4.20). Furthermore, users lose more kinetic energy upon contact with the impulse-augmented penalty virtual constraints than upon contact with the penalty-based constraints (Figure 4.21). The increased contact passivity improves the coupled stability of the haptic interaction.  4.5.3  Experimental validation  In this section, the performance of the impulse-augmented penalty simulation approach is validated experimentally. The controlled interaction depicted in Figure 4.13 is implemented on the planar haptic simulation system described in Chapter 3 using the synchronous communication between the device and the virtual environment (see Figure 3.2(a)). Recall that the user's hand is represented through a constant wrench F / i = (0.32N  —0.4N  0 N m ) , collisions are considered perfectly plastic (e = 0), and the virtual walls have T  stiffness K u = 4000N/m and damping B u — 30N/(m/s). In the experiments, the virtual environment is wa  wa  generated using the penalty and the impulse-augmented penalty approaches . 7  Two implementations of the  i m p l e m e n t a t i o n s of haptic rigid body manipulations of constraint-based virtual environments running synchronously with the control loop are presently not available. T h i s is because constraint-based simulations with guaranteed completion time have been developed only for perfectly plastic collisions between smooth convex objects [6]  4.5 The performance of the impulse-augmented penalty simulation  66  impulse-augmented penalty approach are used: one considering that the device can fully apply the simulated impulses, and one accounting for device limitations. Figures 4.22, 4.23, and 4.24 monitor the device trajectories, the forces applied to users, and the kinetic energy of the device, respectively, when the peg is driven towards the corner, contacts the vertical virtual wall, slides down along it, contacts the horizontal virtual wall, and stops moving. The experimental results in these figures correspond to the cases when the virtual environment is generated using: (i) the impulseaugmented penalty approach and the device fully applies the simulated impulses to users, " I A P B " ; (ii) the impulse-augmented penalty approach and the device applies at most Fu i  m t  = (15N  15N  lNm)  T  to users,  " I A P B L M " ; and (iii) the penalty approach, " P B " . The "ideal" trajectory is generated in M a t l a b ™ through simulating the interaction with a constraint-based virtual environment . 8  Figures 4.22 and 4.24 demonstrate the increased stability of the interaction in the impulse-augmented penalty virtual world compared to the penalty-based one. The peg loses more kinetic energy upon contact, settles into the corner faster, and bounces less when collision impulses are applied to the device. Hence, the impulsive forces applied upon constraint penetration improve the passivity of the virtual world and, thus, the stability of the haptic interaction. Moreover, they imrpove the realism of the interaction, because they generate a device trajectory that is closer to the trajectory imposed by a real rigid corner. As predicted by simulations, the performance of the impulse-augmented penalty-based world decreases when collision impulses are applied over several steps in order to meet the device rendering capabilities. A device that can apply only limited forces and torques dissipates less energy upon impact than that predicted by the chosen coefficient of restitution. However, it dissipates more energy than during penalty-based interaction. Moreover, the perceptual advantage of large forces upon contact is maintained, because the impulsive forces rendered to users are much larger than the penalty-based contact forces despite the fact that they are limited by device capabilities (see Figure 4.23). Note that increased virtual wall damping would also result in larger force transitions upon impact and less bouncing, i.e., more stable contact. However, the virtual damping is limited by the physical damping, the virtual wall stiffness, and the simulation step during 1 D O F interaction with a virtual wall [41]: the maximum allowable virtual damping decreases as the physical damping in the haptic interface decreases and as the wall stiffness and the simulation step increase. Hence, only limited improvements in the perceived rigidity of the virtual contacts could be achieved through increasing the virtual damping if the virtual wall stiffness is large (as needed for a convincingly rigid resting contacts) and the physical damping is small (as needed for imperceptible device dynamics). On the other hand, the impulsive forces provide a physically-based technique for enhancing the realism of the interaction that increases the stability of the interaction (without increasing the kinetic energy of the simulated environment) and is limited only by the device capabilities. 8  N o t e that, due to friction in the interface, the "ideal" and the experimental trajectories have slightly different initial conditions.  A s a result, the virtual peg contacts the horizontal wall earlier in the simulation than in the experiments.  4.5 The performance of the impulse-augmented penalty simulation  Device trajectory  -160  C0M I™]  X  (a)  0.04  Device rotation  Figure 4.22: Device trajectories when the peg is pushed into the corner by a controller force. The virtual environment is generated using the impulse-augmented penalty-based and the penalty approaches ("IAPB" - impulse-augmented penalty-based world, full collision impulses applied to the device; "IAPB L M " - impulse-augmented penalty-based world, limited impulses applied to the device; " P B " - penalty world). The "ideal" trajectory is generated using a constraint-based simulation.  67  4.5 The performance of the impulse-augmented penalty simulation  68  Environment force along y-axis  Environment force along x-axis  iApp L M  i;  ideal  IAPB  0.35  time[s]  0.3  0.4  ^  I 0.35  timefs]  (a)  0.4  (b)  Environment torque  0.25  0.3  0.35jj r lp.4 me  s  0.45  0.5  (c)  Figure 4.23: Environment wrenches applied to the device when the peg is pushed into the corner by a controlled force. The virtual environment is generated using the impulse-augmented penalty-based and the penalty approaches ("IAPB" - impulse-augmented penaltybased world, full collision impulses applied to the device; " I A P B L M " - impulse-augmented penalty-based world, limited impulses applied to the device: " P B " - penalty world). The "ideal" wrenches are generated using a constraint-based simulation.  0.45  0.5  4.6 Discussion  69  Figure 4.24: Kinetic energy of the device while it is pushed towards the corner by a controlled force. The virtual environment is generated using the impulse-augmented penalty-based and the penalty approaches ("IAPB" - impulse-augmented penalty-based world, full collision impulses applied to the device; "IAPB L M " - impulseaugmented penalty-based world, limited impulses applied to the device; " P B " - penalty world). The "ideal" wrenches are generated using a constraint-based simulation.  4.6  Discussion  The impulse-augmented penalty simulation approach introduced in this chapter has been designed to enable users to feel large forces upon contact with virtual objects and, in doing so, to increase the perceived rigidity of the virtual world [90,125]. Compared to existing simulation paradigms, the approach allows haptic rendering of rigid body collisions, hence a more realistic representation of virtual contacts to users.  Simulated and  experimental interactions show that the selected friction model renders the stick-slip phenomenon similar to other models of dry friction used in haptics and that the selected model may potentially be used for changing the feel of sticking contacts. The simultaneous collision resolution never increases the kinetic energy of the virtual environment, thus maintaining the physical accuracy of the simulation and improving the stability of the virtual contacts.  Moreover, it implicitly eliminates the overdetermined constraints, thereby increasing simulation  efficiency.  Inadequate force capabilities of the haptic interface degrade the performance of the approach.  Performance degradation is limited through a novel technique that adapts the coefficient of restitution to the device capabilities without increasing the kinetic energy of the virtual environment. The new technique spreads the collision impulses over several steps of the simulation when the haptic interface cannot dissipate the kinetic energy prescribed by the desired coefficient of restitution during one simulation step.  Chapter  5  M a n i p u l a t i o n of serial linkages The previous chapter has proposed a simulation approach that computes impulsive forces upon contact and penalty-based and friction forces during contact. Virtual world dynamics have been formulated in configuration space, such that distinguishing between contacts of virtual objects and contacts of virtual linkages has not been necessary. This chapter discusses realistic kinesthetic operation of virtual linkages. Realistic haptic manipulation of virtual linkages depends on: (i) the physical accuracy of the linkage simulation; (ii) the user's perception of the linkage contacts; and (iii) the user's perception of the linkage dynamics and of the topological constraints. Since linkage dynamics are challenging to simulate at the haptic rates, prior work focused on developing both application-specific [114] and general-purpose [126-128] techniques that address the first two factors. In particular, efficient dynamics are the main concern in [114,126,128] and accurate collision/contact models are developed in [127]. In all cases, only penalty-like forces are applied to users. Moreover, operation of linkages only from links with redundant degrees of freedom is allowed. In this chapter, simulation and control techniques are introduced that enable realistic operation of linkages from any user-selected link and through singularities. The chapter starts by discussing the control architecture employed to render convincing force interactions between a virtual linkage and a rigid multibody virtual environment.  A coordinate-invariant representation of contact forces/impulses at the user's hand is then  introduced. The control methods that allow users to convincingly perceive the inertia and the topology of the virtual linkage follow. Combined with the impulse-augmented simulation approach presented in the previous chapter, these techniques permit realistic and unrestricted haptic manipulation of serial linkages. Throughout this chapter, it is assumed that users manipulate a virtual linkage, i.e., virtual tool and virtual linkage are used interchangeably. The schematic of haptic manipulation of virtual linkages within an impulse-augmented virtual world is shown in Figure 5.1. For simplicity, the position coordination signals used by the haptic controller are not represented.  In this figure, F  c o n s t r  ,  Fi  n e r t  i , a  Fh,  and  F  e n v  are the wrenches due to the  topological constraints, the linkage inertia, the user, and the contact forces/impulses that act on the virtual tool, respectively, all represented at the user's hand. fj., and f-j are the contact forces/impulses that act between the virtual linkage and other virtual objects, n  t c  represents the configuration space directions of constraint 70  5.1 Control architecture for realistic linkage manipulation  71  imposed on the user's motion by the linkage topology, A^" is the inverse of the operational space inertia of 1  the linkage [80] at the user's hand, and T  env  is the configuration space torque due to the contacts between the  virtual tool and the virtual environment. F  mmi  .+F„  Device Inverse of the linkage inertia at user's hand A controller Impulse-augmented penalty simulation Direction of topological h  • F..fl.  F,  Haptic controller -3  constraint n„ Coordinate invariant F„„ at user's hand  Figure 5.1: Realistic haptic manipulation of virtual linkages in an impulseaugmented penalty virtual world.  5.1  Control architecture for realistic linkage manipulation  Besides the linkage simulation, both the device controller and the haptic controller contribute to realistic force interaction within the virtual environment during linkage manipulation. The device controller [136] uses impedance control [68] to apply the inertia of the linkage to the user's hand. In particular, it changes the impedance of the haptic device to match the impedance of the virtual tool. The haptic controller [136] uses a four channel architecture to apply the contact forces/impulses that act on the virtual tool to the user's hand (see discussion in Section 4.1). A 1 D O F mechanical equivalent of this control architecture (excluding the communication delay) is depicted in Figure 5.2. In this respectively,  fh,  fenv,  figure, and  and Z  vt  are the impedance of the haptic device' and of the virtual tool,  are the hand, the environment, and the coordination force , respectively, and 1  f  pc  S scales the two force channels of the haptic controller (from the device to the virtual environment and viceversa). If the device impedance equals the virtual tool impedance, Zd = Z , and forces are not scaled, S = 1, vt  the haptic controller is transparent. Users perceive the inertia and the topology of the linkage through the virtual tool impedance Z . They perceive the motion constraints imposed by other virtual objects through vt  the environment force  f  e n v  .  Hence, the transparency of the haptic controller and the accurate reflection of the  virtual tool impedance through device control are important for realistic haptic manipulation of linkages. Although not apparent in the 1 D O F mechanical equivalent of the control architecture, the representation at the user's hand of the contact forces/impulses acting on the linkage also affects the realism of the interaction. This representation is discussed in the following section. Haptic rendering of the linkage inertia and of the motion constraints imposed on users by the linkage topology are presented in later sections. ' T h e coordination force is applied by the position coordination channels of the haptic controller. These channels are shown in Figure 5.2 as a spring-damper connection between the device a n d the virtual tool.  5.2 Linkage contacts at user's hand  72  tnp  Figure 5.2: 1 D O F mechanical equivalent of the control architecture.  5.2  Linkage contacts at user's hand  The rigid multibody simulation approach developed in Chapter 4 evolves the virtual environment dynamics in configuration space. In this space, virtual objects are indistinguishable from virtual linkages. However, users interact with the virtual environment in Cartesian (operational [80]) space through a haptic controller that applies the simulated forces/impulses to their hand (see Figure 5.1). Hence, the contact forces/impulses acting on the virtual tool must be suitably represented at the user's hand during virtual linkage manipulation . A 2  coordinate invariant representation of these contact forces/impulses at the user-selected link is proposed in this section. The user-applied wrench Fh is mapped to the hand configuration torque  by the transpose of the Jacobian  of the virtual linkage computed at the user's hand, Jl, herein called the "hand Jacobian": (5.1) Hence, j £ maps the space F of wrenches at the user's hand to the space T of configuration torques. Conversely, the space of configuration torques must be mapped to the space of wrenches at the user's hand in order to allow the haptic controller to apply the contact forces/impulses that act on the virtual linkage to the user. When Jj[ is invertible, the required mapping is provided by the transpose of the inverse of the virtual tool Jacobian at the user's hand: (5.2)  env •  In Equation (5.2), T  env  is the configuration space representation of the contact forces/impulses acting on the  virtual tool (i.e., the environment configuration torque), and F  env  is its wrench representation at the user's  hand (i.e., the environment wrench). According to the derivations in Chapter 4:  (5.3) i=l during resting contact, and: (5.4) i=i Cartesian (operational) space coincides with configuration space during manipulation of virtual objects and the m a p p i n g of contact forces/impulses to configuration space through the contact Jacobians J j suffices. 2  5.2 Linkage contacts at user's hand  73  during colliding contact. In Equations (5.3) and (5.4), c  is the number of contacts between the virtual tool  vt  and the virtual environment (assuming a suitable numbering of the contacts of the contact group to which the virtual tool belongs), and jT  = [Jfni  ...  Jfn,  ...  J^ n i  ] . r  C u (  During unrestricted haptic manipulation of arbitrary linkages,  may not be invertible. This is because:  (i) the linkage may have any number of degrees of freedom (in particular, more or less than 6DOF); (ii) the user may choose to operate it from any link (from the one proximal to the fixed base to the one distal to the base); and (iii) the user may move the linkage through singularities. In all these situations, the configuration torques must be mapped to hand wrenches such that the results are physically significant, i.e., gauge invariant under rigid body transformations [46]. The mapping of hand wrenches to simulated configuration torques given in Equation (5.1) and the mapping of simulated configuration velocities q to simulated hand twists x/^,: *fc„, = J ^ q  (5.5)  form a dual system [46] (see Appendix B for definition). In Equation (5.5), x.h = [ v £ vl  c u ] , where v T  T  h  € TZ  3  and iv S TZ are the linear velocity of the user-selected operational point on the virtual tool and the angular 3  velocity of the virtual link held by the user, respectively. Therefore, as discussed in detail in [46] and briefly summarized in Appendix B , a mapping from T to F is provided by the transpose of the weighted generalized inverse of  [46]: J* =MQ C 1  t  (CMg'C ') 7  1  (F Mi/F) T  - 1  F M .  (5.6)  T  v  In Equation (5.6), J * is a weighted generalized inverse of the hand Jacobian,  = F C is a full-rank decompo-  sition (factorization) [20] of the hand Jacobian J ^ , M Q is the metric on the space Q of configuration velocities, 3  and M y is the metric on the space V of hand twists. Coordinate invariance of the mapping J * is ensured through a suitable choice of metrics on the spaces Q and V [46]. In [46] it is shown that kinetic energy metrics are gauge invariant. In this work, the metric on Q is provided by the configuration space mass matrix of the virtual tool D „ f The metric on V is provided by the mass matrix of the link held by the user M / , j , computed in the link's coordinate frame. Hence, the coordinate invariant weighted generalized inverse of the Jacobian of the virtual tool computed at the user's hand is given by:  J# = D - ^ C  1  (CD- C ) 1  T  (F M F)  _ 1  T  t  U  _ 1  FM . T  hl  (5.7)  The properties of jf are presented in Appendix B . The algorithm proposed for computing a full rank factorization of of  J-  is detailed in the following. The algorithm uses the Singular Value Decomposition (SVD)  h  J  3  A f u l l - r a n k d e c o m p o s i t i o n [46]  k«xd  = F C satisfies  = U6x6S xdVj 6  x d  .  rank (J^) = row-rank ( C ) = column-rank ( F ) .  (5.8)  5.3 Haptic rendering of linkage inertia  74  For a hand Jacobian with rank (J/j) = r, only the first r diagonal elements of the matrix X are non-zero and:  rxr  °rx(d-r)  0(6-r)xr  0(6-r)x(d-r)  S  •'6xd  where I  r  x  r  Irxr  0(6_ ) r  is the r-dimensional unit matrix and  (5.9)  Oj-x^-r)  0(6-r)xr  r  denotes a matrix of zeros with the given dimensions.  X  A full rank factorization of J/j is then given by: F  = USx  (5.10)  V and the coordinate invariant inverse of the hand Jacobian for the chosen set of metrics is: J*  = D7 VSnS2V D- V2^)" (Sj'U M US )~ S?'U M . 1  T  1  t  1  T  t  1  w  r  1  Hence, the environment configuration torque is mapped to the environment wrench F  w  env  env — J/j, T  5.3  e n v  (5.11)  according to:  .  (5.12)  Haptic rendering of linkage inertia  In Cartesian (operational [80]) space, both the inertia and the topological constraints are embedded in the operational space inertia of the virtual linkage computed at the user's hand [80], A^: A =(j D^ jJ')" . 1  h  (5.13)  1  h  Therefore, A/, is incorporated in the desired impedance of the virtual tool Z . Specifically, the device controller vt  changes the dynamics of the device to match the desired dynamics of the virtual tool [136]:  Ax h  + bx + kx  h  d  h  d  h  =F + F h  e  nv H" Fp . C  In Equation (5.14), A/, € 1Z * , 6  S 7Z , and  6  €7c  6x6  and stiffness of the virtual tool, respectively, F  6 x 6  (5.14)  are the desired (i.e., simulated) inertia, damping,  € Tl is the position coordination wrench between the device 6  p c  and the virtual tool (due to the generalized spring-damper shown acting between them in Figure 5.2), and X f t 6 TZ , x/j e TZ , and x^ € TZ are the desired body acceleration, velocity, and position of the device at the 6  6  6  user-selected link, respectively. Equation (5.14) is equivalent to:  x =A h  1 h  (F + F h  env  +F - bx - kx) pc  d  h  d  h  (5.15)  5.3 Haptic rendering of linkage inertia  75  This form is advantageous because the inverse of the operational space inertia of the virtual tool at the user's hand can be computed regardless of the rank of the Jacobian:  (5.16)  Aft = J f c D ^ j J . 1  In contrast, A/i can be computed only when the Jacobian is full row-rank. When  is not full row-rank, A ^  1  becomes infinite along certain directions of the operational space and AJ^ drops rank. 1  Rank deficiency of A^~ indicates that the topology of the virtual tool restricts the simulated instantaneous 1  motion of the user. This may happen when the user holds a link with fewer than 6 DOFs during spatial interaction and fewer than 3 DOFs during planar interaction or when they operate the virtual linkage through a singularity. To illustrate how the virtual motion constraints imposed on users by the linkage topology are represented in the rank deficiency of A^" , consider the example manipulations of the planar linkage depicted 1  in Figure 5.3. In the first example, the user holds the linkage from the centre of mass of the middle link in the position shown in Figure 5.3(a). The held link has insufficient degrees of freedom to allow arbitrary position and orientation at the user-selected operational point. Hence, the linkage topology instantaneously constrains the user's motion according to: Xft = Jhq,  (5.17)  where the hand Jacobian is: -h sin(gi) - l  C2  h  cos (qi) + l  C2  sin (91 + q )  -l sin(q  2  cos (<7i + q ) 2  + q) 2  0  Z cos (qi + q )  0  C2  1  2  C2  1  1  (5.18)  0  Equation (5.18) shows that, due to the topological constraint, rank (J/j) = 2. From Equation (5.16), it follows that rank  (A^" ) = 2. In the second example, the user holds the linkage from the centre of mass of the distal 1  link in the position shown in Figure 5.3(b). In the example position, the Jacobian:  -Zi sin {qi) - l sin (gi) - l 2  C3  Zi cos (gi) + Z cos (qi) + l  C3  2  sin (gi + 93)  ~h  cos (qi + q )  l  3  2  sin (gi) - l  C3  sin (gi + q ) 3  cos (g ) + Z cos (gi + q ) x  1  3  C3  -l  C3  l  C3  sin (qi + q ) 3  cos(<?i + q ) 3  1 (5.19)  is singular, rank (J/,) = 2. Consequently, rank (A7J ) = 2. In other words, A^" drops rank when the linkage 1  1  restricts the user's motion. Since A ^ can be computed regardless of the rank of the hand Jacobian, Equation (5.15) is used for the 1  impedance control of the device. This equation computes the desired body acceleration at the user's hand. The dynamics of the haptic interface are given by:  M Xft + CdXh = F + u, d  h  (5.20)  5.4 Haptic rendering of linkage topology  76  (a) M a n i p u l a t i o n f r o m a l i n k p r o x i m a l t o t h e b a s e .  (b) M a n i p u l a t i o n t h r o u g h a s i n g u l a r i t y .  Figure 5.3: Example manipulations restricted by the topology of the virtual linkage.  where  and C  d  are the inertia and Christoffel matrices of the device and u is the control signal [136]. By  combining Equations (5.15) and (5.20), the impedance control law for the interface is obtained as: u =  {M. Kl - I) F l  d  + MdA^ (F „ + F ) 1  f c  en  pc  + (C - MdA^bd) x - M A ^ ^ x , , . d  ft  d  (5.21)  Equation (5.21) applies the inertia and the topology of the virtual tool to users at the acceleration level through A ^ . Due to numerical drift and limited device stiffness, Equation (5.21) is not sufficient for realistic 1  haptic rendering of the topological constraints. A technique that enables users to convincingly perceive these constraints is presented in the following section.  5.4  Haptic rendering of linkage topology  An analysis of the virtual tool dynamics suggests how the control in Equation (5.21) can be augmented to enable users to accurately perceive linkage topology. Consider that, at the user-selected operational point, the environment wrench is F  and the inverse of the virtual tool operational space inertia is A ^ . This inverse 1  env  maps the environment wrench to the body acceleration at the user's hand: x  = Aft F 1  k  eni)  .  In Equation (5.22), x/, is an element of the space of body accelerations at the user's hand, A; i.e., x^ = (v;[  (5.22) cu ) T  where v/, € TZ and cb £ TZ are the linear and the angular acceleration of the user's hand, respectively. When 3  3  A ^ is full rank, the environment wrench produces body acceleration of the user's hand along all directions of 1  A. In other words, the virtual tool behaves as an inertia along all directions of A. On the control side, Equation (5.21) applies this inertia to users by controlling their acceleration to the value given in Equation (5.22). When A^~ is rank-deficient, x/, lies in a subspace of A. Then, environment wrenches exist that lie in the null 1  space of A ^ , M (A^" ) (hereafter called the space of wrenches of topological constraint). These wrenches are 1  1  directly opposed by the structural stiffness of the virtual linkage and have no effect on the body acceleration at the user's hand. Hence, at the user-selected operational point, the linkage behaves as an infinite stiffness along  77  5.4 Haptic rendering of linkage topology  the directions of topological constraint and behaves as an inertia along the directions of topological freedom of motion . On the control side, the impedance control in Equation (5.21) is unsuitable for rendering the infinite 4  structural stiffness of the virtual tool along the directions of topological constraint. Realistic perception of the topological constraints is enabled through augmenting the control provided by Equation (5.21) with penalties applied along the singular directions of A ^ " , as shown in Figure 5.4. In this 1  figure, vt is the user's position and orientation on the virtual tool (i.e., in the simulation), while h is their real position and orientation. For simplicity, only a one dimensional null space of A^" is depicted. 1  Figure 5.4: Penalty wrench constraining users to the configuration manifold of the virtual linkage that they manipulate.  The directions of topological constraint n  are provided by the null space of AJ^ . Note that the maximum  tc  number of such directions is five during spatial rigid body (6 DOFs) manipulation and two during planar rigid body (3 DOFs) manipulation. The constraint position and velocity are provided by the body position and velocity x^^, of the user-selected operational point on the virtual tool. Then, given the stiffness k  tc  and  damping b of the topological constraints, the impedance control law in Equation (5.21) is modified according tc  to:  u  ^Ji.  =  ( M d A " - I) F + M d A ^ " ( F „ + F ) + (C„ - M ^ A " ^ ) x -  hvi  {kt nj ,i ( h„ x  c  c  1  t  d  tc  ci  - x )) h  ] T {k nj »=i tc  c4  p c  ( x „ , - x ) + b nj h  fe  tc>  1  e n  tc  ci  ft  ( x „ , - x )) n , h  where n is the number of topological constraints, i.e., nt = 6 — tc  hvt  1  h  +M  ft  (k  n }j + C x  M  d  +^  - x ) + b nf  =  c  h  d  -F  h  h  =  M A^k x d  d  h  (5.23)  tC;i  rank(A~ ). 1  Note that no additional numerical effort (beyond computing the coordinate invariant representation of the contact forces/impulses at the user's hand) is necessary to derive the directions of topological constraint n j. Ci  4  T h e directions of topological freedom of motion form a basis on ^(J/,), the range space of J ^ .  5.4 Haptic rendering of linkage topology  78  This is because the singular directions of A ^ are the same as the singular directions of JT. The proof is 1  given in Appendix C. Physically, the proof can be understood by considering the actions of J j and of A^" . 1  In particular, the user-applied torques balanced by the structural stiffness of the linkage are filtered out from the virtual tool configuration space dynamics by Jl.  A^" filters out the same torques from operational space. 1  Hence, the two mappings have equal null spaces and the directions of topological constraint are a byproduct of the full rank factorization of JhThe ability of the control law in Equation (5.23) to enforce linkage topology is validated through simulations in the following section and through experiments in Section 5.4.2.  5.4.1  Simulations  In the simulations, users manipulate a three links planar linkage with rotational joints. Its parameters are given in Table 5.1. Initially, the virtual linkage is at rest in the position q = ( f rad  frad  —|rad)  T  (see  Figure 5.5). The user applies a constant force fh = I N along the x direction. The device is also at rest, but its position differs from that of the user-selected operational point by 5 mm along the x direction and by 3 mm along the y direction. Table 5.1: Parameters of the three links planar virtual linkage operated by the user in the simulations. L i n k l e n g t h (m)  L i n k mass (kg)  Link inertia (kgm )  h = 0.042 m  mi = 3 kg m = 3 kg m = 1 kg  h = 0.015 kgm h = 0.015 kgm h = 0.005 kgm  l = 0.042 m 2  / = 0.03 m 3  (a) M a n i p u l a t i o n  2  3  from third link.  (b)  2  Manipulation  from  2  2  2  second  link.  Figure 5.5: Simulated manipulations of a planar virtual linkage. Initial linkage position is shown in black. Linkage positions during manipulation are shown in grey.  The Simulink™ diagram of the simulated manipulation is shown in Figure 5.6. In the diagram, the virtual environment is represented by the configuration space dynamics of the linkage. The haptic device controlled  5.4 Haptic rendering of linkage topology  79  according to Equation (5.21) is represented by the operational space dynamics of the linkage. The haptic device controlled according to Equation (5.23) is represented by the operational space dynamics of the linkage augmented with the control signal:  Uadd  {k n[c,i ( ft»< - h) + btcnj {x - x )) n x  = M ^  x  tc  d  c<i  hvt  h  t C i i  ,  where M d is considered constant throughout the workspace and equal to unity. The coordination between the device and the virtual tool has stiffness K (70 N/(m/s) straints are k  70 Nm/(rad/s) tc  p c  = (100 N / m  100 N / m  0.5 N m / r a d )  T  and damping B  p c  =  0.375 N / ( m / s ) ) , while the stiffness and damping of the topological conr s  = 10 N / k g / m and bt = 7 N/kg/(m/s), respectively. c  1 >ciiii. under impi'iLinci: rii:urol  Figure 5.6: Simulink diagram of the haptic manipulation of a planar virtual linkage.  Figure 5.7 depicts the results for the first simulated interaction, in which the user manipulates the virtual tool from the C O M of the distal link, as shown in Figure 5.5(a). The user's hand trajectories on the device and in the simulation are shown in Figure 5.7(a) for the case when topology is imposed at the acceleration level. They are shown in Figure 5.7(b) for the case when user's departure from the constraint manifold is penalized according to Equation (5.24). Since the user holds a link with 3 DOFs, the topological penalties are applied only intermittenly, when they user moves through a singularity. Note that the device drifts from the virtual tool when linkage constraints are imposed at the acceleration level. In contrast, the drift is substantially reduced through penalizing violations of the virtual topology. The results for a second simulated interaction are shown in Figure 5.8. In this interaction, the user manip5  T h e y are chosen to match the values implemented on the haptic device used in the experiments [136].  80  5.4 Haptic rendering of linkage topology  Hand trajectory  (a) T o p o l o g y i m p o s e d a t t h e a c c e l e r a t i o n l e v e l .  Hand trajectory  (b) T o p o l o g y i m p o s e d t h r o u g h p e n a l t i e s .  Figure 5.7: Simulated planar manipulation of a three links virtual tool held from the COM of the distal link. ulates the virtual tool from the COM of the middle link, as shown in Figure 5.5(b). This link has only 2 DOFs. Hence, the virtual topology constrains user's motion throughout the interaction. Once more, the trajectories in Figures 5.8(a) and 5.8(b) illustrate that the drift between the user and the virtual tool is significant unless user's departure from the configuration manifold is appropriately penalized through device control. Hand trajectory  Hand trajectory  x rcml  x rcml  (a) T o p o l o g y i m p o s e d a t t h e a c c e l e r a t i o n l e v e l .  (b) T o p o l o g y i m p o s e d t h r o u g h p e n a l t i e s .  Figure 5.8: Simulated planar manipulation of a three links virtual tool held from the COM of the second link.  5.4 Haptic rendering of linkage topology  5.4.2  81  Experiments  The control in Equation (5.23) imposes linkage topology on users through employing the virtual tool directions of topological constraint, n ; . In addition, the physically motivated mapping of the virtual dynamics from tCj  configuration space to operational space involves the computation of the null space of the linkage Jacobian at the user's hand, J . Both the computation of the directions of topological constraint and the computation of k  the null space are based on the S V D of J/,, which may be numerically expensive for the speed requirements of the haptic control loop. Suitable approximations of these directions, n ,i, and of the virtual linkage dynamics tc  are derived in the following chapter, in Section 6.3. Herein, they are used to illustrate the ability of the control in Equation (5.23) to limit users' motion as required by the virtual linkage that the users operate. The controlled experiments presented in this section mimic linkage manipulation from various user-selected links. These experiments compare two types of trajectories.  Trajectories obtained by enforcing topological  constraints through penalizing the user's departure from the configuration manifold are contrasted to trajectories obtained by enforcing virtual topology through the inverse of the operational space inertia of the linkage (i.e., through controlling to zero user's acceleration along directions orthogonal to the configuration manifold). In the experiments, a constant wrench F^ = (0.4N ON 0 N m ) represents the user . In the first experiT  6  ment, the wrench is applied to the C O M of the distal link of the linkage shown in Figure 5.9 (dimensional and inertial properties were given in Table 5.1). In the second experiment, the wrench is applied to the C O M of the second link of the same linkage. In both experiments, the linkage is initially at rest, in the configuration space position shown in Figure 5.9 (i.e., qo = (Orad  f rad  — ^ r a d ) ) . Asynchronous communication is used. The T  linkage dynamics are approximated by: D ^ q = JlF ,  (5.25)  h  i.e., the linkage moves freely under the controlled hand wrench. In Equation (5.25), D„t is the approximate inertia matrix of the virtual tool and is computed as described in Section 6.3.1. Furthermore, the control in Equation (5.24) is approximated by:  Uadd  = M  d  ^2  { tcS-J ,i k  c  ( fc„« x  _  x  /i) +  tc^f ,i (xft„, - x )) n  b  c  ft  t C i i  ,  i=l  where nt j approximates i-ih topological constraint direction, as described in the following chapter, in SecCi  tion 6.3.2. The gains of the control in Equation (5.26) are k  tc  = 200N/(kg-m) and b = 50N/(kg-(m/s)). tc  The user's hand trajectories on the haptic device ("HD") and in the simulated virtual environment ("VE") plotted in Figure 5.10 correspond to manipulation from the distal link. The trajectories shown in Figure 5.11 represent manipulation from the second link. Figures 5.10(a) and 5.11(a) demonstrate that the control in Equation (5.24) effectively imposes the topological constraints on users regardless of the link they operate. The user's hand trajectory on the device follows the simulated user's hand trajectory in the virtual environment s  F o r the impedance type haptic interface employed in these experiments, a constant wrench represents a worst case approxi-  mation for stability because it eliminates the adaptive d a m p i n g of the user's hand  [64].  (5.26)  5.4 Haptic rendering of linkage topology  82  Base (ground) join;  Figure 5.9: Testbed virtual environment for controlled linkage manipulation from any user-selected link.  within the steady state error due to the limited stiffness of the penalties applied to users. On the other hand, Figures 5.10(b) and 5.11(b) show that the drift between the user's hand trajectories on the haptic device and in the virtual environment simulation is large when virtual topology is enforced only at the acceleration level. The impedance controller alone cannot enforce the virtual topology. The user's hand freely moves away from the configuration manifold until it reaches the mechanical boundary of the workspace of the haptic interface. In Figures 5.10(b) and 5.11(b), the sudden changes occuring in the trajectories of the user's hand on the haptic device after approximately 0.7s represent hard constraints due to device workspace limitations.  time [sec]  time [sec]  (a) Topology enforced through penalties.  (b) Topology enforced at acceleration level.  Figure 5.10: Haptic manipulation from the distal link of the virtual linkage shown in Figure 5.9. Users apply a constant wrench F = (0.4N ON 0 N m ) . T  h  5.4 Haptic rendering of linkage topology  83  Figure 5.11: Haptic manipulation from the second link of the virtual linkage shown in Figure 5.9. Users apply a constant wrench F = (0.4N ON 0 N m ) . T  h  The experiments described in this section have demonstrated realistic haptic manipulation of serial linkages from various user-selected links. In these experiments, motion constraints due to linkage topology have been imposed on users through penalties applied along the linkage directions of infinite structural stiffness (i.e., orthogonal to the linkage configuration manifold). A virtual coupler controller [3,42] would also allow users to perceive the linkage topology. However, the control proposed in this thesis allows larger gains to be used for enforcing the virtual topology. This is because the proposed control is applied only when users depart from the configuration manifold and its gains are limited only by the device itself and by the sampling rate. On the other hand, the virtual coupler acts throughout the interaction and its gains must be chosen so as to achieve a compromise between free and constrained motion. While large gains are necessary during constrained motion, low gains are desirable during free motion. Since the proposed control enforces the linkage topology through larger gains than the virtual coupler, it enables users to perceive stiffer topological constraints. Hence, it allows more transparent interaction.  This chapter has proposed simulation and control techniques that enable convincing and unrestricted haptic manipulation of virtual linkages from any user-selected link. Physically motivated mappings of the virtual dynamics from configuration space to operational space have been derived. A control penalizing users' departure from the configuration manifold of the linkage has been proposed for limiting users' motion as required by the virtual linkage that they operate. Simulations and experiments performed using a planar haptic interaction system have demonstrated that the techniques described in this section can be used to enable unrestricted operation of virtual linkages. In the experiments, the real time performance of these techniques has been  5.4 Haptic rendering of linkage topology  84  guaranteed through efficiently implementating them in a local model of interaction. This implementation is detailed in the following chapter.  Chapter  6  Efficient haptic rendering of rigid b o d y motion w i t h constraints This chapter presents a local model of interaction that is suitable for adding realistic forces and torques to the manipulation of virtual tools within rigid multibody virtual environments.  The local model efficiently  implements the earlier developed simulation and control techniques that facilitate transparent haptic rigid body interaction within virtual worlds. After presenting a model synopsis, the chapter details the geometric and the dynamic approximations that enable the proposed local model to achieve haptic speed requirements. The chapter ends with an experimental comparison between the proposed local model of rigid body interaction and prior local models. The schematic of haptic virtual tool manipulation using the proposed local model is shown in Figure 6.1. In this figure, dashed lines indicate low bandwidth, asynchronous communication between the virtual environment and the local model of interaction. Solid lines indicate high bandwidth, synchronous communication between the local model and the device. For simplicity, position coordination signals are not represented. Local model of rigid body interaction  Controller  Virtual D i r e c t i o n o f the topological constraint n  lc  (at update)  Inverse o f the linkage inertia at the user's hand A ' (at update* h  Geometry Proxy state (position, velocity)  F  +F.  env u Figure 6.1: Realistic haptic manipulation of virtual tools using the local model of interaction proposed in this thesis.  85  environment  6.1 Synopsis of the local model  6.1  86  Synopsis of the local model  To meet the need of haptic devices for physically motivated forces and torques provided at high fixed frequencies, the local model of interaction implements a less complex simulation. The local model approximates the interaction between the virtual tool and the virtual environment through the interaction between the virtual tool and nearby objects only (see Figure 6.2). A proxy of the virtual tool and constraints imposed on the motion of the virtual tool by nearby objects comprise the model. The quality of the approximation is maintained by updating the local model at each step of the virtual environment simulation. Coordination between the proxy and the haptic device is achieved using force and state (i.e., position and velocity) information, while coordination between the proxy and the virtual tool in the virtual environment is achieved using proxy state (i.e., position and velocity) information.  Haptic Ijyrenches, motions! Microcontroller device .. Local model . hundreds of H z  ^geometry, topology proxy state  Host VE  tens of H z  Figure 6.2: Communication between the virtual environment (VE), the local model of interaction, and the haptic device.  Four features distinguish this local model of interaction from existing models. First, the model considers not only the virtual objects already in contact with the virtual tool, but also nearby virtual objects. As shown in Section 6.2.3, prediction of virtual tool contacts with nearby virtual objects effectively alleviates the delay in updating local geometry when users operate the virtual tool through small clearances.  Second, the local  model computes impulsive forces upon contact in addition to penalty and friction forces during contact. By dissipating the kinetic energy of the user's hand as prescribed by the coefficient of restitution and allowed by the device capabilites, the impulsive forces increase the stability of the virtual contacts. By enforcing user's hand accelerations that are much larger upon contact than the user's hand accelerations imposed by the penalty forces during contact, the impulsive forces enhance the perceived rigidity of the virtual contacts. Third, the model includes a dynamic proxy of the virtual tool in addition to local geometry. The dynamic proxy allows a four channel haptic controller [136] to be used that both transmits wrenches between the device and the proxy and coordinates positions between them.  Compared to the virtual coupler [3,30], the four channel  controller transmits the impulsive wrenches computed by the local model of interaction to the user's hand. Hence, it enables users to feel collisions upon impact and thereby increases the perceived rigidity of the virtual  6.2 Local geometry  87  objects [90,125]. Compared to directly coupling an impedance device to an admittance simulation [21], the four channel controller allows stiffer contacts to be stably rendered to users. This is demonstrated through experiments in Section 6.4.3. Lastly, the new local model can be used to add realistic haptic feedback to interactive virtual environments generated using any commericial simulation package. This is because the model makes no assumptions about the data structures or the collision detection and the dynamic response algorithms used by the virtual environment simulation. This feature is illustrated in Section 6.4 by using the proposed local model of interaction to allow haptic interaction within a testbed virtual world generated using V o r t e x ™ , a physics engine developed by CMLabs Simulations Inc (www.criticalmasslabs.com). The features of the local model of interaction are presented in the following sections, starting with the local geometry in Section 6.2. The local dynamics are detailed in Section 6.3.  6.2  Local geometry  The forces and torques that provide a convincing haptic experience arise at the contacts between the virtual tool and other virtual objects, collectively referred to as the virtual environment. Hence, the local representation of the virtual world geometry is important for realistic force feedback.  Complete geometry is desirable for  accurate computation of interactions, since it allows all contact transitions of the virtual tool to be resolved locally. Nevertheless, complete geometry cannot be used in the local model due to the complexity and the variability across simulation packages of the data structures and of the collision detection algorithms associated with complex geometric models. Partial (i.e., some subset of the) virtual environment geometry may result in delayed updating of the virtual tool contacts in the local model of the interaction. In turn, this computational delay may cause perceptual artifacts or unstable interaction due to undesirable force discontinuities at model updates. A n example interaction during which undesirable force discontinuities arise at model updates is depicted in Figure 6.3. Successive time steps of the virtual environment are shown in Figures 6.3(a), 6.3(b), and 6.3(c), and the corresponding local model updates are shown in Figures 6.3(d), 6.3(e), and 6.3(f). The example illustrates the case where only virtual environment geometry in contact with the rectangular virtual tool is sent to the local model. Note that the virtual tool contacts arrive late, i.e., with significant penetration, to the local model and cause force discontinuities at updates. In turn, these force discontinuities may destabilize the interaction. The selection of salient geometry to be used in the local model and the choice of a method for transitioning between the old and the new local models are critical for a stable and convincing virtual kinesthetic experience. The techniques employed for selecting appropriate local geometry and for transitioning between successive local models are presented in the following subsection.  6.2.1  Active geometry  While the geometry of the virtual world is important for convincing haptic feedback, its importance varies temporally and spatially. Thus, the instantaneous interactions of the virtual tool are independent of the  6.2 Local geometry  88  (a) V i r t u a l e n v i r o n m e n t at t i m e t.  (b) V i r t u a l e n v i r o n m e n t a t t i m e t+ 1.  (c) V i r t u a l e n v i r o n m e n t a t t i m e t + 2.  (d) L o c a l m o d e l u p d a t e d a t t i m e t.  (e) L o c a l m o d e l u p d a t e d a t t i m e t+1.  (f) L o c a l m o d e l u p d a t e d at t i m e t + 2 .  Figure 6.3: Example haptic interaction during which the computational delay of the virtual environment results in delayed updating of the virtual tool contacts in the local model of interaction.  geometry of virtual objects that are not currently in the same contact group with the virtual tool (i.e., in contact with it either directly or through other objects). For example, the geometry of body C in Figure 6.4 has no influence on the forces presently acting on the virtual tool. Furthermore, the virtual tool interactions are influenced more by the geometry of virtual objects in direct contact with the virtual tool than by the geometry of virtual objects in indirect contact with it. For example, consider that a force  is applied to A and a force  fe is applied to B in Figure 6.4. Both these forces are transmitted to the virtual tool via the contact between the virtual tool and A . Hence, the magnitude of the force applied to the virtual tool depends both on  and  on is- However, the direction of the interaction force between the virtual tool and the virtual environment depends on the geometry of A and is independent of the geometry of B. Thus, virtual environment geometry in direct contact with the virtual tool at simulation updates, called active geometry, is of primary importance for realistic user interaction within the virtual world. The active geometry is chosen initially to approximate the virtual environment in the local model of the interaction.  Figure 6.4: The influence of the geometry of the virtual environment on the interactions of the virtual tool.  6.2 Local geometry  89  Coupling of a device to an arbitrary virtual world is enabled by restricting the active geometry to information generally available in physically-based virtual environments. Typical simulation packages represent rigid bodies as collections of convex polyhedra and compute contacts between pairs of these polyhedra [85]. They approximate rigid body contact through a finite number of vertex-face and/or edge-edge contacts [152]. For each contact, they provide a contact point £  a constraint normal direction n, and a penetration depth.  The constraint normal direction is the direction of motion of the virtual tool contact point that is disallowed by the contact. The penetration depth is equal and opposite to the separation distance s between the bodies; i.e., it is positive when the virtual tool and the virtual environment overlap. The  active geometry consists of  all contacts of the virtual tool with the virtual environment. In the local model of interaction, the active geometry is encapsulated in the local proxy constraints. A local  proxy constraint is defined through the identifiers of the two contacting objects and the constraint geometric, kinematic, and dynamic properties. Typically, the identifiers of the contacting objects are the addresses of their data structures in the simulation. The geometric properties of a local proxy constraint embed the geometric information provided by typical simulation packages for each contact of the virtual tool (see Figure 6.5). As depicted in Figure 6.5(d), these properties consist of: • the  local constraint position £ . This is the position of the face point closest to the vertex when the  virtual tool is in vertex-face contact with the virtual environment (Figure 6.5(a)); it is the position of the vertex when the virtual tool is in face-vertex contact with the virtual environment (Figure 6.5(b)); and, it is the position of the point on the environment edge that is closest to the virtual tool edge when the virtual tool is in edge-edge contact with the virtual environment (Figure 6.5(c)). • the  local constraint normal direction n. This is the constraint normal direction supplied by the simulation.  • the  local contact point P. This is the vertex itself when the virtual tool is in vertex-face contact with  the virtual environment (Figure 6.5(a)); it is the face point closest to the vertex when the virtual tool is in face-vertex contact with the virtual environment (Figure 6.5(b)); and, it is the point on the virtual tool edge closest to the environment edge when the virtual tool is in edge-edge contact with the virtual environment (Figure 6.5(c)). Note that constraints are considered to extend infinitely in the local model of interaction. Moreover, a proxycentric view is adopted in describing the active geometry: constraint normals are directed from the virtual environment to the virtual tool and the contact point is taken on the proxy. It is assumed that both the local constraint position c  and the local constraint normal direction n are  provided by the virtual environment simulation in world coordinates . The local contact point P is in proxy 1  coordinates and is computed locally by:  (6.1) 1 If t h i s i s n o t t h e c a s e , t h e y c a n b e t r a n s f o r m e d i n t o w o r l d c o o r d i n a t e s t h r o u g h a n a d d i t i o n a l r i g i d b o d y t r a n s f o r m a t i o n a n d an additional rotation, respectively.  6.2 Local geometry  90  Virtual tool  Virtual tool  Virtual tool An  Proxy  m  ®  / r Virtual environment (a) V i r t u a l t o o l i n v e r t e x face c o n t a c t w i t h t h e v i r tual environment.  Virtual environment  Virtual environment (b) V i r t u a l t o o l i n f a c e vertex contact w i t h the virtual environment.  (c) V i r t u a l t o o l i n e d g e edge c o n t a c t w i t h the v i r tual environment.  Local constraint (d) L o c a l c o n t a c t .  Figure 6.5: Geometric information provided by typical simulation packages for each type of virtual tool contacts (Figures 6.5(a), 6.5(b), and 6.5(c)) and its representation in the local contact geometry (Figure 6.5(d)).  In Equation (6.1), R is the rotation from the world to the proxy coordinates and c is the position of the ~COM C O M of the proxy. The kinematic properties of a local proxy constraint are provided by the local constraint velocity v  , which is the projection along the constraint normal n of the velocity of the virtual environment  constr  point in contact with the virtual tool v  : Vconstr =n V  •  T  c  (6.2)  The local constraint velocity is used to predict the constraint position at the next simulation update, thus diminishing local geometry discontinuities at model updates. The dynamic properties of a local proxy constraint consist of the constraint stiffness K  ,  contact  the constraint damping  B  c o n t a  ct,  the coefficient of restitution e, and  the coefficients of static and kinetic friction, / J and ^k, respectively. While all local proxy constraints have the s  same stiffness, damping, and coefficient of restitution, they may have different coefficients of friction. The local proxy constraints are the only geometry that the local model is aware of. Hence, only partial virtual tool and virtual environment geometry is available locally. This is illustrated in Figure 6.6, where the local representation of an example virtual world geometry is depicted. Partial geometry makes the proposed approach compatible with simulation packages regardless of the data structures that they use for representing the virtual objects.  Moreover, it simplifies local collision detection, which becomes an iteration through all  local constraints in order to compute separation distances according to:  s= n  T  (c  +R" P - c 1  ).  In Equation (6.3), s is the separation between the proxy and the local constraint of interest and R rotation from proxy to world coordinates.  (6.3)  _ 1  is the  6.2 Local geometry  o  91  Virtual environment  n,  Local constraints (a) E x a m p l e c o n t a c t g e o m e t r y i n t h e v i r tual environment simulation.  (b) L o c a l r e p r e s e n t a t i o n o f t h e v i r t u a l environment geometry shown in F i g ure 6.6(a).  Figure 6.6: Partial virtual environment geometry available in the proposed local model of interaction.  The active geometry shifts the computational delay of the simulation from delay in computing interaction forces to delay in updating local geometry. As a result of faster force computation, the users' haptic experience is improved in two ways: (i) they can manipulate much stiffer objects, due to the high control rate that can be achieved; and (ii), they can feel physical phenomena that rely on fast force transitions, such as collisions and stick-slip friction. However, at model updates, undesirable discontinuities may arise in the locally computed forces that may destabilize the interaction. Two techniques are proposed to avoid such discontinuities: local proxy deformation and the use of e-active local geometry.  These techniques are detailed in the following  Sections 6.2.2 and 6.2.3.  6.2.2  Local proxy deformation  A n update of the local model may result in: (i) significant proxy penetration into new local constraints (this may happen when the user moves the virtual tool very quickly into a constraint); and (ii) existing local constraints with significant discontinuities in the penetration depth (this may happen when the virtual environment is generated using a penalty-based simulation, because the proxy interactions are approximations of the virtual tool interactions and the proxy penetration into the local constraints may differ from the virtual tool penetration into the virtual environment). Discontinuities in the local constraint penetration lead to discontinuities in the setpoint of the force control loop and may produce unacceptable force spikes or may destabilize the interaction. This difficulty is addressed by maintaining the proxy penetration into the local constraints continuous through proxy , i.e., the new local contact point is computed by (see Figure 6.7):  (6.4)  92  6.2 Local geometry  When the virtual environment sends an existing constraint to the local model, the local contact point is recomputed by: (6.5) The proxy is expanded back towards the actual geometry of the virtual tool whenever a local constraint becomes inactive. This approach is independent of the stiffness of the local constraints.  Figure 6.7: Local proxy deformation due to violation of new constraints.  Local proxy deformation maintains penetration continuity and eliminates force discontinuities at model updates. However, it allows the proxy to drift relative to the virtual tool. This drift is noticeable during fast motions in locally cluttered virtual environments, when significant proxy deformation may occur in order to maintain penetration continuity. Therefore, the free space perceived haptically may be much larger than the free space perceived visually, e-active local geometry is used in the proposed local model to alleviate drift.  6.2.3  e-active local geometry  To diminish drift and allow realistic haptic rendering of small clearances (in virtual environments generated using constraint-based simulations) and of tight constraints (in virtual worlds generated using penalty-based simulations), the local model is augmented by including constraints within some neighborhood of the virtual tool. The relevant neighborhood is defined by sweeping a sphere of radius e over the volume of the virtual tool, as described in [88] and depicted in Figure 6.8(a). If the virtual tool has concave surfaces, then typical simulation packages represent it as a collection of convex polyhedra [59, 85] and the augmented virtual tool (hereafter called e-active virtual tool) is obtained by sweeping a sphere of radius e over the volume of each component polyhedron, as illustrated in Figure 6.8(b). This technique reduces the drift by adding prediction capabilities to the proposed model: constraints are sent to the local model before they become active. Moreover, the approach is simple to implement and general enough to be applicable to any virtual environment. In addition to diminishing drift, the e-active virtual tool selects a unique constraint normal at degenerate contacts, as illustrated in Figure 6.9 for a vertex-vertex contact  6.2 Local geometry  93  Figure 6.8: Neighboring constraints included in the local model through the use of an e-active virtual tool, obtained by sweeping a sphere of radius e over the volume of the virtual tool.  of the virtual tool with the virtual environments. Hence, it alleviates the singularity in the constraint normal computation at degenerate contacts. Virtual tool  Figure 6.9: The e-active virtual tool eliminates the singularity in the constraint normal computation at a vertex-vertex contact between the virtual tool and the virtual environment by selecting a unique constraint normal direction.  At simulation updates, e-active geometry is included in the local model by computing the local contact position for new local constraints according to: P = R ( C \  +( - ) -c s  P  £  n  V COMJ  (6.6)  6.2 Local geometry  94  and by updating the local contact position for existing constraints according to:  p = R  c  if s„ < s — e <  0 or 0 < s - e  if s — e < s  0  V~P  P = R P = R  + sn  — c  <  p  if  s  —  e <  (6.7)  0<s  v  K~p  In Equations (6.6) and (6.7), s is the separation between the e-active virtual tool and the virtual environment as reported by the simulation, and s is the separation between the proxy and the local constraint in the local p  model of the interaction. The ability of the e-active geometry to improve user's perception of tightly constrained virtual tools is illustrated through a M a t l a b ™ simulation of a planar (3 D O F ) peg-in-hole manipulation. In this simulation, the commercial package generating the virtual environment is implemented as a penalty-based simulation of a l k g peg in a tight hole that fits the peg exactly (see Figure 6.10(a)). The computational delay of the package is chosen 20ms. Furthermore, the user's action is implemented as a horizontal force fh — 5sin ( f § ) N . In x  the simulation, the effect of the e-active geometry is illustrated separately from the effect of the impulsive forces and the effect of the haptic controller. The effect of the impulsive forces is eliminated by implementing penalty-based local constraints (with stiffness  K  t  contac  =  lOOOON/m and damping  B  c o n t  act  = 15N/(m/s)). As  depicted in Figure 6.10(b), the effect of the haptic controller is eliminated by directly coupling the local model to the haptic device (i.e., by directly applying the environment forces to the user's hand). The trajectory of the user's hand is plotted in Figure 6.11(a) for the case where only active geometry is included in the local model, and in Figure 6.11(b) for the case where e-active geometry within a 5mm neighborhood of the virtual tool is sent to the local model. As shown in the plots, the drift is significant without e-active geometry and is eliminated when an e-active virtual tool is employed. Note that the motion of the user's hand is due only to the limited stiffness of the local constraints in the latter case. 1 2  ms  Local model  Virtual environment (a) V i r t u a l t a s k .  (b) D y n a m i c m o d e l o f t h e t a s k .  Figure 6.10: One dimensional peg-in-hole task used to illustrate the influence of local geometry on user's perception of locally cluttered virtual environments.  6.2 Local geometry  95  Hand position  0  Hand position  5 time [sec]  1  0  5 time [sec]  (a) A c t i v e g e o m e t r y .  (b)  e-active g e o m e t r y ( w i t h i n 5 m m  10 of  the v i r t u a l tool).  Figure 6.11: Hand position for the simulated peg-in-hole task depicted in Figure 6.10(a).  6.2.4  Contact clustering and constraint coherence  Local proxy deformation and an e-active virtual tool smooth forces at model updates and improve the perception of tight virtual spaces provided that constraint continuity is ensured at these updates. Constraint continuity requires individual constraints to be identified and their local geometry to be updated according to Equations (6.5) and (6.6). In the local model of interaction, constraint continuity is maintained through contact clustering [85] and through maintaining constraint coherence based on geometric coherence [38]. Contact clustering is necessary in order to improve the stability of the haptic interaction in nonconvex regions of penalty-based virtual environments [85]. Typical collision detection algorithms employed in interactive simulations decompose nonconvex polyhedra into collections of convex pieces and approximate the penetration depths between pairs of overlapping convex pieces using algorithms such as proposed in [33,59,83,147]. In nonconvex regions of the virtual world, the convex decomposition of the virtual objects may result in a large number of convex pieces around concavities. In turn, this may lead to a large number of contacts and, hence, a larger resultant environment stiffness which may induce instability in certain configurations [85]. Contact clustering aims to avoid instability caused by large environment stiffnesses arising from the large number of contacts reported by the collision detection algorithms in concave regions of the virtual world. Contact clustering is performed in the virtual environment simulation using a translational proximity threshold. Specifically, a cluster is a collection of n  cc  contact point c  contacts of the virtual tool that satisfy the condition that their  is within the translational threshold c  tvE  from the contact points c  of all other contacts  in the cluster: l£  \<it  VE  P>*  Vj = l...n . cc  (6.8)  P,3  A cluster of virtual tool contacts in the virtual environment is represented by one local contact (see Figure 6.12).  6.2 Local geometry  96  The separation distance of the local contact is computed by:  s =  min  Sj  j  Vj = 1... n ,  (6.9)  cc  thus ensuring that the local contact exists as long as there are contacts in the cluster. Furthermore, the normal of the local contact is obtained averaging the normals of all cluster contacts:  E  n  ce  3=1  n  (6.10)  i  iE;=in ;  Note that the translational threshold needs to be chosen based on the stiffness of the contacts in the simulation. The more compliant the virtual environment, the larger this threshold needs to be.  Virtual environment Local constraint  /  s  s  °  \  Proxy  Virtual tool  (a) Cluster of contacts in the virtual en-  (b)  vironment.  cluster in Figure 6.12(a).  Local  constraint  representing the  Figure 6.12: Contact clustering in the virtual environment.  Besides contact clustering in the virtual environment, geometric coherence [38] is used in the local model to maintain the temporal coherence of the local constraints at model updates . In other words, constraints 2  existing at two consecutive updates, t — 1 and t, are considered the same local constraint if their contact points and normals are within some translational e  thM  and rotational t  ThM  proximity thresholds, respectively, from  each other:  |c  (6.11)  — c  (6.12) Time coherence breaks down if the virtual objects move more than prescribed by e  tLM  and t  TLM  during one  T h i s is necessary because the geometry of the local constraints does not include connectivity information. L a c k of connectivity ensures that less information is sent to the local model, that local collision detection is computationally inexpensive, and that the local model is compatible with any commercial simulation package.  6.2 Local geometry  97  simulation step. Therefore, these thresholds must be chosen based on the maximum speed of the objects in the virtual environment and the average time step of the simulation, At . avg  The translational proximity threshold  is computed by: e*LM  —  (6.13)  At gV , aV  where the maximum speed of the virtual objects v  max  max  is set by the user in the virtual environment and  At  avg  is computed in the local model by averaging the last five time steps of the simulation. Similarly, the rotational proximity threshold e  rLM  is computed using the maximum angular speed of the virtual objects u  max  set by  the user in the virtual world. In the local model, the rotational proximity threshold is also used to distinguish between discrete approximations of round virtual objects through polyhedra and polyhedral virtual objects (see Figure 6.13). Therefore, the rotational proximity threshold is lower bounded by the cosine of the maximum angle between two adjacent faces of a polyhedral object, cos(d ), as set by the user in the virtual environment max  simulation: r  £  LM  =max(At w ,cos(9 )). avg  max  (a) P o l y g o n a l a p p r o x i m a t i o n o f a r o u n d  rnax  (6.14)  (b) P o l y g o n a l v i r t u a l o b j e c t ,  virtual object.  Figure 6.13: Typical commercial simulation packages use polygonal representations (in a planar virtual world) for both rounded objects and polygonal objects.  The local geometry described in this section approximates the virtual environment geometry in close proximity to the virtual tool. The approximation is derived such that it eliminates the collision detection bottleneck from the local model of interaction and improves the user's perception of tight clearances. Furthermore, it is able to interface the haptic device to a virtual environment regardless of the algorithms and data structures employed to generate the rigid multibody virtual world.  6.3 Local dynamics  6.3  Local  98  dynamics  Besides local geometry, local dynamics affect the realism of the haptic feedback provided to users. They allow the local model to account for the inertia of the virtual tool during the computation of the virtual interactions. Furthermore, the local dynamics allow the local model to compute physically-motivated impulsive forces that are guaranteed not to increase the kinetic energy of the proxy during collisions. In essence, the local dynamics provide an efficient implementation of the simulation and the control techniques developed in Chapters 4 and 5. Their presentation in this section emphasizes the approximations that they involve. Local dynamics comprise the dynamics of a proxy of the virtual tool constrained by the e-active local geometry presented in Section 6.2. As described in the following section, the e-active local geometry is enforced through the impulse-augmented penalty approach introduced in Chapter 4. In the local model, this approach is based on a model of proxy contact derived from the rigid body contact model discussed in Section 4.2. The proxy contact model has three states: free motion, colliding contact, and resting contact. The proxy is said to be in free motion if it has no contacts (i.e., the separation distances between the proxy and all local constraints are strictly positive). The proxy is said to be in colliding contact if it has at least one new contact with negative normal contact velocity. Finally, the proxy is said to be in resting contact if it is neither in free motion nor in colliding contact. Using this proxy contact model, the constraints imposed on the virtual tool by the virtual environment are applied to users as detailed in the following section.  6.3.1  Local contact interactions applied t o users  Haptic manipulation of linkages is incorporated by defining the virtual tool to be the entire articulated structure when the user holds one of its links (see Figure 6.2) and by computing the proxy dynamics in configuration 3  space.  If the virtual tool is a rigid object, then its Cartesian space dynamics are its configuration space  dynamics. The computational load of the local model of interaction is much reduced compared to that of the virtual environment simulation. This is because local geometry makes local collision detection trivial and because only the virtual tool dynamics are simulated locally. During manipulation of linkages, numerical efficiency is further increased by approximating the configuration space linkage inertia D t, the configuration space gravitational v  terms G t, and the topological constraints with their values at the moment of the update. Furthermore, the v  Coriolis and centripetal effects are neglected, similar to work in [126,128]. Thus, the proxy dynamics during resting contact are simulated by: c  D tq+G„ = £ j f F - r - j £ F 1 )  t  i  h l  (6-15)  i=l  while its dynamics during colliding contact become: D „ q = D „ q + J p. t  3  t  0  c  B y t h i s d e f i n i t i o n , n e a r b y c o n t a c t s o f a l l l i n k s m u s t b e sent t o t h e l o c a l m o d e l .  (6.16)  6.3 Local dynamics  99  In Equations (6.15) and (6.16), D  v t  and G  VT  are the values of D  V T  and G„t, respectively, computed by the  simulation at the moment of the update. Furthermore, c is the number of proxy contacts, F j are penalty and friction forces, p is the collision impulse (computed as discussed in Section 4.4), Jj is the proxy Jacobian at the i-th proxy contact, Jc —  is the proxy Jacobian computed at the user's hand, F is the user-applied wrench, k  (JT i • • • Jf t • • • Jj c) > and rij is the local constraint normal at the i-th proxy contact. Passivity of n  n  n  T  proxy dynamics during collisions is guaranteed by assuming static constraints when resolving local collisions. This assumption prevents the simulation from inputting energy into the local model. Finally, the proxy contact interactions are applied to users using the approximate SVD of the hand Jacobian Jh, computed by the virtual environment simulation at the moment of the update: J =USiE V. fc  In Equation (6.17), S i € lZ  6xr  elsewhere, £  € W  xd  2  elsewhere, and r =  (6.17)  2  is a matrix that has the singular values of J on its main diagonal and zeros k  is a matrix whose elements are equal to one on the main diagonal and equal to zero  rank(Jh)-  Hence, the local coordinate invariant inverse of the linkage Jacobian at the user's  hand is:  J* = D ^ V E i " ( s V D 7 T  2  1 t  VE^)"  1  (sfU M UEi)~ T  w  1  £7v M , T  (6.18)  hl  where M u is the mass matrix of the link held by the user. Users feel an environment force:  upon proxy contact, and they feel:  Fenv = j f X X F j  ( - °) 6  2  i=l  during proxy contact.  6.3.2  Local control  As discussed in Chapter 5, linkage inertia and the motion constraints due to linkage topology are imposed on users through device control. When users operate virtual linkages through the local model of interaction, the configuration space inertia of the linkage D „ and the singular directions of j £ are approximated through their t  values in the virtual environment at the moment of the update. It follows that the inverse of the operational space inertia of the virtual tool at the user's hand A^" and the directions of topological constraint 1  n ,i tc  are also  approximated in the local model of interaction. Specifically, the local inverse of the operational space inertia of the virtual tool at the user's hand A^" is given by: 1  A^  1  -  JftD^Jh,  (6.21)  6.4 Experiments  100  and the control signal imposing the approximate topology on users is:  Uadd =  M  { tc^[ ,i k  d  c  ( /»„« ~ ft) + t n[c,i ( ^< - >0) Stc.i, x  x  X  b  (6-22)  x  c  i=l  where n  tc  and Htc,i are the number of topological constraints and the i-th topological constraint direction,  respectively, computed by the virtual environment simulation at the moment of the update. recall from Chapter 5 that Md is the mass matrix of the haptic device, kt and b c  tc  Furthermore,  are the stiffness and  damping of the topological constraints, x^ and x^ are the body position and body velocity of the user's hand, and x / i  v(  and x/j  vt  are the simulated body position and simulated body velocity of the virtual tool at the user's  hand in the local model of interaction. Equations (6.21)-(6.15) provide significant computational savings in the local model of the interaction and allow users to manipulate linkages with larger number of links. Realistic operation of virtual linkages using these approximations is demonstrated experimentally in the following section.  6.4  Experiments  In this section, the performance of the proposed local model of interaction is compared to the performance of the intermediate representation proposed in [21,22]. The intermediate representation comprises the geometry of the contacts of the virtual tool at the moment of the update. Furthermore, stable interaction within a virtual environment connected to the haptic interface through the intermediate representation is ensured through smoothing forces at updates using constraint interpolation, as proposed in [100]. Five experiments are performed in this section, four of which are controlled experiments. All experiments are implemented on the planar haptic simulation system described in Chapter 3 using asynchronous communication. The first experiment validates the passivity of the proxy dynamics during collisions in the impulse-augmented penalty local model. The second experiment demonstrates that the e-active geometry decreases the drift caused by local proxy deformation between the user-perceived constraints and the constraints in the virtual environment. The third experiment shows that the e-active geometry improves the haptic perception in locally cluttered virtual environments. The fourth experiment establishes that higher contact stiffnesses are achievable in the local model when a four channel haptic controller [136] renders the virtual interactions to the user then when the impedance haptic device is directly coupled to the admittance local model, as done in [21]. The last experiment suggests that further improvements are required to increase the realism of the contact with dynamic objects.  6.4.1  Passive collision dynamics i n the impulse-augmented penalty local model  The passivity of the proxy dynamics during collisions in the impulse-augmented penalty local model is investigated through user manipulations within the virtual environment shown in Figure 6.14(a) in its initial configuration. In the experiments described in this section, the user manipulates the distal link of the virtual linkage,  6.4 Experiments  101  as schematically depicted in Figure 6.14(b). The contact stiffness and damping are K  t  contac  and B  c o n t a c t  = 15000N/m  — 100N/(m/s), respectively. The contacts are perfectly plastic (e = 0) during a first example  interaction, and perfectly elastic (e = 1) during a second example interaction. The number of simultaneous collisions, the kinetic energy of the proxy, and the environment wrenches are monitored during interaction. The number of simultaneous collisions and the kinetic energy are plotted for both coefficients of restitution in Figure 6.15. Note that the kinetic energy of the proxy decreases during perfectly plastic collisions with the virtual environment (Figure 6.15(a)) and remains constant during perfectly elastic collisions with the environment (Figure 6.15(b)).  These experimental results validate the passivity of the proposed collision  resolution approach in the local model of interaction regardless of the contact restitution properties and of the number of simultaneous collisions. The passive collision resolution enhances the passivity of a simulation that computes impulsive and penalty interactions compared to a simulation that computes only penalty interactions. Therefore, the coupled stability of the haptic manipulation of rigid virtual tools increases within an impulseaugmented virtual environment compared compared to a penalties-only virtual world.  (a) Testbed virtual environment.  (b) Schematic of the interaction. Initial position is shown in black.  Sam-  ple later positions are shown in shades of grey (the lightest shade indicates the latest configuration).  Figure 6.14: Experimental interaction used to investigate the passivity of the collision resolution approach proposed in this thesis.  The environment wrenches during the two experimental interactions are shown in Figure 6.16 for perfectly plastic collisions, and in Figure 6.17 for perfectly elastic collisions. Because the wrenches applied to users during collisions are much larger than those applied to users during contact (by up to four orders of magnitude), these figures display the environment wrenches using a logarithmic scale. Recall that the perceived rigidity of a virtual environment can be enhanced by allowing users to feel large forces upon contact [90,125]. Hence, in addition to validating the passivity of the collision resolution, the experiments presented in this section confirm that the impulsive forces improve users' perception of rigidity within virtual worlds interfaced to the haptic  6.4 Experiments  102  Number of simultaneous collisions  Number of simultaneous collisions  20  0  30  10  20  Change of KE  30  40  50  60  Change of KE  -10  5  (a) P e r f e c t l y p l a s t i c c o n t a c t s , e = 0 .  10  15  20 25 time [s]  30  35  40  45  (b) P e r f e c t l y e l a s t i c c o n t a c t s , e — 1.  Figure 6.15: Passivity of proxy dynamics in an impulse-augmented penalty local model of interaction. The user operates the virtual linkage in Figure 3.3(b) from its distal link.  device through the proposed local model of rigid body interaction.  6.4.2  Transparency i n the local model of interaction  The first experiment described in this section investigates the effect of local e-active geometry on the transparency of haptic rendering of static rigid contacts. In this experiment, the user's action is represented by a constant wrench F ^ = (ON 0.5N 0 N m )  T 4  pushing the rectangular virtual tool shown in Figure 6.18(a) towards  the bottom-most horizontal wall of the rigid enclosure, as schematically depicted in Figure 6.18(b). The stiffness and the damping of the local contacts are K  c o n t a c  t —  4000N/m and B  contact  = 30N/(m/s), respectively.  These values represent the maximum contact impedance for which the interaction is stable in the intermediate representation (i.e., the virtual tool can be inserted in the tight-fitting hole at the bottom of the virtual world and remains stable upon being left there). Furthermore, collisions are considered perfectly plastic, i.e., e = 0. The device trajectories are plotted in Figure 6.19. These trajectories are obtained by interfacing the haptic device to the virtual world through: (i) the intermediate representation ("IR"); (ii) the proposed local model 5  including active geometry ("LMo"); and (iii) the proposed local model including e-active geometry within 5mm from the virtual tool ( " L M " ) . 5  The trajectories depicted in Figure 6.19 show the user's motion toward the virtual wall and the bounce once the virtual tool contacts the wall. Note that users perceive the virtual wall at different locations when the device 4 T h e f o r c e is l i m i t e d b y t h e p e r f o r m a n c e o f t h e v i r t u a l e n v i r o n m e n t ( i . e . , t h e v i r t u a l t o o l m o v e s f a s t e n o u g h t o p o p t h r o u g h the w a l l for larger h a n d forces). 5  S t a b l e i n t e r a c t i o n h a s b e e n e n s u r e d t h r o u g h s m o o t h i n g f o r c e s at u p d a t e s u s i n g c o n s t r a i n t i n t e r p o l a t i o n , as p r o p o s e d i n [100].  6.4 Experiments  „  103  F o r c e on hand along the x - a x i s  8  6  61  4  4  I  2 :  D) O  o  F o r c e on hand along the y - a x i s  2  o  -2 -2 -4 -4 -6 10  20  timefs]  30  40  -6 10  (a)  20  time[s]  30  (b)  Torque on hand  Figure 6.16: Environment wrenches applied to users during user manipulation of the distal link of the articulated object in Figure 6.14(a) in an impulse-augmented penalty local model of interaction with perfectly plastic contacts, e = 0.  40  6.4 Experiments  104  F o r c e on hand along the x - a x i s  F o r c e on hand along the y - a x i s  Torque on hand  E 2-2 O)  -5-4  -6 10  20  time[s]  30  40  (c)  Figure 6.17: Environment wrenches applied to users during user manipulation of the distal link of the articulated object in Figure 6.14(a) in an impulse-augmented penalty local model of interaction with perfectly elastic contacts, e = 1.  6.4 Experiments  105  (a) Testbed virtual environment.  (b) Schematic of the interaction. Initial position is shown in black.  Sam-  ple later positions are shown in shades of grey (the lightest shade indicates the latest configuration).  Figure 6.18: Experiment used to investigate the haptic rendering of static contact.  (a)  (b)  Figure 6.19: Device trajectories obtained when the virtual tool shown in Figure 6.18(a) is pushed with a constant force towards the bottommost horizontal wall of the virtual environment, as depicted in Figure 6.18(b). The intermediate representation ("IR"), the local model including active geometry ( " L M " ) , and the local model including e-active geometry within 5mm from the virtual tool ( " L M 5 " ) interface the haptic device to the virtual environment simulation. n  6.4 Experiments  106  is interfaced to the virtual world through the intermediate representation, the local model with active geometry, and the local model with e-active geometry. In particular, users perceive the virtual world at its location in the virtual environment, ycOM  = 116mm, when e-active geometry within 5mm of the virtual tool is included in  the local model. Users perceive the virtual wall at ycoM = 87mm and at ycOM  = 80mm in the intermediate  representation and in the local model with active geometry, respectively . The drift between the perceived 6  location of the wall and its location in the virtual environment is caused by the smoothing techniques used to avoid force discontinuities at updates in the intermediate representation and in the local model, constraint interpolation and proxy deformation, respectively. Since users perceive no drift when interacting with the local model with e-active geometry, the experiment illustrates that no proxy deformation occurs in this case. In turn, this means that constraints are sent to the local model before they become active. Hence, the experiment demonstrates that e-active geometry diminishes the drift between the user-perceived position of the constraint and its position in the virtual environment. While the drift is dependent on the user-applied forces and the virtual world geometry, the experiment illustrates that it can be significantly reduced (by approximately 3.5cm in this case) by including in the local model e-active constraints within a relatively small neighborhood of the virtual tool (5mm in this experiment). By alleviating this drift, the e-active geometry increases the transparency of the interaction. The second experiment demonstrates the role of the e-active geometry in improving the perception of tight virtual clearances. In this experiment, the controlled interaction represents a peg-in-hole manipulation. The user rotates the rectangular virtual tool shown in Figure 6.20(a) by 90° and inserts it into the hole at the bottom of the rigid enclosure that exactly fits the peg. The user then releases the peg and their hand is replaced by the wrench F^, = (lsin(57r)N ON 0 N / m ) , i.e., the peg is shaken horizontally by a sinusoidally varying force, T  as schematically represented in Figure 6.20(b). As before, the stiffness and damping of the local contacts are K t ct con a  = 4000N/m and B ct conta  — 30N/(m/s), respectively. The experimental device trajectories are  shown in Figure 6.21. They are obtained by interfacing the haptic device to the virtual environment through: (i) the intermediate representation ("IR"); (ii) the present local model including active geometry ("LMrj"); and (iii) the local model including e-active geometry within 5mm of the virtual tool ( " L M 5 " ) . The experimental trajectories illustrate that, depending on the technique employed to haptically render the virtual hole, the device travels different distances along the x direction at different locations along the y direction in the virtual world. The differences in the three trajectories result from two combined factors. First, the constraints (representing both the lateral and the top virtual walls) arrive with delay in the intermediate representation and in the local model with active geometry. Hence, not all constraints representing the hole geometry are avaialable in these models throughout the interaction and these models render only approximations of the hole geometry to the user. The constraints arrive before becoming active in the local model with e-active geometry. Therefore, all constraints representing the hole geometry are available in this model throughout the interaction and this model renders the exact hole geometry to the user. Approximate 6  T h e different locations are due to the slightly different s t a r t i n g positions of the device. T h e s e differences are caused by f r i c t i o n  i n t h e d e v i c e a n d t h e l i m i t e d s t i f f n e s s of t h e c o n t r o l l e r u s e d t o d r i v e t h e d e v i c e t o t h e i n i t i a l p o s i t i o n .  6.4 Experiments  107  Figure 6.20: Experiment used to demonstrate the role of e-active geometry in improving user's perception of tight virtual clearances.  Device trajectory '  (a)  1  •  1  1  , „  Device rotation  (b)  Figure 6.21: Device trajectories obtained when the virtual tool shown in Figure 6.20(a) is rotated by 90°, inserted into the tight-fitting hole at the bottom of the virtual world, and shaken horizontally with a sinusoidally varying force, as depicted in Figure 6.20(b). The intermediate representation ("IR"), the local model including active geometry ("LMo"), and the local model including e-active geometry within 5mm from the virtual tool ( " L M 5 " ) interface the haptic device to the virtual environment simulation.  6.4 Experiments  108  (i.e., incomplete) hole geometry accounts for the larger device rotation and for the larger free space haptically perceived by users (represented through larger device motion along the x direction) in the IR and the LMn trajectories in Figure 6.21. Second, the constraints are imposed on the device using only penalty forces in the intermediate representation and using impulsive and penalty forces in the local model of interaction. Approximate (incomplete) hole geometry together with the contact model explain the variances along the y axis in the device trajectories. In particular, when the virtual constraints are imposed through penalty forces (the IR trajectory), the device bounces and loses contact with the constraints. During free motion, the interface drift along the y direction is unbalanced and the device moves away from the horizontal wall of the hole, as illustrated in Figure 6.21(a). When the virtual constraints are imposed through impulsive and penalty forces (the L M  0  trajectory), the device maintains contact with the lateral constraints longer and friction in the lat-  eral walls balances the interface drift. The interface motion towards the bottom of the hole is caused by the intermittent loss of contact with the bottom wall, which results in the horizontal constraints arriving late to the local model with active geometry. In turn, the delayed horizontal constraints cause local proxy deformation along the y-direction and allow the interface to drift towards the bottom of the hole.  6.4.3  The effect of the teleoperation controller on the achievable contact stiffness  To demonstrate the effect of the teleoperation controller on the achievable contact stiffness, the controlled pegin-hole experiment depicted in Figure 6.20 is repeated using only the local model of the interaction and increasing the stiffness and the damping of the local contacts to  K  c o n t a  ct  =  15000N/m and B  c o n t a  ct =  100N/(m/s),  respectively . This stiffness is almost four times larger than the contact stiffness used in the previous exper7  iments. Moreover, the controlled wrench applied on the device is F = (3sm(57r)N ON 0 N / m ) . The device T  trajectories obtained using both active ( " L M " ) and predicted ( " L M " ) geometry in the local model of the 0  5  interaction are plotted in Figure 6.22. The experimental trajectories show that much higher contact stiffness can be rendered to users when a four channel controller applies the virtual interactions to the user's hand then when the impedance haptic device is directly coupled to the admittance local model of interaction. Furthermore, these trajectories validate the positive effect of predicted geometry on the perception of tight clearances. Note that the device travels approximately 2cm more along both the x and the y axes and rotates more around the horizontal orientation of the virtual peg when the local model incorporates active geometry than when it incorporates predicted geometry within 5mm of the virtual tool. However, due to overshoot, the device penetrates into the constraints more than 5mm and the neighborhood selected for geometry prediction is not sufficient to eliminate the impact of the virtual environment delay on the interaction represented by the L M 5 trajectory. This impact is represented in the larger device motion along the x axis than permitted by the contact stiffness. 7  O n l y t h e l o c a l m o d e l o f i n t e r a c t i o n is u s e d t o i n t e r f a c e t h e d e v i c e t o t h e s i m u l a t i o n b e c a u s e t h e l a r g e r c o n t a c t s t i f f n e s s r e n d e r s  t h e i n t e r a c t i o n u n s t a b l e i f t h e d e v i c e i s d i r e c t l y c o u p l e d t o t h e i n t e r m e d i a t e r e p r e s e n t a t i o n p r o p o s e d i n [21]  6.4 Experiments  109  Device trajectory  (a)  Device rotation  (b)  Figure 6.22: Device trajectories obtained when the virtual tool is inserted in the tight-fitting hole at the bottom of the virtual world shown in Figure 6.20(a) and shaken horizontally with a sinusoidally varying force, as depicted in Figure 6.20(b). The local model including active geometry ("LMn") and the local model including predicted geometry within 5mm from the virtual tool ( " L M 5 " ) interface the haptic device to the virtual environment simulation.  6.4.4  Dynamic contact  The controlled manipulation described in this section investigates the realism of the haptic interaction with moving virtual objects using the virtual environment depicted in Figure 6.23(a). To allow comparison between the local model and the intermediate representation, the virtual contacts have stiffness 4000N/m and damping  B tact con  =  K  c o n  t ct a  =  30N/(m/s). Moreover, the user's hand is replaced by a constant wrench  F/j = (ON 0.1N 0 N / m ) that pushes the virtual tool until the last link of the virtual linkage hits the top-most T  horizontal wall of the enclosure (placed at ycoM = 181mm), as schematically depicted in Figure 6.23(b). Figure 6.24 depicts the device trajectories obtained by interfacing the haptic device to the virtual world through: (i) the intermediate representation ("IR"); (ii) the local model including predicted geometry within 5mm of the virtual tool and static constraints ("LMsc")i and (iii) the local model including predicted geometry within 5mm of the virtual tool and moving (kinematic) constraints ("LM^c?"). In the local model, the kinematic constraints move along the contact normal with the velocity that they have in the virtual environment at the moment of the update. Several remarks result from the experimental trajectories. First, the haptic device (i.e., the user's hand) loses contact with the virtual link several times during pushing, regardless of the technique used to interface the device to the simulation. This non-physical behavior is demonstrated by the bounces recorded in all trajectories in Figure 6.24. Second, the device bounces more when connected to the virtual world through the intermediate representation then when connected through the local model of interaction. Third, the bounces are similar when the local model of interaction interfaces the device to the virtual environment, regardless of whether  6.4 Experiments  l  110  Vjticulatec i strix<;ttoe  Rigid ^nclogura (tllllKWJKWLi.  Virtual tool Moving  irtual  Object |  tdol  ^77777777777  (a) Testbed virtual environment.  (b) Schematic of the interaction. Initial position is shown in black.  Sam-  ple later positions are shown in shades of grey (the lightest shade indicates the latest configuration).  Figure 6.23: Experiment used to investigate the realism of the haptic rendering of dynamic contact.  Device rotation  Device trajectory  175 170 165  |  160  §155  o  ^150 145  i  j if  fj  140  24  26  CO M  X  (a)  8  r i mm  30  2 time [s]  3  0>)  Figure 6.24: Device trajectories obtained when the virtual tool pushes the last link of the virtual linkage in Figure 6.23(a) with a constant force, as shown in Figure 6.23(b). The intermediate representation ("IR"), the local model including predicted geometry and static constraint ( " L M s c " ) , and the local model including prdicted geometry and moving (kinematic) constraints ( " L M x c " ) interface the haptic device to the virtual environment simulation.  6.4 Experiments  111  the local constraints are static or kinematic. Two conclusions can be drawn from these experimental results. First, further improvements are required both in the intermediate representation and in the local model of interaction to increase the realism of dynamic contact, i.e., to maintain contact while interacting with moving virtual objects. Second, the local contact model enhances the stability of dynamic contact compared to the intermediate representation. This improvement is due to the local modeling of collisions and is not influenced by the use of kinematic local constraints.  This chapter has proposed the first local model of rigid body interaction that has been used to constrain both the translation and the rotation of a haptic device. This model comprises geometry and dynamics techniques that enable efficient and transparent haptic rendering of rigid body motion with constraints. Local geometry maintains force continuity at model updates through implementing novel techniques like local proxy deformation and e-active geometry. Local dynamics efficiently implements the simulation and control methods developed in Chapters 4 and 5. The model has been validated by comparing its performance to the performance of existing approaches through experiments performed using a planar haptic interaction system. The experiments confirm that the local model proposed in this thesis increases the stability and the transparency of rigid body interaction within multibody virtual environments. This model allows users to feel collisions upon contact, stiffer constraints during contact, and topological constraints during linkage manipulation. The model also allows haptic manipulation of virtual linkages from any user-selected link, as demonstrated by the experiments presented in Chapter 5.  In addition, the proposed local model of interaction provides a general purpose  methodology for interfacing a haptic device to a virtual environment simulation, as it imposes no restrictions on the data structures and the algorithms used to generate the rigid multibody virtual world.  Chapter  7  Conclusion This thesis has proposed a novel approach to applying realistic forces to users while they operate rigid virtual tools within multibody virtual worlds. In this approach, an original model of unilateral contact significantly increases the stability and the perceived rigidity of virtual contacts. A new model of bilateral contact and a coordinate invariant mapping of interactions from configuration space to operational space enable users to freely manipulate virtual linkages in addition to virtual objects. Fast local approximations of rigid body interaction interface a haptic device to a virtual environment simulation in a manner transparent to the application developer.  As a result, convincing force feedback can easily be added to virtual prototyping and design  applications, in a modular fashion and without a detailed consideration of the haptic simulation and control algorithms. The specific contributions to haptic rendering of rigid body motion with constraints of the thesis are summarized in the following section.  7.1  Summary of Contributions  • Rigid contact modeling A novel model of unilateral rigid body contact has been introduced that enables users to feel collisions. In the proposed model, contacts are infinitely stiff when they arise and have limited stiffness thereafter. To achieve this effect, the model combines constraint-based collision resolution with penalty-based contact resolution. This allows interaction forces and impulses to be computed and applied to users at the high fixed rates imposed by the haptic controller. The approach eliminates the need for constraint stabilization following collision resolution. A new collision resolution method has been developed that allows multiple collisions to be resolved simultaneously such that the kinetic energy of the multibody system never increases during collisions. The method is shown to apply equally to independent and to overdetermined constraints, as it implicitly eliminates constraint overdeterminancy. A physically motivated approach has been proposed to incorporate the limitations of the haptic device 112  7.1 Summary of Contributions  113  when impulses are applied to the user. The approach is consistent with conservation of energy principles and is made possible by the fact that constraints need not be stabilized after impact. Furthermore, the approach allows graceful degradation of performance when device capabilities diminish by maximizing the user's hand kinetic energy dissipated during impacts subjected to device capabilities. A friction model developed in computational mechanics [61] has been used for the first time for haptic rendering of dry friction. Simulations and experiments have shown that this model renders a richer set of sticking contact characteristics than prior methods for haptic rendering of Coulomb friction. These characteristics are easily adjusted through a parameter with physical significance. • Haptic manipulation of linkages Realistic manipulation of virtual linkages from any user-selected link and through singularities has been enabled through the development of new simulation and control methods. To increase computational efficiency and to eliminate the need for constraint stabilization, the dynamics of the multibody system have been simulated in configuration space. To enable compelling and unrestricted operation of virtual linkages, the configuration space dynamics have been mapped to operational space and applied to users through control. A four channel teleoperation controller imposes the unilateral constraints on users, while an augmented impedance-controller applies the linkage inertia and the bilateral (topological) constraints to users. The coordinate invariant mapping of the linkage dynamics from configuration space to operational space and the representation of the topological constraints at the user's hand through control contribute new developments to haptic manipulation of virtual linkages. The coordinate invariant representation of linkage contacts at the user's hand is based on the theory of weighted generalized inverses developed in robotics [46] for purposes of hybrid control. While the mapping itself is not a new theoretical development, it has not been used in prior simulation research. The topological constraints are imposed on users through penalizing their departure from the configuration manifold of the virtual linkage through device control. The penalties embody the directions along which the virtual linkage has infinite structural stiffness at the user's hand. Simulations and experiments have been used to demonstrate that the topological penalties successfully impose the bilateral constraints on users. • Fast local approximations of rigid body interaction Local models of rigid body interaction have been considered in prior research [21]; however, this work introduces the first local model that has been used to constrain both the translation and the rotation of the haptic device . Five new features distinguish the proposed local model of rigid body interaction 1  from existing models. First, the proposed local model uses local proxy deformation to avoid destabilizing force discontinuities at model updates. Experiments performed in a planar testbed virtual environment 1 D e v i c e r o t a t i o n is t y p i c a l l y m o r e c h a l l e n g i n g t o c o n s t r a i n b e c a u s e , i n t y p i c a l s i m u l a t i o n p a c k a g e s , r i g i d b o d y c o n t a c t i n c o n c a v e a r e a s o f t h e v i r t u a l e n v i r o n m e n t is r e p r e s e n t e d t h r o u g h a l a r g e n u m b e r o f v e r t e x - f a c e c o n t a c t s w h o s e s t i f f n e s s e s c o m p o u n d s u c h t h a t t h e y r e s u l t i n l a r g e r i n c r e a s e s i n t h e r o t a t i o n a l s t i f f n e s s t h a n i n t h e t r a n s l a t i o n a l s t i f f n e s s o f t h e v i r t u a l e n v i r o n m e n t [85].  7.2 Future Work  114  demonstrate that local proxy deformation successfully maintains interaction stability in contact configurations that present challenges for prior local models (like peg-in-hole manipulations [21]). Second, the model comprises e-active geometry, which has been demonstrated to increase the realism of the haptic manipulation of the virtual tool through small clearances. Third, the proposed local model implements the impulse-augmented penalty model of rigid contact developed in this thesis. Hence, it computes impulsive forces upon contact and penalty and friction forces during contact. The impulsive forces increase the stability of the local contacts. Fourth, the proposed local model includes a dynamic proxy in addition to local geometry. The dynamic proxy allows a four channel haptic controller [136] to be used for transmitting wrenches between the device and the proxy and for coordinating positions between them. Compared to the virtual coupler [3,30], the four channel controller applies the locally computed impulsive forces to the user's hand, thus enabling users to feel collisions and improving the perceived rigidity [90,125] of the virtual contacts. Compared to directly coupling an impedance device to an admittance simulation [21], the four channel controller allows users to feel stiffer contacts. Lastly, the proposed local model of interaction imposes no restrictions on the data stuctures and the dynamic response algorithms used to generate the virtual environment. Therefore, the model can be used to add force feedback to interactive virtual worlds generated using any commercial simulation package.  7.2  Future Work  The results of this work can be extended and its limitations addressed in several ways: • Improved contact models Compared to prior models used for haptic rendering of rigid contact, the model developed in this thesis constrains the user's hand trajectory closer to the trajectory imposed by a real rigid wall. However, the performance of this model is little improved during contact with dynamic virtual environments. This is expected, since the model is derived from a penalty representation of static contact. More realistic models of dynamic contact can be developed that include inertial properties in addition to stiffness and damping. Such impedance-type representations of rigid contact may be developed that require minimal additional computational effort and more accurately embody dynamic contact. • Collisions with friction and arbitrary coefficients of restitution The present collision resolution approach assumes frictionless collisions, all with the same coefficient of restitution. This assumption has successfully increased the stability and the perceived rigidity of the virtual contacts. However, it cannot be used to render rolling, sticking, or reversed sliding that might occur after collision. The development of a collision resolution method compatible with the fixed haptic rate requirements would allow such phenomena to be rendered to users, thus increasing the realism of the haptic manipulation of rigid objects within multibody virtual worlds.  7.2 Future Work  115  • Local geometry prediction based on user's intent Local models of interaction are necessary for allowing haptic collaboration over internet. This work has been the first to allow rigid body manipulation using such a model. However, only crude geometry has been included in the proposed local model. As environment geometry is significant during the interaction with rigid virtual objects, improved prediction of local geometry is likely to result in enhanced realism of the interaction. The enhancement will allow the local model to suitably render more complex virtual environments, i.e., virtual environments with larger computational delays and more complex contact geometries. One promising approach to enriching the local geometry is contact prediction based on user intent rather than based on spatial proximity. Various mechanisms can be used to predict user intent, from simple ones (like user velocity) to others more complex (like Hidden Markov Models). • Virtual geometry for haptics A wider spread of haptic applications requires the availability of general purpose interfaces between haptic devices and rigid multibody simulations like the local model proposed in this thesis. Such interfaces allow application developers to add force feedback to virtual environments without requiring them to have detailed knowledge of the haptic simulation and control algorithms. However, the quality of the haptic feedback provided by such interfaces depends on the suitability of the data structures used in the virtual environment simulation. Present data structures result in perceptual artifacts during haptic simulations, since they are typically not able to maintain contact coherence during contacts that involve concave surfaces. In this work, such limitations have been alleviated through contact clustering and through maintaining temporal coherence based on spatial coherence. However, this technique breaks down if the user moves fast in highly concave regions. Geometric data structures that can be used to maintain rigid body contact coherence must be developed to eliminate such artifacts. • Alternative collision rendering techniques and user studies In this work, collisions have been rendered to users through impulsive forces that dissipate the kinetic energy predicted by the coefficient of restitution over one simulation step. This approach is based on the assumption that impacts are instantaneous events. Nonetheless, alternative rendering approaches might suitably approximate the instantaneous nature of collisions while allowing users to perceive other contact characteristics in addition to rigidity. Sinusoidally varying impulsive forces might allow users to distinguish between wood on wood or wood on metal collisions, for example. Future user studies need to be performed to identify salient features of the impulsive forces that can be manipulated in order to allow users to perceive such differences. Furthermore, the model used in this work to render dry friction has been shown adequate for displaying various stiction behaviors. Whether these behaviors can be associated to meaningful contact characteristics must be investigated through future user studies.  7.2 Future Work  116  • User studies The haptic rendering techniques developed in this work have been shown to enable more transparent manipulations of rigid virtual objects than prior methods.  The transparency of the interaction has  been assessed via comparing physical outcomes of virtual manipulations implemented through varying modeling approaches to physical outcomes of real manipulations. No attempt has been made to assess the performance of the proposed methods via perceptual studies. This is partly because the motivation for some of the new techniques (such as the use of impulsive forces for haptic rendering of rigid contact) has been provided by existing user studies [125], [90]. Nevertheless, human perception is key to evaluating any haptic rendering approach. Therefore, future user studies are necessary to further validate the advantages of the methods developed in this work. • Guaranteed stable interaction Lastly, the present work does not guarantee stable interaction within arbitrary rigid multibody virtual environments. This is because contacts are resolved using traditional penalty techniques. These techniques are conditionally stable and render environment stiffnesses that depend on contact geometry. Hence, their stability cannot be guaranteed for arbitrary virtual environments. This limitation needs to be addressed by developing a new concept of stiffness for rigid body constraints that breaks away from combining point stiffnesses into a resultant rigid body contact stiffness. The successful development of such a concept hinges on a physically meaningful metric for defining distances in configuration space. At the end of this work, rigid body interaction within arbitrary virtual environments continues to be a challenging endeavor. The developments presented in this thesis have led to increased stability and enhanced realism of the force feedback applied to users while they manipulate virtual objects or virtual linkages. However, they do not guarantee that the interaction will be free of perceptual artifacts or stable regardless of the complexity of the virtual world. Further work is required to provide such guarantees and to transform haptic devices into widely used computer interfaces.  Bibliography [1] M . Abadie. Dynamic Simulation of Rigid Bodies: Modeling of Frictional Contact. In B . Brogliato, editor, Impacts in Mechanical Systems: Analysis and Modeling, pages 62-144. Springer Verlag, Berlin, Lecture Notes in Physics, 2000. [2] Y . Adachi, T. Kumano, and K . Ogino. Intermediate Representation for Stiff Virtual Objects. In Proceedings of the IEEE Virtual Reality Annual International Symposium, pages 203-210, Research Triangle Park, N C , 1995. [3] R . J . Adams and B . Hannaford. Stable Haptic Interaction with Virtual Environments. IEEE Transactions on Robotics and Automation, 15(3):465-474, June 1999. [4] S. Ahmed, H . M . Lankarani, and M.F.O.S. Pereira. Frictional Impact Analysis in Open-Loop Multibody Mechanical Systems. ASME Journal of Mechanical Design, 121:119-127, March 1999. [5] M . Anitescu and G . D . Hart. A Hard-Constraint Time-Stepping Approch for Rigid Multibody Dynamics with Joints, Contact, and Friction. In Proceedings of the 2003 Conference on Diversity in Computing, pages 34-41, Atlanta, G A , 2003. [6] M . Anitescu and G.D. Hart. A Constraint-Stabilized Time-Stepping Approach for Rigid Multibody Dynamics with Joints, Contact and Friction. International Journal for Numerical Methods in Engineering, in press 2004. [7] M . Anitescu and F A . Potra. Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems. ASME Nonlinear Dynamics, 14:231-247, 1997. [8] M . Anitescu, D.E. Stewart, and F.A. Potra. Time-Stepping for Three-Dimensional Rigid Body Dynamics with Friction. Computer Methods in Applied Mechanics and Engineering, 177(3-4):183-197, 1999. [9] 0 . Astley and V Hayward. Multirate Haptic Simulation Achieved By Coupling Finite Element Meshes Through Norton Equivalents. In Proceedings of the 1998 IEEE International Conference on Robotics and Automation, pages 989-994, Leuven, Belgium, 1998. [10] D. Baraff. Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies. Computer Graphics (Proceedings of the 1989 SIGGRAPH), 23:223-231, 1989. [11] D. Baraff. Issues in Computing Contact Forces for Non-Penetrating Rigid Bodies. Algorithmica, 10(24):292-352, 1993. [12] D. Baraff. Fast Contact Force Computation for Non-penetrating Rigid Body Simulation. Computer Graphics (Proceedings of the 1994 SIGGRAPH), 29:23-34, 1994. [13] D. Baraff. Interactive Simulation of Solid Rigid Bodies. IEEE Computer Graphics and Applications, 15:63-75, 1995. [14] F. Barbagli, D. Prattichizzo, and J . K . Salisbury. Dynamic local models for stable multi-contact haptic interaction with deformable objects. In Proceedings of 11th Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems, pages 109-116, Los Angeles, C A , 2003. [15] A . A . Barhorts and L . J . Everett. Contact/Impact in Hybrid Parameter Multiple Body Mechanical Systems. ASME Journal of Dynamic Systems, Measurement, and Control, 117:559-569, 1995. 117  BIBLIOGRAPHY  R. Barzel and A . H . Barr. A Modeling System Based on Dynamic Constraints. (Proceedings of the 1994 SIGGRAPH), 22:179-188, 1988.  118  Computer Graphics  C. Basdogan, C H . Ho, and M . A . Srinivasan. A Ray-Based Haptic Rendering Technique for Displaying Shape and Texture of 3D Objects in Virtual Environments. In Proceedings of ASME Dynamics Systems and Control Division, volume DSC-61, pages 77-84, 1997. J.A. Battle and S. Cardona. The Jamb (Self-Locking) Process in Three-Dimensional Collisions. ASME Journal of Applied Mechanics, 65:417-423, June 1998. J. Baumgarte. Stabilization of Constraints and Integrals of Motion in Dynamical Syatems. Computer Methods in Applied Mechanics and Engineering, 1:1-16, 1972. A. Ben-Israel and T . N . E . Graville. Generalized Inverses: Theory and Applications. New York:Wiley, 1974. J . M Berkelman. Tool-Based Haptic Interaction with Dynamic Physical Simulations using Lorentz Magnetic Levitation. P h D thesis, Carnegie Mellon University, 1999. P.J. Berkelman, R . L . Hollis, and D. Baraff. Interaction with a Realtime Dynamic Environment Simulation using a Magnetic Levitation Haptic Interface Device. In Proceedings of the 1999 International Conference on Robotics and Automation, pages 3261-3266, Detroit, Michigan, 1999. V . Bhatt and J . Koechling. Partitioning the Parameter Space According to Different Behaviors During Three-Dimensional Impacts. ASME Journal of Applied Mechanics, 62(3):740-746, 1995. V . Bhatt and J . Koechling. Three-Dimensional Frictional Rigid-Body Impact. ASME Journal of Applied Mechanics, 62(4):893-898, 1995. C S . Bonaventura and K . W . Jablokow. A n O(N) Modular Algorithm for the Dynamic Simulation of Robots Constrained by a Single Contact. The IEEE Transactions on Systems, Man, and Cybernetics: Part C, 32(4):406-415, 2002. R . M . Brach. Rigid Body Collision. ASME Journal of Applied Mechanics, 56:133-138, 1989. B. Brogliato. Nonsmooth Mechanics, Second Edition. Springer-Verlag, London, 1999. J . M Brown. Passive Implementation of Multibody Simulation for Haptic Display. P h D thesis, Nothwestern University, 1995. J . M Brown. A Theoretical and Experimental Investigation into the Factors Affecting the Z-width of a Haptic Display. Master's thesis, Nothwestern University, 1998. J . M . Brown and J . E . Colgate. Physics-based Approach to Haptic Display. In Proceedings of the Topical Workshop on Virtual Reality, International Symposium on Measurement and Control in Robotics, volume 1, pages 101-106, Houston, T X , 1994. J . M . Brown and J.E. Colgate. Minimum Mass for Haptic Display Simulations. In Proceedings of the 1998 ASME International Mechanical Engineering Congress and Exhibition, pages 85-92, Nashville, T E , 1997. M . Buck and E . Schomer. Interactive Rigid Body Manipulation with Obstacle Contacts. Visualization and Computer Animation, 9(4):243-257, 1998.  Journal of  S. Cameron. Enhancing G J K : Computing Minimum and Penetration Distance between Convex Polyhedra. In Proceedings of the IEEE International Conference on Robotics and Automation, pages 3112-3117, 1997. B. Chang and J . E . Colgate. Real-Time Impulse-Based Simulation of Rigid Body Systems for Haptic Display. In Proceedings of the 1997 ASME International Mechanical Engineering Congress and Exhibition, pages 1-8, Dallas, Texas, 1997. A. Chatterjee and A . Ruina. A New Algebraic Rigid Body Collision Law based on Impulse Space Considerations. ASME Journal of Applied Mechanics, 65(4):939-951, 1998.  BIBLIOGRAPHY  119  J. Chen, C. DiMattia, R . M . Taylor II, M . Falvo, P. Thiansathaporn, and R. Superfine. Sticking to the Point: a Friction and Adhesion Model for Simulated Surfaces. Proceedings of the ASME Dynamic Systems and Control Division, DSC-61:167-171, 1997. M . Cline and D. Pai. Stabilization for Rigid Body Simulation with Contact and Constraints. In Proceedings of the IEEE International Conference on Robotics and Automation, pages 3744-3751, Taipei, Taiwan, September 2003. J.D. Cohen, M . C . Lin, D. Manocha, and M . K . Ponamgi. I-COLLIDE: A n Interactive and Exact Collision Detection System for Large-Scale Environments. In Symposium on Interactive 3D Graphics, pages 189196, Monterey, C A , August 1995. J.E. Colgate. Coupled Stability of Multiport Systems - Theory and Experiments. Transactions of the ASME, Journal of Dynamic Systems, Measurement, and Control, 116(3):419-428, 1994. J.E. Colgate and J . M . Brown. Factors Affecting the Z-width of a Haptic Display. In Proceedings of the 1994 IEEE International Conference on Robotics and Automation, pages 3205-3210, San Diego, C A , 1994. J.E. Colgate and G . G . Schenkel. Passivity of a Class of Sampled-Data Systems: Application to Haptic Interfaces. Journal of Robotic Systems, 14(l):37-47, 1997. J.E. Colgate, M . C . Stanley, and J . M . Brown. Issues in the Haptic Display of Tool Use. In Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 140-145, Pittsburgh, PA, 1995. D. Constantinescu, I. Chau, L . Filipozzi, S.P. DiMaio, and S.E. Salcudean. Haptic Rendering of Planar Rigid-Body Motion using a Redundant Parallel Mechanism. In Proc 2000 IEEE Int Con} Rob Autom, pages 2440-2445, Berkeley, Ca, 2000. D. Constantinescu, S.E. Salcudean, and E . A . Croft. Haptic Feedback using Local Models of Interaction. In 11th International Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems, pages 416-421, Los Angeles, Ca, 2003. P.R. Dahl. A Solid Friction Model. Technical Report TOR-158(3107-18), The Aerospace Corporation, E l Segundo, C A , May 1968. K . L . Doty, C. Melchiorri, and C. Bonivento. A Theory of Generalized Inverses applied to Robotics. The International Journal of Robotics Research, 12(1):1-19, February 1993. P.E. Dupont. The Effect of Coulomb Friction on the Exstence and Uniqueness of the Forward Dynamics Problems. In The IEEE Proceedings of the 1992 International Conference on Robotics and Automation, pages 1442-1447, Nice, France, 1992. R. Featherstone. Robot Dynamics Algorithms. BostomKluwer, 1987. R. Featherstone. A Divide-and-Conquer Articulated-Body Algorithm for Parallel o(logn) Calculation of Rigid-Body Dynamics: Part LBasic Algorithm. The International Journal of Robotics Research, 18(9):867-875, 1999. R. Featherstone. A Divide-and-Conquer Articulated-Body Algorithm for Parallel o(logn) Calculation of Rigid-Body Dynamics: Part 2: Trees, Loops, and Accuracy. The International Journal of Robotics Researc/i,18(9):876-892, 1999. J  J. Garcia de Jalon and E . Bayo. Kinematics and Dynamics Simulation of Multibody Systems: The Real-Time Challenge. Springer-Verlag, 1993. B. Gillespie. Haptic Display of Systems with Changing Kinematic Constraints: The Virtual Piano Action. PhD thesis, Stanford University, 1996. R . B . Gillespie. Kane's Equations for Haptic Display of Multibody Systems. Haptics-e, 3(2):l-20, 2003.  BIBLIOGRAPHY  120  [54] R . B . Gillespie and M . R Cutkosky. Interactive Dynamics with Haptic Display. Advances in Robotics, Mechatronics, and Haptic Interfaces, ASME, DSC-49:65-72, 1993. [55] B . G . Gilmore and D.A. Streit. A Rule Based Algorithm to Predict the Dynamic Behavior of Mechanical Part Orienting Operations. In Proceedings of IEEE Conference on Robotics and Automation, 1988. [56] B . J . Gilmore and R . J . Cipra. Simulation of Planar Dynamic Mechanical Systems with Changing Topologies - Part 1: Characterization and Prediction of Kinematic Constraint Changes. Journal of Mechanical Design, Transactions of the ASME, 113:70-76, March 1991. [57] S. Goyal, E . N . Pinson, and F . W . Sinden. Simulation of Dynamics of Interacting Rigid Bodies Including Friction I: General Problem and Contact Model. Engineering with Computers, 10:162-174, 1994. [58] S. Goyal, E . N . Pinson, and F . W . Sinden. Simulation of Dynamics of Interacting Rigid Bodies Including Friction II: Software System Design and Implementation. Engineering with Computers, 10:175-195,1994. [59] A . Gregory, A . Mascarenhas, S. Ehmann, M . Lin, and D. Manocha. Six Degree-of-Freedom Haptic Display of Polygonal Models. In IEEE Proceedings on Visualization, pages 139-146, 2000. [60] P . J . Hacksel and S.E. Salcudean. Estimation of Environment Forces and Rigid-Body Velocities using Observers. In Proceedings of the 1994 IEEE International Conference on Robotics and Automation, pages 931-936, San Diego, C A , May 1994. [61] D.A.Jr. Haessig and B . Friedland. On the Modeling and Simulation of Friction. Transactions of the ASME, Journal of Dynamic Systems, Measurement, and Control, DSC-113:354-362, September 1991. [62] J . K . Hahn. Realistic Animation of Rigid Bodies. Computer Graphics (Proc. SIGGRAPH), 22:299-308, 1988. [63] I. Han and B . J . Gilmore. Impact Analysis for Multiple Body Systems with Friction and Sliding Contact. ASME Journal of Mechanical Design, 115:412-422, September 1993. [64] B . Hannaford and R . J . Anderson. Experimental and Simulation Studies of Hard Contact Force Reflecting Teleoperation. In Proceedings of the 1988 IEEE International Conference on Robotics and Automation, pages 24-29, Scottsdale, A Z , May 1988. [65] B . Hannaford and J . H . Ryu. Time-Domain Passivity of Haptic Interfaces. IEEE Transactions on Robotics and Automation, 18(1): 1-10, February 2002. [66] E . J . Haug, S.C. Wu, and S.M. Yang. Dynamics of Mechanical Systems with Coulomb Friction, Stiction, Impact and Constraint Addition-Deletion - I: Theory. Mechanism and Machine Theory, 21(5):401-406, 1986. [67] V . Hayward and B . Armstrong. A new computational model of friction applied to haptic rendering. In P. Corke and J . Trevelyan, editors, Experimental Robotics VI. Lecture Notes in Control and Information Sciences, volume L N C S 250, pages 404-412. Springer: New-York, 2000. [68] N . Hogan. Impedance Control: A n Approach to Manipulation. Theory, Implementaion, Application. Transactions of the ASME, Journal of Dynamic Systems, Measurement, and Control, 107:1-24, March 1985. [69] N . Hogan. Controlling Impedance at the Man/Machine Interface. In Proceedings of the 1989 IEEE International Conference on Robotics and Automation, pages 1626-1631, 1989. [70] Y . Hurmuzlu and D . B . Marghitu. Three Dimoensional Rigid Body Collisions with Multiple Contact Points. ASME Journal of Applied Mechanics, 62:725-732, 1995. [71] J . D . Hwang, M . D. Williams, and G . Niemeyer. Toward Event-Based Haptics: Rendering Contact Using Open-Loop Force Pulses. In The 12th International Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems, pages 14-21, Chicago, IL, March 2004.  BIBLIOGRAPHY  121  K . W . Jablokow and C S . Bonaventura. Alternate Forms of the Operational Space Inertia Matrix for the Dynamic Simulation of Complex Robot Systems. In Proceedings of the 2003 International Conference on Systems, Man, and Cybernetics, pages 16-21, Washington, D C , October 2003. A. Jain. Unified Formulation of Dynamics for Serial Multibody Systems. Journal of Guidance, Control, and Dynamics, 14(3):531-542, 1991. A. Jain and G . Rodriguez. Diagonalized Lagrangian Robot Dynamics. IEEE Transactions on Robotics and Automation, ll(4):571-584, 1995. D . E . Johnson and P. Willemsen. Six Degree-of-Freedom Haptic Rendering of Complex Polygonal Models. In 11th International Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems, pages 229-235, Los Angeles, Ca, 2003. D . E . Johnson and P. Willemsen. Six Degree-of-Freedom Haptic Rendering of Complex Polygonal Models. In 12th International Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems, pages 18-23, Chicago, IL, 2004. D. Karnopp. Computer Simulation of Slip-Stick Friction in Mechanical Dynamic Systems. Transactions of the ASME, Journal of Dynamic Systems, Measurement, and Control, 107:100-103, March 1985. J.B. Keller. Impact with Friction. ASME Journal of Applied Mechanics, 53:1-4, 1986. A . J . Kelley and S.E. Salcudean. The Development of a Force Feedback Mouse and its Integration into a Graphical User Interface. In Proceedings of the 1994 ASME International Mechanical Engineering Congress and Exposition, volume DSC-55-1, pages 287-294, Chicago, IL, 1994. O. Khatib. A Unified Approach for Motion and Force Control of Robot Manipulators: The Operational Space Formulation. IEEE Transactions on Robotics and Automation, 3(l):43-53, 1987. O. Khatib. Inertial Properties in Robotic Manipulation: A n Object-Level Framework. International Journal of Robotics Research, 13(l):19-36, 1995. S.S. K i m and M . J . Vanderploeg. Q R Decomposition for State Space Representation of Constrained Mechanical Systems. ASME Journal on Mechanism, Transmissions and Automation in Design, 108:183188, 1986. Y . J . K i m , M . C . L i n , and D. Manocha. Incremental Penetration Depth Estimation between Convex Polytopes Using Dual-Space Expansion. IEEE Transactions on Visualization and Computer Graphics, 10(2):152-163, March/April 2004. Y . J . K i m , M . A . Otaduy, M . C . L i n , and D. Manocha. Six-Degree-of Freedom Haptic Display Using Localized Contact Computations. In Tenth Symposium on Haptic Interfaces For Virtual Environment and Teleoperator Systems, pages 209-216, Washington, D C , March 2002. Y . J . K i m , M . A . Otaduy, M . C . L i n , and D. Manocha. Six-Degree-of Freedom Haptic Display Using Incremental and Localized Contact Computations. Presence: Teleoperators and Virtual Environments, 12(3):277-295, June 2003. P.R. Kraus and V . Kumar. Compliant Contact Models for Rigid Body Collisions. In Proceedings of the IEEE International Conference on Robotics and Automation, pages 1382-1387, Albuquerque, N M , 1997. H . M . Lankarani and P.E. Nikravesh. A Contact Force Model with Histeresis Damping for Impact Analysis of Multibody Systems. ASME Journal of Mechanical Design, 112:369-376, 1990. E. Larsen, S. Gottschalk, M . C . Lin, and D. Manocha. Fast Proximity Queries with Swept Sphere Volumes. Technical Report TR99-018, University of North Carolina, Chapel Hill, N C , 1999. D.A. Lawrence. Stability and Transparency in Bilateral Teleoperation. IEEE Transactions on Robotics and Automation, 9(5):624-637, 1993. D.A. Lawrence, L. Pao, A . M . Dougherty, M . A . Salada, and Y . Pavlou. Rate-Hardness: A New Performance Metric for Haptic Interfaces. IEEE Trans Rob Autom, 16(4):357-371, 2000.  BIBLIOGRAPHY  122  K . W . Lilly and D.E. Orin. Efficient Dynamic Simulation of Multiple Chain Robotic Mechanisms. Journal of Dynamic Systems, Measurement, and Control, 116:223-231, 1994. P. Lotstedt. Coulomb Friction in Two-Dimensional Rigid Body Systems. Zeitschrift fur Angewandte Mathematik und Mechanik, 61(12):605-615, 1981. P. Lotstedt. Numerical Simulation of Time Dependent Contact and Friction Problems in Rigid Body Mechanics. SIAM Journal of Scientific Statistical Computing, 42(2):281-296, 1982. P. Lotstedt. Numerical Simulation of Time-Dependent Contact and Friction in Rigid Body Mechanics. SIAM Journal of Scientific Statistical Computing, 5(2):370-393, 1984. S. Mahapatra and M . Zefran. Stable Haptic Interaction with Swistched Virtual Environments. In Proceedings of the 2003 IEEE International Conference on Robotics and Automation, Taipei, Taiwan, September 2003. M . Mahvash and V . Hayward. High Fidelity Haptic Synthesis of Contact With Deformable Bodies. IEEE Computer Graphics and Applications, 24(2):48-55, 2004. T. Maneerwan, B . Hannaford, D . Sorti, and M . Ganter. Haptic Rendering for Internal Content of an Implicit Object. In ASME Winter Annual Meeting Haptics Symposium, Nashville, T N , November 1999. N . K . Mani, E . J . Haug, and K . E . Atkinson. Application of Singular Value Decomposition for Analysis of Mechanical System Dynamics. ASME Journal on Mechanism, Transmissions and Automation in Design, 107:83-87, 1985. D.B. Marghitu and Y . Hurmuzlu. Three Dimensional Rigid Body Collisions with Multiple Contact Points. The ASME Journal of Applied Mechanics, 62(3):725-732, 1995. W.R. Mark, S.C. Randolph, M . Finch, J . M . Van Verth, and R . M . Taylor II. Adding Force Feedback to Graphics Systems: Issues and Solutions. In Haptic Virtual Reality for Blind Computer Users, Assistive Technologies, pages 92-99, 1998. M . Mason. Mechanics and Planning of Manipulator Pushing Operations. Robotics Research, 5(3):53-71, 1885.  International Journal of  M . Mason and Y . Wang. On the Inconsistency of Rigid-Body Frictional Planar Mechanics. In IEEE 1988 International Conference on Robotics and Automation, pages 524-528, Philadephia, P A , 1988. W . A . McNeely, K . D . Puterbaugh, and J.J. Troy. Six Degree-of-Freedom Haptic Rendering using Voxel Sampling. In Proceedings of the 1999 AMC SIGGRAPH, pages 92-99, 1999. A. Miller and H.I. Christensen. Implementation of Multi-rigid-body Dynamics within a Robotic Grasping Simulator. In Proceedings of the IEEE International Conference on Robotics and Automation, pages 2262-2268, Taipei, Taiwan, September 2003. B. E. Miller, J . E . Colgate, and R . A . Freeman. Guaranteed Stability of Haptic Systems with Nonlinear Virtual Environments. IEEE Transactions on Robotics and Automation, 16(2):712-719, December 2000. J.K. Mills and C . V . Nguyen. Robotic Manipulator Collisions: Modeling and Simulation. The AMSE Jounral of Dynamic Systems,, Measurement, and Control, 114:650-658, December 1992. B. Mirtich. Impulse-Based Dynamic Simulation of Rigid Body Systems. PhD thesis, University of California, Berkeley, C A , 1996. B. Mirtich. V-Clip: Fast and Robust Polyhedral Collision Detection. ACM Transactions on Graphics), 17(3):177-208, 1998. B. Mirtich. Timewarp Rigid Body Simulation. In Proceedings of the 2000 ACM SIGGRAPH, New Orleans, L A , 193-200 2000.  BIBLIOGRAPHY  123  [110] B . Mirtich and J . Canny. Impulse-based Dynamic Simulation. In K . Goldberg, P. Halperin, J.-C. Latombe, and R. Wilson, editors, Workshop on Algorithmic Foundations of Robotics, pages 407-418. A.K.Peters, Boston, M A . , 1994. [ I l l ] M . Moore and J . Wilhelms. Collision Detection and Response for Computer Animation. Computer Graphics (Proceedings of ACM 1988 SIGGRAPH Conference), 22:289-298, 1988. [112] J.J. Moreau. Unilateral Contact and Dry Friction in Finite Freedom Dynamics. In J . J . Moreau and P.D. Panagiotopoulos, editors, Nonsmooth Mechanics and Applications. CISM Courses and Lectures, volume 302, pages 1-82. International Centre for Mechanical Sciences, Springer-Verlag, 1988. [113] A . Nahvi, J . M . Hollerbach, R. Freier, and D.D. Nelson. Display of Friction in Virtual Environments based on Human Finger Pad Characteristics. In Proceedings of the ASME Dynamics Systems and Control Division, volume DSC-64, pages 179-184, 1998. [114] A . Nahvi, D.D. Nelson, J . M . Hollerbach, and D.E. Johnson. Haptic Manipulation of Virtual Mechanisms from Mechanical C A D Designs. In Proceedings of the 1998 IEEE International Conference on Robotics and Automation, pages 375-380, Leueven, Belgium, May 1998. [115] P.L. Nash. A New First Integral for a Binary Rigid Body Collision of Arbitrarily Short Duration. Jounral of Mathematical Physics, 40(6):2816-2829, 1999. [116] N . Orlandea, M . A . Chace, and D.A. Calahan. A Sparsity-Oriented Approach to the Dynamic Analysis and Design of Mechanical Systems - Part I. Journal of Engineering in Industry, 99:773-784, 1977. [117] M . A . Otaduy and M . C . Lin. Sensation Preserving Simplification for Haptic Rendering. J 4 C M Transactions on Graphics. Special issue: Proceedings of the ACM SIGGRAPH 2003, 22(3):543-553, 2003. [118] J.-S. Pang and J.C. Trinkle. Complementarity Formulations and Existence of Solutions of Multi-Body Contact Problems with Coulomb Friction. Mathematical Programming, 73:199-226, 1994. [119] J . G . Park and G . Niemeyer. Haptic Rendering with Predictive Representation of Local Geometry. In Proceedings of the 12th International Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems, pages 331-338, Chicago, IL, 2004. [120] F. Pffeifer and C. Glocker. Multibody Dynamics with Unilateral Contacts. Wiley, New York, 1996. [121] J.C. Piatt and A . H . Barr. Constraint Methods for Flexible Models. Computer Graphics (Proc. SIGGRAPH), 22:279-288, 1988. [122] D. Redon, A . Kheddar, and S. Coquillart. Gauss' Least Constraints Principle and Rigid Body Simulations. In Proceedings of IEEE International Conference on Robotics and Automation, pages 517-522, Washington, D C , 2002. [123] M . Renz, C. Preusche, M . Potke, H.P. Kriegel, and G . Hirzinger. Stable Haptic Interaction with Virtual Environments Using an Adapted Voxmap-Pointshell Algorithm. In Proceedings of the Eurohaptics Conference, pages 149-154, 2001. [124] C. Richard and M . Cutkosky. The Effects of Real and Computer Generated Friction on Human Performance in a Targeting Task. In Proceedings of the ASME/IMECE Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems, Nashville, T N , 2000. [125] L . B . Rosenberg and B.D. Adelstein. Perceptual Decomposition of Virtual Haptic Surfaces. In Proceedings of the 1993 IEEE Symposium on Research Frontiers in Virtual Reality, pages 46-53, San Jose, C A , 1993. [126] D. Ruspini and 0 . Khatib. Dynamic Models for Haptic Rendering Systems. Kinematics: ARK98, pages 523-532, Strobl/Salzburg, Austria, 1998.  In Advances in Robot  [127] D. Ruspini and O. Khatib. Collision/Contact Models for Dynamic Simulation and Haptic Interaction. In Proceedings of the 9th International Symposium on Robotics Research, pages 185-195, Snowbird, U T , 1999.  BIBLIOGRAPHY  124  [128] D . Ruspini and 0 . Khatib. Haptic Display for Human Interaction with Virtual Dynamic Environments. Journal of Robotic Systems, 18(2):769-783, 2001. [129] D . Ruspini, K . Koralov, and 0 . Khatib. The Haptic Display of Complex Graphical Environments. In Proceeding of the SIGGRAPH 97, pages 345-352, Los Angeles, C A , 1997. [130] S.E. Salcudean. Control for Teleoperation and Haptic Interfaces. In B Siciliano and K . P . Valavanis, editors, Lecture Notes in Control and Information Sciences 230 - Control Problems in Robotics and Automation, pages 51-65. Springer-Verlag, 1997. [131] S.E. Salcudean and T . D . Vlaar. On the Emulation of Stiff Walls and Static Friction with a Magnetically Levitated Input/Output Device. ASME J Dyn Meas Contr, 119:127-132, March 1997. [132] S.E. Salcudean, N . M . Wong, and R . l . Hollis. Design and Control of a Force-Reflecting Teleoperation System with Magnetically Levitated Master and Wrist. IEEE Transactions on Automatic Control, 6:844858, December 1995. [133] J . K . Salisbury and C. Tarr. Haptic Rendering of Surfaces Denned by Implicit Fuctions. In Proceedings of the ASME Dynamics Systems and Control Division, volume DSC-61, pages 61-67, 1997. [134] H . Schmidl and V . J . Milenkovic. A Fast Impulsive Contact Suite for Rigid Body Simulation. IEEE Transactions on Vizualization and Computer Graphics, 10(2):189-197, March/April 2004. [135] K . Shimoga. A Survey of Perceptual Feedback Issues in Dexterous Telemanipulation: Part I. Finger Force Feedback. In Proceedings of the 1993 IEEE Annual Virtual Reality International Symposium, pages 263-270, Seattle, W A , 1993. [136] M . R . Sirouspour, S.P. DiMaio, S.E. Salcudean, P. Abolmaesumi, and C. Jones. Haptic Interface Control - Design Issues and Experiments with a Planar Device. In Proc 2000 IEEE Int Conf Rob Autom, pages 789-794, San Francisco, Ca, 2000. [137] C . E . Smith. Predicting Rebounds using Rigid-Body Dynamics. Journal of Applied Mechanics, 58:754758, September 1991. [138] W . Son, K . K i m , and N . M . Amato. A n Interactive Generalized Motion Simulator (GMS) in an ObjectOriented Framework. In Proceedings of the IEEE Computer Animation 2000, pages 176-182, Phyladelphia, P A , May 2000. [139] W . Son, K . K i m , N . M . Amato, and J.C Trinkle. Interactive Dynamic Simulation using Haptic Interaction. In Proceedings of the 2000 IEEE International Conference on Robotics and Automation, pages 145-150, San Francisco, C A , May 2000. [140] W . Son, K . K i m , N . M . Amato, and J.C Trinkle. A Generalized Framework for Interactive Dynamic Simulation for MultiRigid Bodies. IEEE Transactions on Systmes, Man, and Cybernetics - Part B: Cybernetics, 34(2):912-924, April 2004. [141] W . Son, J.C Trinkle, and N . M . Amato. A n Interactive Generalized Motion Simulator (I-GMS) in an Object-Oriented Framework. In Proceedings of the 2001 IEEE International Conference on Robotics and Automation, pages 3789-3794, Seoul, Corea, May 2001. [142] D . E . Stewart. Rigid Body Dynamics with Friction and Impact. SIAM Review, 42(l):3-29, 2000. [143] D . E . Stewart and J.C. Trinkle. A n Implicit Time-Stepping Scheme Rigid Body Dynamics with Inellastic Collisions and Coulomb Friction. The International Journal of Numerical Methods in Engineering, 39(15):2673-2691, 1996. [144] D . E . Stewart and J.C. Trinkle. Convergence of a Time-Stepping Scheme for Rigid Body Dynamics and Resolution of Painleve's Paradoxes. Archives of Rational Mechanics, 1998. [145] W . J . Stronge. Rigid Body Collisions with Friction. Proceedings of the Royal Society of London, A43L169181, 1990.  BIBLIOGRAPHY  125  146] F . E . Udwadia and R . E . Kabala. Explicit Equations of Motion for Mechanical Systmes with Nonideal Constraints. Transactions of the ASMEj Journal of Applied Mechanics, 68:462-467, May 2001. 147] G . van der Bergen. Proximity Queries and Penetration Depth Computation on 3D Game Objects. In Game Developers Conference, 2001. 148] A . van der Schaft. L2-gain and Passivity Techniques in Nonlinear Control. Springer, London, U . K . , 2000. 149] S. Vedula and D. Baraff. Force feedback in interactive dynamic simulation. In Proceedings of the First PHANTOM Users Group Workshop, pages 54-57, 1996. 150] S. Walker and J . K . Salisbury. Large haptic topographic maps and the proxy graph algorithm. Proceedings of the ACM 2003 Symposium on Interactive 3D Graphics, pages 83-92, 2003.  In  151] Y . Wang and M . T . Mason. Two-Dimensional Rigid-Body Collisions with Friction. ASME Journal of Applied Mechanics, 59:635-642, 1992. 152] A . Witkin and D. Baraff. Physically Based Modeling: Principles and Practice. S I G G R A P H '97 Course Notes. 153] S.C. Wu, S.M. Yang, and E . J . Haug. Dynamics of Mechanical Systems with Coulomb Friction, Stiction, Impact and Constraint Addition-Deletion - II: Planar Systems. Mechanism and Machine Theory, 21(5):407-416, 1986. 154] S.C. Wu, S.M. Yang, and E . J . Haug. Dynamics of Mechanical Systems with Coulomb Friction, Stiction, Impact and Constraint Addition-Deletion - III: Spatial Systems. JSME International Journal, Series C: Mechanical Systems, Machine Elements and Manufacturing, 43(5):387-393, June 2000. 155] C B . Zilles and J . K . Salisbury. A Constraint-based God Object Method for Haptic Display. In ASME Haptic Interf Virt Envir Teleop Syst, pages 146-150, Chicago, IL, 1994.  Appendix A  Proof of implicit elimination of constraint overdeterminancy This appendix proves that, if J is row rank deficient, then: c  jj  (jeD" Jj) 1  J  ]  where J is full row rank, and j j = [j£  = J „ (jnBT  c  JD'  1  1  (A.l)  J  n  jj\.  n  Proof: Let  J  c  rank(J ) r  be given by Jj  =  [j„  J j Y  w  n  e  r  Jn  e  is full row rank, i.e., rank(J )  = rank(J )  n  c  = n, and  = r. By elementary row operations:  In  0  A  I  r  (A.2)  Cfr  0  where A € TZ , and I„ and I are identity matrices of rank n and r, respectively. rxn  r  Then:  /r (JcV-'J?)'  =  V  -| -1  In A  0  I  •  Jn  D  1  In  0  -A  I  r  Cfn  0  o  Jl  0  r  D"  In  A  0  I  T  In  1  Jl  *  T  r  —A  T  T  0  I,  t In  A  JnD^JI  0  In  0  0  I  0  0  A  I  T  r  126  r  (A.3)  127  The last algebraic manipulation is based on the fact that:  (xx ) = (x ) xt. T  f  r  f  (A.4)  To show that Equation (A.4) holds, let the SVD of X be given as X = U E V . Then: T  (XX ) T  f  =  (UEEU ) T  =UE'S'U  f  = UE'V VE'U  T  T  T  = (X ) X , T  t  t  (A.5)  where: 1 Si  E'  (A.6)  0 and E i , • • •, E are the singular values of X . n  jj  (j D c  _  1  Jj)^  Jc  can now be computed using Equations (A.2) and (A.3):  jI{JcV- jJ)^ l  J  = Jc  2  c  In  A  0  I  T  JnB-'jI  0  In  0  0  0  A  I  r  (^D-ijJ)"  1  0  Jc  =  r  \Jn  0 — J  n  (tTnD  J  n  )  J  n  (A.7)  Appendix B  Weighted generalized inverses for linkage manipulation This appendix briefly overviews weighted generalized inverses applied to robotics, following the in-depth discussion of their properties and their correct application in robotics introduced in [46].  B.l  Weighted generalized inverses of matrices  Consider a linear transformation A that maps a real vector space X with metric M , into a real vector space U with metric M : u  u = Ax.  (B.l)  If A is singular, then a solution x to Equation (B.l) may be derived as: s  x =A#u,  (B.2)  s  where A * is called the weighted generalized inverse of A [20] and is computed by: A  #  = M~ C l  T  (CM^C ) T  _ 1  (F M„F) T  _ 1  F M . T  U  (B.3)  In Equation (B.3), A = F C is a full rank factorization of A , i.e., F is full row rank and C is full column rank. A * is a weighted generalized inverse of A because, for any symmetric positive definite weights M , and M „ , it satisfies the following properties [20]: 1. A A * A = A 2. A#AA# = A# 3. ( M A A # ) U  T  = M AA* n  128  129  B.2 Weighted generalized inverses applied to robotics  4. (M A*A)  T  X  = MsA^A  Additionally, A# is unique, (A#)  = A and (A#)  #  = (A )  T  r  [20].  #  The weighted generalized inverse A * computes the minimum M^-norm: WM,  :  > ,  =< ^ . M A  (B.4)  of all solutions that produce the minimum M -least squares error: u  |u - A X „ | M  :=< (u - A x ) , M  u  s  (u - Ax,) > ,  u  (B.5)  where <, > denotes the inner product on a metric space. As discussed in [46], the weighted generalized inverse is invariant to the choice of metric M  x  on X if A = C  has full column rank: A  #  = M "  ^ (CM.-'C )"  1  7  and is invariant to the choice of metric M A*  u  = M~ A  1  1  (AMJ A )"  T  1  T  (B.6)  1  on U if A = F is full row rank:  = (F M„F) T  _ 1  F M T  = (A M„A) T  U  A M .  _ 1  (B.7)  T  U  If A is invertible, then its weighted generalized inverse is invariant to both metrics:  A  B.2  #  = (A M A) T  U  _ 1  A M„ = A T  _ 1  M~ A 1  - T  A" M r  = A" . 1  u  (B.8)  Weighted generalized inverses applied to robotics  A linear transformation A (that maps X into U) and its transpose A  T  (that maps W into T):  u  =  Ax  (B.9)  t  =  A w  (B.10)  P:=<u,v>  (B.ll)  T  and the scalar product:  are called a dual system [46] if all terms in P have the same physical units. The importance of dual systems lies in the fact that their solution is simplified by the fact that M " and M " are metrics on T and W whenever 1  M  x  and M  u  1  and metrics on X and U. Consequently, the dual system may be solved using the weighted  B.2 Weighted generalized inverses applied to robotics ,  ...  130  generalized inverse A # and its transpose ( A * ) : T  x  =  A  w  =  (A#) t.  #  u  (B.12) (B.13)  T  Until now, the discussion of linear transformations and their weighted generalized inverses has been general, but it has implicitly relied on the metric structure of X and U and the inner products defined on these spaces through the metrics M , and M . The relevance of this theory to robotics (and to virtual linkage manipulation) u  is apparent when considering suitable choices of metrics on the linear spaces of interest. In robotics, the linear transformations of interest are the Jacobian matrix of a linkage J and its transpose J . The linkage Jacobian maps the space Q of joint velocities to the space V of twists (rigid body velocities): T  x = Jq.  (B.14)  In Equation (B.14), q is the vector of configuration (linear and angular joint) velocities of the linkage, and x = (v  T  w) T  T  is the twist at the point of interest, with v S TZ being the linear velocity at the point of 3  interest and cu S TZ being the angular velocity of the link on which the point of interest lies. The transpose 3  of the linkage Jacobian maps the space of wrenches F to the space of configuration (joint) torques T: T = J F.  (B.15)  T  In Equation (B.15), T is the vector of configuration (joint) torques, and F = ( f point of interest, with f £ K  3  T  r) T  T  is the wrench at the  being the force at this point and T € TZ being the torque. According to the 3  discussion above, J , J , and the scalar product (representing power): T  P:=<x,F>  (B.16)  form a dual system, i.e., x and F (and q and T) are said to be duals of each other. This means that finding suitable metrics on Q and V is sufficient for inverting both Equation (B.14) and Equation (B.15). Therefore, only the mapping of configuration velocities to twists is considered in the following. When choosing metrics on these two spaces, two conditions must be satisfied in order for the metrics to have physical significance: 1. The inner products on the two spaces must have the same units in all the terms (i.e., they must be physically consistent). 2. The metric must be gauge invariant (i.e., invariant to coordinate, or rigid body, transformations). As shown in [46], kinetic energy metrics satisfy these two properties. This thesis is concerned with enabling realistic and unrestricted (i.e., from any user-selected link and through singularities) haptic manipulation of arbitrary virtual linkages. Therefore, two kinetic energy metrics are chosen  B.2 Weighted generalized inverses applied to robotics  131  to compute the weighted generalized inverse of the Jacobian of the virtual tool at the user's hand, according to: J# = D ^ C  T  (CD;/C ) T  _  1  ( F M  W  F  T  )  -  1  FM . T  hl  (B.17)  In Equation (B.17), ~D is the configuration space inertia of the virtual tool and M u is the inertia of the link vt  held by the user, computed in the world coordinate frame. J * is used to represent the contact interactions of the virtual tool at the user's hand in a coordinate-invariant manner when the user holds a link whose motion is restricted by the linkage to which it belongs. Example such interactions include a user rotating a virtual crank or opening a virtual door. Note that, if the Jacobian of the virtual tool at the user's hand is full column rank, i.e., J/» = C , its weighted generalized inverse reduces to the dynamically consistent inverse of the Jacobian [80]:  J^D-^nJfcD- ^)- . 1  1  (B.18)  When the Jacobian looses rank, i.e., when the linkage topology restricts the rigid body motion at the user's hand, the inertia matrix of the link held by the user M M selects the rigid body motion that best approximates, in the M^-least-squares sense, the rigid body motions not in the range of J/,.  Appendix C 1  Equality of null spaces of  rri  and  This appendix proves that the SVD of the virtual tool Jacobian computed at the user's hand  can be used to  derive the singular directions of the inverse of the virtual linkage operational space inertia at the user's hand Ar . 1  Let the SVD of  be given by: (C.l)  J^exe — U 6 x 6 £ 6 x < f V j . x d  Then:  hL  A  6  -  U  6x6S xdVj 6  x  D„^  d  Vd  x d  (C.2)  S r>T d UTTT 6 . x 6  x 6  Since D ^ is symmetric positive definite, it can be reduced to the diagonal form A through a suitable rotation 1  R: 0  01  A =  RD^R 0  Moreover, Y = R V = (yi  ...  0  (C.3)  T  adj  yd) is orthogonal (since V R R V T  T  = I), and the product V D , V is r  1  t  t  diagonal and equal to A : V  T  D ^ V  V R ARV T  T  LA (aj yi  0  yi 0  \yl)  0  • • • yd  ad)  (yf A 1  1  aiyi  •• •  \y J T  d  132  = A.  ay d  d  (CA)  133  Furthermore, let rank (J^) = r, i.e.  ^Vxr ^6xd  °rx(d-r)  0(6-r)xr  where  S  R X R  .  (C.5)  — 0(6-r)x(d-r)  is the r-dimensional diagonal matrix having the r non-zero singular values S i , . . . , E of r  on its  main diagonal. Then, from Equation (C.2) and Equations (C.4) and (C.5), it follows that:  (C.6) ®(d—r)xr  In Equation (C.6), A  TrXT  0(d-r)x(d-r)  is the diagonal matrix having a i ,  aiS S'AE'  , a on its main diagonal, and: r  2  0  =  (C.7) o £  2  r  Finally, after substitution from Equation (C.7), Equation (C.2) becomes:  A/Jve  with  U  6  X  6  orthogonal, and £ A  6 x 6  —  U 6 X 6 S A  6  X  6  U ^  X  6  ,  (C.8)  diagonal and having the first r diagonal elements non-zero. In other words,  Equation (C.8) gives the S V D of A r . 1  Appendix D  The dynamics and control of the planar haptic interface The planar haptic interface used in the experiments described in this thesis is shown in Figure D . l .  The  interface comprises two identical pantograph mechanisms mounted in parallel. As depicted in Figure D.2, each pantograph consists of four carbon fibre links and five rotational joints. The pantograph endpoints are connected through a link that freely rotates with respect to each mechanism. This mechanical arrangement endows the haptic interface with three DOFs, translation along any horizontal axis and unlimited rotation about the vertical axis.  Figure D . l : The twin pantograph planar haptic interface.  D.l  Interface actuation and dynamics  Joint angles 9\ and #2 are measured by digital encoders with a resolution of 0.0225 degrees (16,000 counts per revolution). These joints are actuated by 90W D C motors and are considered torque sources. A l l other joints are passive. The four idential motors that actuate the haptic interface are driven by constant current sources  134  D.l Interface actuation and dynamics  135  Figure D.2: The parameters of one pantograph linkage. The base joints are shown in black. The pantograph endpoint P is the joint between the two not grounded links. e  and produce torques T proportional to the armature current i through the motor torque constant kta  T = ki . t  (D.l)  a  The motors have a motor torque constant k — 0.0525NmA , a continuous current rating of 2A, and a peak -1  t  current rating exceeding 10A. They provide a peak force of 23N and a peak torque of 115Ncm at the ahndle. Peak accelerations of 6g have been obtained at half the maximum coil current. Joint velocities #1 and # and user-applied wrenches F/j are estimated from joint angle measurements and 2  applied motor torques using a force observer [60] and accurate interface dynamics [136]. In terms of the actuated joint variables 0 = (9i  0 ) , the configuration space dynamics of each pantograph are derived as [136]: T  2  D (0) f)> C ( M ) 0 = T - J j f , p  where D T  p  p  p  p  p  (D.2)  and C are the configuration space inertia and Christoffel matrices of one pantograph, respectively, p  is the vector of applied actuator torques, J  is the pantograph Jacobian computed at its endpoint, and f  p  p  is the external force applied at its endpoint. The configuration space dynamics of the individual pantographs are combined into the operational space dynamics of the haptic interface [136]: Mi  d h  + C x =F +3^T d  h  h  = F + u. h  (D.3)  In Equation (D.3), x^,, x^, x/j are the body acceleration, velocity, and position of the user's hand (i.e., device handle), respectively, interface, respectively,  and C  d  are the operational space inertia and Christoffel matrices of the haptic  is the user-applied wrench, J is the Jacobian matrix of the interface, T is the vector e  D.2 Interface control  136  of applied actuator torques, and u is the operational space control signal. Since the haptic interface has four motors and three DOFs, the actuation redundancy is used to minimize the force along the link connecting the two pantograph endpoints: (D.4) where j £ = J (J^J ) 0  0  is the right pseudo-inverse of J J. In turn, Jj maps the forces applied by the endpoints  of the two pantograph mechanisms to the wrench applied by the haptic interface:  fyl  (D.5)  1x2  D.2  Interface  control  The interface is controlled by considering that the haptic interaction within virtual environments is a teleoperation task. According to this view, the haptic interface is the master and the virtual tool is the slave. High performance interaction [89,130] is achieved by controlling the haptic device using a four channel teleoperation architecture [89]. The advantage of the four channel architecture is that it allows the virtual environment to be designed independently from the haptic interface control. The four channel architecture is illustrated in Figure D.3 for one D O F interaction. In this figure, fh and f  env  are the user-applied and the environment force, respectively, ih and ±h  are the user's hand and the  vt  virtual tool velocities, respectively, and Zh, Z , Z , d  vt  and Z  env  are the impedances of the user's hand, the  device, the virtual tool, and the virtual environment, respectively. Furthermore, C\ and CA are the position channel controllers, Ci and C 3 are the force channel controllers, and C& and C t are the position controllers of v  the local device and the virtual tool, respectively. Transparent interaction is difficult to achieve if the master and the slave are kinematically and dynamically dissimilar [130]. Therefore, an impedance controller shapes the impedance of the haptic device to match the desired impedance of the virtual tool, i.e., Z  d  Equation (5.21).  = Z. vt  Impedance shaping is performed according to  D.2 Interface control  137  Virtual environment Virtual tool & local control (slave) Communication channel Haptic device (master) Human operator Figure D.3: The four channel architecture for 1 D 0 F haptic interaction within virtual environments. Simulated blocks are shaded.  In this thesis, the following controller parameters were used:  c  d  =  35+^2  =  35+A22  =  c  =  1  c  =  1  2  3  = The adaptive damping B dp = Kbfenv + B i„ a  B  min  = Okg/s) [132].  s  Ci  m  Ct v  -c  d  is added to C t in order to reduce chattering (Ki = 60s/m, v  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0064992/manifest

Comment

Related Items