UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Properties of robot forward dynamics algorithms with applications to simulation Cloutier, Benoit P. 1995

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

Item Metadata

Download

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

Full Text

P R O P E R T I E S OF R O B O T F O R W A R D D Y N A M I C S A L G O R I T H M S W I T H A P P L I C A T I O N S TO S I M U L A T I O N By Benoit P. Cloutier B. Eng. Ecole Polytechnique de Montreal A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E DEGREE O F M A S T E R OF SCIENCE in T H E FACULTY OF GRADUATE STUDIES COMPUTER SCIENCE We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA March 1995 © Benoit P. Cloutier, 1995 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Computer Science The University of British Columbia 2366 Main Mall Vancouver, Canada V6T 1Z4 Date: Apr,'/, mr Abstract This thesis presents issues related to the dynamic simulation of robot manipulators. The numerical simulation problem is usually treated as two separate problems: the forward dynamics problem for computing system accelerations, and the numerical integration problem for advancing the state in time. Using a compact and. unifying notation, ex-isting methods for computing robot forward dynamics are presented and compared. It is shown that the articulated-body method [19, 20] (ABM) is better suited to deal with certain types of numerical problems than the composite rigid-body method [58] (CRBM). Simulation results are presented and the practical implications of these results are con-sidered. In particular, it is shown that the fastest forward dynamics methods are not necessarily best when considered in conjunction with the popular adaptive step-size in-tegration methods. Finally, this thesis also reviews issues related to the difficult task of incorporating contact simulation in the robot simulation process. ii Table of Contents Abstract ii List of Figures vi 1 Introduction 1 1.1 Thesis concerns 1 1.2 Thesis Overview 2 2 Fundamentals of Robot Manipulator Simulation 3 2.1 General notation conventions 3 2.2 Robot model 5 2.3 The formulation problem 7 2.4 Dynamics of a single rigid body 10 2.4.1 Spatial vectors 10 2.4.2 Rigid body transformation operator 11 2.4.3 Coordinate transformation operator 12 2.4.4 Newton and Euler's equations of motion 13 2.5 The inverse dynamics problem 14 2.5.1 Generalized spatial operators 16 2.6 The motion integration problem 18 3 Solving the Forward Dynamics Problem 20 3.1 0(n 3) Algorithms: inertia matrix inversion methods 20 iii 3.1.1 The composite rigid-body method 21 3.1.2 Orthogonal complement, based methods 22 3.2 0(n 2) Algorithms: conjugate gradient methods 23 3.3 0(n) Algorithms: linear recursive methods 24 3.3.1 The articulated-body method 25 3.4 Summary and comparisons 27 4 Numerical Properties of Forward Dynamics Algorithms 29 4.1 Review of issues in numerical computation 29 4.2 Numerical properties of the C R B M 30 4.3 Numerical properties of the A B M 32 4.4 Impact on motion integration: formulation stiffness 34 5 Implementation and Experimental Verification 35 5.1 PRRRobot simulation software 35 5.1.1 Assumptions and definitions 35 5.1.2 Forward dynamics test program 37 5.1.3 Off-line simulation program 38 5.1.4 Animation program 40 5.1.5 Implementation issues 41 5.2 Results and observations: simulation of ill-conditioned systems 44 5.2.1 Performance of the forward dynamics 44 5.2.2 Performance of the motion integration 48 5.3 More simulation results 50 5.3.1 "Bad" implementations of the A B M 50 5.3.2 Conservation of energy 50 5.3.3 Simulation of long chains 53 iv 6 Extensions: Contact Simulation 57 6.1 Notation 57 6.2 Contact detection 59 6.3 Friction modelling 60 6.4 Impact simulation 61 6.4.1 Penalty methods 61 6.4.2 Analytical methods 62 6.4.3 Finite element methods 65 6.5 Simulation of continuous contact 66 6.5.1 Penalty methods 66 6.5.2 Analytical methods 66 6.5.3 Impulse based methods 68 7 Conclusion 69 Bibliography 71 v List of Figures 2.1 Notation for link k 6 5.1 N-link PRRRobot manipulator 37 5.2 PRRRobot off-line simulator interface and options 39 5.3 Effect of configuration on the condition number of 2 link chain. Note that the condition number is largest at the singularities in this case. The relative difference for joint 1 accelerations (Ar#i) is also shown; results for joint 2 accelerations are similar 45 5.4 Effect of link length 46 5.5 Comparison of calculated joint 1 accelerations (with l[2]=m[2]=le5). Note the smoothness of the A B M solution 47 5.6 Number of forward dynamics evaluations required to simulate 10 seconds of motion for various two link chains of increasing condition number. . . 49 5.7 Joint accelerations for "good" and "bad" implementations of the A B M (the smooth curve is the "good" implementation of the ABM) 51 5.8 Evolution of the total energy. 52 5.9 Summary of results for the simulation of long regular chains 54 5.10 Simulation of long chains 55 6.1 Contact simulation notation 58 vi Chapter 1 Introduction Computer simulation is playing an increasingly important role in the field of robotics. Simulation of robot dynamics has been used in off-line tasks such as robot design testing, robot programming, and controller validation. Recently, the need for real-time perfor-mance has become important for on-line applications such as predictive displays for time-delayed teleoperation [53, 9], virtual environments for operator training and development of advanced robot control schemes. In robotics, the simulation problem requires that the joint trajectory of a robot be determined given prior knowledge of the applied forces and the initial robot state. The simulation problem can be divided into two subproblems: • The forward dynamics problem: computing the joint accelerations given the actu-ator forces/torques and the current robot joint positions and velocities. • The motion integration problem: computing the joint trajectory (joint positions and velocities) given the joint accelerations and initial conditions. 1.1 Thesis concerns In recent years, the need for increased performance in robot simulation has led to the development of many different algorithms for solving the forward dynamics problem. The algorithms range in computational complexity from 0(n 3) to 0(n) where n is the number of degrees of freedom of the robot. Detailed counts of the total number of operations 1 Chapter 1. Introduction 2 per forward dynamics evaluation is a common means of comparing the performance of the different algorithms. However, the computational complexity of a forward dynamics method is only half the story; accurate numerical integration of the resulting acceleration is an important factor in the total efficiency of a simulation method. The overall efficiency of a simulation procedure is best measured in terms of the computational effort necessary to achieve a desired accuracy in the solution. In this light, it will be shown how the numerical properties of forward dynamics algorithms can play an important role in the selection and implementation of algorithms for the simulation of robot manipulators. These issues will be first discussed in a theoretical framework and verified experimentally in a number of computer simulations. Another important issue in the development of robot simulators is the need to simulate contacts of the robot with its environment. Recently, the field of contact dynamics has attracted a significant amount of research efforts. In this thesis, the problems related to incorporating such capabilities into a robot simulation system will be discussed. 1.2 Thesis Overview This thesis is organized as follows. Chapter 2 presents the problem statement, assump-tions and notation conventions. Chapter 3 gives an overview of the algorithms for solving the forward dynamics problem. Chapter 4 explores the numerical properties of two of the most commonly used forward dynamics methods: the composite-rigid body method (CRBM) and the articulated body method (ABM). Chapter 5 presents some numeri-cal results obtained with regards to the forward dynamics problem and the simulation problem. Chapter 6 discusses possible extensions to include contact simulation. Finally, chapter 7 presents some concluding remarks. Chapter 2 Fundamentals of Robot Manipulator Simulation This chapter introduces and defines some of the basic concepts that will be explored within the context of this thesis. The first section of the chapter presents an overview of the general notation conventions. As will be demonstrated, a good notation plays a key role in the formulation of clear and intuitive equations of robot dynamics. The second section presents the robot model and some of its associated notation. The third section looks at the problem of formulating the equations of motion of a robot manipu-lator. The formulation problem is denned and some of the key high level issues related to this problem are introduced. The fourth section of this chapter introduces the spatial notation through the study of the dynamics of a single rigid body. The fifth section deals with multibody dynamics. In particular, a solution to the inverse dynamics problem is presented (the forward dynamics problem is reviewed in detail in chapter 3). Finally, the last section presents a brief overview of the motion integration problem. Notation is in-troduced gradually throughout each section of this chapter. This notation is the result of a combination of several previously published notations including Craig's, Featherstone's and Jain's [15, 20, 32]. 2.1 General notation conventions The following notation conventions are observed throughout this thesis: • Uppercase script letters (Latin) are used to designate points, sets of points and coordinate systems (e.g. Og designates an origin point on a rigid body B). 3 Chapter 2. Fundamentals of Robot Manipulator Simulation 4 • Lowercase letters (Latin or Greek) are used to represent vector and scalar quantities. • Uppercase letters are used to represent linear operators (matrix quantities). • Overlined letters designate a three dimensional vector quantity or linear operator (e.g. OB is the 3 x 1 dimensional angular velocity vector of a body B). • Letters topped with a caret symbol ( * ) designate six dimensional spatial vectors or linear operators (the spatial notation will be defined in a later section). • Underlined letters designate n-dimensional generalized vectors or linear operators (e.g. 9 is the n x 1 dimensional vector of joint angular velocties of a n degree of freedom robot manipulator). • Overlines, carets and underlines can be combined to designate 3n or 6n dimensional block entities (e.g. v is the 6n x 1 dimensional vector of link spatial velocities of an n degree of freedom robot manipulator, where ujt is the 6 x 1 dimensional spatial velocity of link k). • Trailing (right) subscripted numbers or lowercase letters are used to designate com-ponents of a vector or linear operator (e.g. J4JJ is component (ij) of the linear operator A). • Trailing subscripted uppercase letters can be used to provide an additional descrip-tion of a quantity (e.g. / T is the generalized total spatial force acting on the links of a robot manipulator). • Trailing superscripts are used to define special operators (e.g. A"1 is the inverse of A and AT is the transpose of A). Chapter 2. Fundamentals of Robot Manipulator Simulation 5 • Leading (left) superscripted letters or numbers are used to designate the coordinate frame with respect to which the components of a quantity are expressed. A quantity without a leading superscript is said to be coordinate free. However, we may often assume that lVi = tij in order to maximize clarity (such usage should be clear from context). • Leading subscripts are not used. • The identity matrix is designated by 1 and the null matrix is designated by 0. The dimensions of these quantites should be clear from the context. • The 3 x 3 skew-symmetric matrix (cross-product matrix) corresponding to a 3 x 1 vector v = [ vi v-i vj, ] T is denoted v. 0 V3 -v2 V = -t>3 0 U 2 -Vi 0 In summary, lowercase and uppercase letters are used to distinguish between vector quantities and linear operators. Overlines, carets and underlines are used to distinguish between three dimensional, six dimensional (spatial) and n-dimensional (generalized) entities. 2.2 Robot model The robot manipulators considered within the context of this study are composed of a serial chain of n links having no closed-loops (open-chain configuration) and connected by single degree of freedom (DOF) rotational joints. The proximal link of the chain Chapter 2. Fundamentals of Robot Manipulator Simulation 6 Figure 2.1: Notation for link k. is connected to a fixed base element. The distal link of the chain is called the end-effector link. The end of the distal link is called the robot manipulator tip. The links are numbered 1 through n from the base to the tip. Each link k is considered to be a single rigid body having a uniformly distributed mass rrifc. The link center of mass is located at the point C\ and Iq is the link's inertia matrix about CV To the proximal end of each link is attached a fixed coordinate frame. Link coordinate frames are associated to the robot manipulator according to the Denavit-Hartenberg convention [15] (see figure 2.1). The vector from the location of the origin Ok of link coordinate frame k and the origin Ok+\ of the successive link coordinate frame is lk,k+i- The vector from Ok to the link center of mass Ck is pk. Xk+i,k is defined as the coordinate tranformation matrix from frame Ok+i to frame Ok 1-Each joint k is characterized by a unit joint axis vector, for a revolute joint, Zk — [ 0 0 1 ] T in the direction of the z axis of the links coordinate frame. Joints are transformation operators will be denned later in this chapter. Chapter 2. Fundamentals of Robot Manipulator Simulation 7 equipped with actuators which apply torques 77. about the joint axis based upon the robot's particular control law.2 The masses of the joints and actuators are considered to be incorporated into the description of the link. Kinematic joint limits and dynamic actuator limits will not be considered. The manipulator is assumed to be in a gravity field characterized by a constant linear acceleration vector g. At first, the effects of joint damping friction3 and external contact forces will not be considered. However, the topic of incorporating contact simulation into a robot simulation will be briefly introduced in a later chapter of this thesis. Overall, the type of multibody serial chain considered in this thesis allows one to build basic models of many of today's popular robot manipulator designs. Serial multibody chains also form the basic components of many of the more complex multibody robotic system designs. Hence the importance of developing the capabilities to efficiently simulate these basic types of robot manipulators. 4 2.3 The formulation problem A preliminary step in solving the simulation problem for an n-link robot manipulator is to formulate its equations of motion. These equations relate forces and torques to accelerations given the current state of the robot and its inertial properties. This step is often referred to as solving the formulation problem. There are two approaches to formulating the equations of motion of multibody sys-tems: elimination methods and augmentation methods [52]. Elimination based methods (or formulations in configuration space) are based on defining the equations of motion in 2Issues related to robot control will not be discussed within the context of this thesis. Many of our analyses will consider the robot manipulator to be an unactuated multibody chain controlled by a gravity law. 3For more information on the impact of friction on robot simulation see [17] 4The simplicity of the proposed robot model will allow us to better concentrate on specific issues related to robot simulation. Chapter 2. Fundamentals of Robot Manipulator Simulation 8 terms of generalized coordinates (e.g. joint angles) which span the constraint manifold. For open-chain problems, these methods result in a compact set of ordinary differential equations (ODE). However, these methods are not well designed to incorporate con-straints that are not integrable (non-holonomic constraints) or discontinuous constraints (e.g. contacts with the environment). Augmentation based methods (or formulations in descriptor space) start from a for-mulation of the dynamics in Cartesian space and restrict the motion of the bodies with constraint forces (Lagrange multipliers) applied at each joint. These methods result in a relatively large set of differential algebraic equations (DAE) (differential equations of motion with algebraic constraint equations). These equations are generally more difficult to integrate than ODEs and can be viewed as very stiff ODEs. One of the advantages of these methods is their capability of handling both holonomic and non-holonomic con-straints and allowing the possibility of adding and removing such constraints on-line (without a complete reformulation of the equations). The main focus of this thesis is the study of algorithms for simulating the motion of robot manipulators in open-chain configurations, for example robots without closed kinematic loops and robots whose motion is not constrained by contacts with the envi-ronment 5 . For the simulation of robots in open-chain configurations, elimination based methods are usually preferred since they lead to compact sets of equations which can be directly integrated. There are two popular elimination based methods for formulating the equations of motion in configuration space: Lagrangian methods and Newton-Euler methods. Lagrangian methods are based on formulating the equations of motion of the robot using an analysis of the total energy of the system. Newton-Euler methods are based on a detailed analysis of the forces at each link. Both of these methods have been 5Chapter 6 will introduce issues related to the simulation of contacts between a robot and its environment. Chapter 2. Fundamentals of Robot Manipulator Simulation 9 shown to lead to equivalent sets of equations [54]. They lead to equations having the following configuration space form: T = K{9)9 +b(9,9) where: • r is the n x 1 vector of torques (forces) applied by the joint actuators, • 9 is the n x 1 vector of joint variables (9 and 9 are the joint velocities and accelerations), • H is the n x n joint-space inertia matrix (JSIM) or generalized inertia matrix, • b is the n x 1 bias vector representing the torques (forces) due to gravity, centrifugal and Coriolis accelerations, and any external moments and forces acting on the end effector of the manipulator. For complex manipulators having a large number of degrees of freedom, manually writing out the equations of motion can be a very complicated and time consuming task. Many researchers have developed algorithms that automate the generation of the equations of motion (e.g. [44]) and several commercial software packages are available today (e.g. SD/FAST, ADAMS, DADS). Once the equations of motion are formulated, they can either be solved for joint torques (the inverse dynamics problem) or for joint accelerations (the forward dynamics problem). Inverse dynamics can be easily solved in 0(n) time as will be demonstrated later in this chapter. Forward dynamics is a more difficult problem and will be discussed in detail in the following chapter. Chapter 2. Fundamentals of Robot Manipulator Simulation 10 2.4 Dynamics of a single rigid body In this section, the equations of motion of a single rigid body are developed. These equations form the basis of the following sections which develop the equations of motion for a serial chain of rigid bodies. This section also serves to introduce the spatial notation which will be used throughout this study. 2.4.1 Spatial vectors A coordinate free spatial notation [20, 32] will be used to describe the motion of rigid bodies and the applied forces. In general, spatial vector quantities can be viewed as 6 x 1 column vectors which group the linear (translational) and angular (rotational) components of motion and force, while a spatial operator is a linear transformation of spatial vectors, with a 6 x 6 coordinate matrix. The spatial velocity VQ and the spatial acceleration do of a point O on a rigid body B are defined as vo = vo wo do * cio — vo Where u<n and ao are the 3 x 1 dimensional linear velocity and acceleration of the point O on the rigid body, and u>o and do are the 3 x 1 dimensional angular velocity and acceleration of the rigid body. Unless otherwise specified, all velocities and accelerations are assumed to be relative to a global fixed inertial frame. Similarly, the net spatial force / 0 acting on a rigid body B at a point O is defined as So TO Where f 0 is the 3 x 1 dimensional force (linear) acting on the rigid body B and TO is the 3 x 1 dimensional moment (torque) acting about O on the rigid body B. Chapter 2. Fundamentals of Robot Manipulator Simulation 11 These definitions of spatial velocity, acceleration and force differ slightly from the definitions given by Featherstone. In essence, our definitions specify two different types of spatial vectors: spatial motion vectors (a combination of a free vector with a line vector) and spatial force vectors (a combination of a line vector with a free vector). Mappings between two spatial vectors of the same type are discussed in the following subsection. Mappings from/to spatial motion vectors to/from spatial force vectors are commonly referred to as equations of motion. 2.4.2 Rigid body transformation operator The rigid body transformation $o,c 1 S defined as a spatial operator relating spatial ve-locities/forces at a point O on a rigid body to spatial velocities/forces at a second point C on the body 6 . This transformation is a function of the vector lo,c connecting the two points . r -1 'o,c 0 1 Important properties of this operator include 1 0 (lo,e)T 1 o,c CO, *o,c = $o,v$v,c-where V is a third point on the rigid body. Using this operator, the relationship between spatial forces at C and O is Jo = *c,ofc-Which corresponds to the following well known relationships between forces and torques on a rigid body Jo = Jc> f0 = fc + To,cJc-6This transformation should not be confused with the coordinate transformation operator which will be discussed in the following subsection. 7The coordinates of lo,c are assumed to be expressed with respect to frame O. Chapter 2. Fundamentals of Robot Manipulator Simulation 12 Similarly, $o,c can be used to express the relationship between the spatial velocities at C and O Which corresponds to the following well known relationships between linear and angular velocities vo = vc + &clc,o, &o = The transformation of spatial accelerations is a little more complicated. It is given by do = ^ofio-c + Pc> where pc is the velocity dependent Coriolis and centrifugal spatial acceleration term at the center of mass Pc = $o,cvc = 2.4.3 Coordinate transformation operator The above notation is said to be coordinate free since the quantities are not bound or referenced to a specific coordinate frame. This type of notation allows us to obtain a more compact expression of the equations of motion. However, in order to solve the equations of motion, the components of these quantities must be compared and computed with respect to common frames of reference. Therefore, coordinate frame transformation operators (rotation matrices) must be applied to the quantities considered. These operators are defined as follows: X o,c = $ofiRo,c, where Ro,c = 0 Eo,c 0 0 Eo,c Chapter 2. Fundamentals of Robot Manipulator Simulation 13 where Eo,c 1S the standard orthogonal rotation transformation matrix between the coor-dinate frame O and the coordinate frame C [15]. 2.4.4 Newton and Euler's equations of motion If C is the location of the center of mass on the rigid body, then the spatial inertia of the body about this point is defined as Mc = where m is the total mass of the body and Ic is the body's moment of inertia about its center of mass. The relationship between the net force applied at the center of mass and the resulting linear acceleration of the body is given by Newton's equation: f c — d(mvc)/dt = mac-The relationship between the net torque and angular acceleration is given by Euler's equation: fc = d(Ic&c)/dt = Ic&c + u>clc&c-These equations result in the following spatial relationship fc = Mcac + ic, where 7 C is the velocity dependent spatial gyroscopic force acting at the center of mass: 7C = Mcvc. Using the previously denned spatial rigid body transformation, the relationship be-tween the spatial forces and accelerations can be rewritten for any point O on the rigid body as follows: fo = M 0 a 0 + loi ml 0 0 Ic Chapter 2. Fundamentals of Robot Manipulator Simulation 14 where and M 0 = X c 0M CX c,o ~ t lo = Xc,o(McPo + 7 c ) = 2.5 The inverse dynamics problem uolo^o mQoCJolojc In order to demonstrate the extension of our notation to multibody chains, this section presents an 0(n) recursive algorithm [39] for solving the inverse dynamics problem. As will be shown, inverse dynamics can play a direct role in solving the forward dynamics problem. This section also serves to introduce and define certain generalized spatial operators. The notation for developing the dynamics of multibody chains is a natural extension of the notation used in the analysis of the dynamics of a single rigid body. The following is a summary of the notation associated with multibody chains: Link properties: • /j,t+i is the vector from the origin of link t to the origin of link t -f-1, • pi is the vector from the origin of link i to the center of mass of link i, • Mi is the spatial inertia matrix of link i, • Xi+iti is the transformation matrix from fink i + 1 to link i, Joint properties: • d ,9 ,9 are the generalized relative joint angles, rates and accelerations, where 9i,9i,9i are the position, velocity and acceleration of joint i. Chapter 2. Fundamentals of Robot Manipulator Simulation 15 • T is the generalized joint torque, where r$ is the torque applied to link i about its joint axis. • si is the unit axis of the joint between link i — 1 and link i (for a revolute joint Si=[0zi}T), Generalized velocities and accelerations: • i is the generalized spatial velocity, where t>( is the spatial velocity of the origin of link i. • a is the generalized spatial acceleration, where is the spatial acceleration of the origin of link i. • p is the generalized spatial Coriolis and centrifugal acceleration, where p t is the spatial Coriolis and centrifugal acceleration of the origin of link i. Generalized forces: • / is the generalized spatial inter-link force, where f { is the spatial force exerted on link i by link i — 1 at joint i. • 2 is the generalized spatial gyroscopic force, where 7^  is the spatial gyroscopic force acting on link i. • fT is the generalized spatial net force, where fTi is the net spatial force exerted on link i. The inverse dynamics algorithm is based on work done by Luh, Walker and Paul [39]. The algorithm proceeds in two recursive iterations: an outward iteration from base to tip to compute the spatial velocities and accelerations of each link in the chain, and an inward iteration from tip to base to compute the link forces and joint torques. The Chapter 2. Fundamentals of Robot Manipulator Simulation 16 proposed notation leads to a compact formulation of the familiar recursive Newton-Euler formulaton of inverse dynamics. The outward recursive relationships between the spatial velocity/acceleration of link i - h l and the spatial velocity/acceleration of link i are di+i — Xi+ijcii + S j + i # i + i + p i + 1 , where r ~ - = = -for i = 0 to n — 1 with VQ = 0 and do = 0 (note: the effect of gravity can be included by setting do = <?)• The inward recursive relationship to compute inter link spatial forces from the tip to the base is fi = Midi + 7^  + xf+lJi+i, where the joint torques are given by the projection Pi+i = Xi+ijVi + § i + i 9 i + i = Ti = Sjfit for t = n to 1 with / n + 1 = 0. 2.5.1 Generalized spatial operators Rodriguez, Jain and Kreutz-Delgado [50] introduced a powerful generalized spatial oper-ator algebra which allows the equations of motion of multibody chains to be expressed in a much more compact high-level form. Generalized spatial operators are n x n gener-alized matrices whose elements are spatial entities.8 The following is a list of some key Generalized spatial operator are similar to the "super matrices" defined by Ellis [18] Chapter 2. Fundamentals of Robot Manipulator Simulation 17 generalized spatial operators: M_ = diag Mi M2 .. • Mn S — diag S2 ••• 1 0 0 . . . 0 $2,1 1 0 . . . 0 i = $3,1 $3,2 1 . . . 0 _ * n , l $n,2 $n,3 1 1 0 0 0 -^2,1 1 0 0 X. ^3,1 ^3,2 1 0 Xn,2 XN,3 1 Generalized spatial operators can be interpreted as encapsulating recursive operations. Three structures are of particular interest: block lower triangular operators represent recursions from base to tip, block upper triangular operators represent recursions from tip to base and block diagonal operators are non recursive or memoryless operations. For example, y = X_ x is equivalent to the recursion: y i + l — Xi+i-\-Xi+nyi, and y = X_ x is equivalent to the recursion: yi = z$ +. X i + X i y i + l . Generalized spatial operators may also be referred to as composite spatial operators since they act on the whole chain at a given configuration. For example, $ is often referred to as a composite rigid body transformation. These generalized spatial operators can be used to express the previously described inverse dynamics algorithm: v = X S 9 Chapter 2. Fundamentals of Robot Manipulator Simulation 18 o =X_(S.9 +p) / = i r ( i « +7 ) ~ T ~ which can be rewritten to obtain the familiar state space expression for the equations of motion of robot manipulators: r = H 9 + b where H = STX TM X S b =STX_ T(M X. p + 7 ) 2.6 The motion integration problem Although the motion integration problem is not the focus of this study, we introduce here some of the basic issues related to this part of the simulation process. In robot simulation, the motion integration problem is defined as follows: given the current state of the robot (joint positions and velocities) and the current joint acceler-ations at time t compute the state of the robot at a time t + h, where h is called the integration step size. For numerical simulation, the second order differential equations of motion are converted to a coupled first order differential equation: d dt ' 9(t)^ K9(t) , ( 9 (t) ^ = f(9,9,t) \ 9 ( t ) J Numerical integration is not an exact procedure. It introduces truncation errors which are a function of the size of the time step and the smoothness of the derivative function. It is therefore desirable for an integrator to exert control over the step size; adjusting the Chapter 2. Fundamentals of Robot Manipulator Simulation 19 time step so as to maintain a desired level of accuracy. Such integration schemes are said to be adaptive step size methods. In this thesis, we are principally interested in the interaction of the forward dynamics and the motion integration with respect to adaptive step size integration. In the context of an adaptive step size integrator, a forward dynamics function will be said to exhibit operational stiffness if there is a sharp increase in the number of derivative function evaluations necessary to reach a stated accuracy [30]. These issues will be discussed in chapter 4 and illustrated in chapter 5. For a more general introduction to numerical integration see [25] and for further discussions of the relationship between motion integration and robot simulation see [41, 2]. Chapter 3 Solving the Forward Dynamics Problem Existing algorithms for solving the forward dynamics problem for open chain robot ma-nipulators can be classified into three categories based on their computational complexity: 1. 0(n 3) algorithms based on inertia matrix inversion, 2. 0(n2) conjugate gradient methods, 3. 0(n) linear recursive methods. A good overview of the relationship between the different types of forward dynamics algorithms can be found in [32]. Jain uses a spatial operator algebra to analyze and compare the different classes of forward dynamics algorithms from a unified perspective.1 3.1 0(n 3) Algorithms: inertia matrix inversion methods The 0(n 3) methods for solving the forward dynamics problem are said to be based on inertia matrix inversion. The starting point of these methods is the following state space linear equation: Kl =L -b (3.1) These methods follow three basic steps. 1. Calculate the bias vector (b ). This generalized vector represents the action of forces due to gravity, centrifugal and Coriolis accelerations, and any external forces 1 T h e algorithms reviewed in this chapter are for sequential computation. A review of parallel algo-rithms for robot forward dynamics can be found in [22]. 20 Chapter 3. Solving the Forward Dynamics Problem 21 acting on the end effector of the manipulator (e.g.: contact forces). The bias vector is a function of the current joint configuration 9 and joint velocities 9 . It can be calculated by solving the inverse dynamics problem for 9 =0 (b are joint torques necessary to cause zero acceleration). For this purpose, an efficient 0(n) inverse dynamics algorithm such as the one presented in the previous chapter can be used [39]. 2. Calculate the joint space inertia matrix H_ (or generalized inertia matrix). It is in the method by which this is performed that most 0(n 3) algorithms differ from one another. One of the most efficient methods for computing H_ is known as the Composite Rigid Body Method (CRBM) [58]. This method has a computational complexity of 0(n2). It will be described below. 3. Solve the linear equation for the joint accelerations (inertia matrix inversion step). This step is usually accomplished by such techniques as LU decomposition or Cholesky decomposition [28, 56]. The fact that the joint space inertia matrix is symmetric positive definite means that solving this step should not require any pivoting. The computational complexity of this step is 0(n 3). Although inertia matrix inversion methods have an overall computational complexity of 0(n 3), the coefficient of n 3 is very small (e.g. 1/3) making these methods efficient for small values of n. 3.1.1 The composite rigid-body method Walker and Orin [58] present three methods for computing the joint space inertia matrix (step 2 of the above procedure). Their third method, known as the composite rigid-body method (CRBM), is the most efficient of the three. In this method, a composite rigid-body j is defined as a single rigid body formed by the chain of links from j to n, where Chapter 3. Solving the Forward Dynamics Problem 22 joints j ; + 1 to n are considered fixed. The composite rigid body n is the distal link of the chain. The C R B M proceeds in three steps. 1. For j'< = n to 1, consider the composite rigid-body j. Compute the position of its center of mass, its total mass and inertia matrix with respect to joint j (using the parallel axis theorem). These calculations depend only on the current geometric configuration of the manipulator. 2. For j = 1 to n, consider joint j as the only joint in motion (set 9j = 1 and 9i = 0, for all i ^ j), compute the total forces and moments applied to the center of mass of each composite rigid-body j. 3. For j = n to 1, for i = j to 1, calculate the forces and moments applied to link i — 1 by the composite rigid-body i. If joint i is rotational, then the (i, j)th component of H is equal to the z component of the computed moment. This step takes advantage of the symmetry of H_ . Using the generalized spatial notation, the joint space inertia matrix can be denned as follows [32]: K = lTKTKX_S . (3.2) 3.1.2 Orthogonal complement based methods Two 0(n3) methods for solving the forward dynamics were developed by Angeles and Ma [3]. In these methods, the concept of a natural orthogonal complement f_ (a6nxn matrix of unit twists) is used to construct the generalized inertia matrix H_ of the manipulator. The use of a natural orthogonal complement describing the constraints of the manipulator allows the authors to obtain a factorization of the generalized inertia matrix into H_ = Chapter 3. Solving the Forward Dynamics Problem 23 T -P_ P_ , where P_ is obtained through a factorization of the generalized extended mass matrix. Their first method computes the generalized inertia matrix explicitely by calculating the product H_ = P_ P_ and uses Cholesky decomposition to solve the linear equations of motion for joint accelerations. Their second method does not explicitely form H_ . Instead, the following system of equations is solved: P is a 6N x N block lower triangular matrix; in brief, this method solves an underdeter-mined system followed by an overdetermined system (using the minimum-norm solution and the least squares approximation). They observe that this second method, although less efficient than their first has better numerical properties and is not subject to as potentially high condition numbers as other 0(n 3) methods. The problems related to ill-conditioning of the generalized inertia matrix will be further discussed in the following chapter. 3.2 0{n2) Algorithms: conjugate gradient methods Conjugate gradient methods provide a general 0(n2) means for solving a set of linear equations of the form: Ax = b. Conjugate gradient methods are guaranteed to converge to the solution in, at most, n iterations. In each iteration the matrix (A) is accessed only through multiplication by a vector (y). The product Ay can be obtained in 0(n) time using inverse dynamics. Walker and Orin (method 4 [58]) describe an 0(n2) algorithm for solving the forward dynamics problem based on a conjugate gradient method. The procedure follows the pTy = T-b (3-3) Pq = y (3-4) Chapter 3. Solving the Forward Dynamics Problem 24 same basic steps as described in previous section except that in the last step, a conjugate gradient method is used to solve the linear equation. A similar 0(n 2 ) procedure is also described by Fijany and Scheid [23]. 3.3 0(n) Algorithms: linear recursive methods In general, 0(n) methods are based on propagating the coefficients of motion and force along the robot manipulator in an inward (tip to base) and outward (base to tip) recursive procedure. These methods can be viewed as forming an equation of the following form: I =F(9,r -b) (3.5) where F (the forward dynamics function) is expressed without reference to the joint space inertia matrix. One of the first such procedures was developed by Vereshchagin [57]. Later Armstrong [4] also demonstrated the existence of a linear recursive relationship between the motion of and forces applied to link i, and the motion of and forces applied to links i — 1 and i + 1 of a robot manipulator. This recursive relationship allows the joint accelerations to be determined in 0(n) time. However, the coefficient of n is rather large. One of the most familiar linear recursive methods is Featherstone's articulated-body method (ABM) [19, 20]. Featherstone obtains a reduction in the number and size of equations by the introduction of a spatial notation for accelerations, inertias and forces which combines the translational and rotational components of motion. Featherstone also defines the concept of articulated-body inertias: "the relationship between a spatial force applied to a particular member of an articulated body and the spatial acceleration of that member, taking into account the effect of the rest of the articulated-body". In other terms, the spatial articulated-body inertia A; is: the effective spatial inertia at the Chapter 3. Solving the Forward Dynamics Problem 25 kth. joint of links k to n taken as a single non-actuated articulated-body. The articulated body method will be presented in detail later in this section. Another 0(n) method is presented by Brandl, Johanni and Otter [11]. This method is very similar to Featherstone's algorithm, but uses a more general joint model to allow for the modeling of multiple degree of freedom joints. 2 It also obtains a better performance than Featherstone's by optimizing vector and matrix coordinate transformations. Finally, Rodriguez, Jain and Kreutz-Delgado [49, 50] also present a recursive O(n) method. This method is based on developing an analogy between robot dynamics and techniques from linear filtering and smoothing theory (Kalman filtering and Bryson-Frazier smoothing). The method uses a specially developed spatial algebra based on a spatial notation similar to the one introduced by Featherstone. Although the resulting computations are very similar to the ones of Featherstone, the means by which the equations are derived are very different. This method does not require an initial step to calculate the bias terms (inverse dynamics step), the bias term is naturally incorporated in the formulation in such a way that only a single inward filtering and outward smoothing iteration is necessary to solve the forward dynamics problem. Featherstone in [20] also developed an AB algorithm that did not require pre-computation of the bias term. 3.3.1 The articulated-body method The articulated-body method proceeds in three steps. 1. Calculate the bias vector (b ). This is accomplished using an inverse dynamics procedure (as was done in step 1 of the 0(n3) procedure). 2 For multibody dynamics formulations supporting multiple degree of freedom joints see also [48, 36]. Chapter 3. Solving the Forward Dynamics Problem 26 2. Inward recursion. Moving from the tip to the base, calculate the inertias of each articulated-body / , and the associated correction forces pi (these forces incorpo-rate the effects of non-zero forces in the joints of the articulated-bodies). These quantities enable the calcidation of the spatial forces applied at each joint. 1 r £ _ V . , .T(fi • ~ fr+l) ~ S i + 1 p i + l ) Pi- * t+i,. [Pi+i H — T ^ ),p n = 0 Si+l / i + 1 S i + i (3.7) 3. Outward recursion. Calculate the link accelerations a* and joint accelerations 9{ by a recursive propagation from the base to the tip of the robot. O j = Xit-iiii-i + SiOi, oo = 0 (3.8) §i = ( t > ' b i ) ~ fo^-A-i - SJPJ (3 9) SiTIi Si Although this algorithm has a computational complexity of 0(n) it is shown to be less efficient than Walker and Orin's composite-rigid body method for n < 9 [20]. Physical interpretation: force and acceleration propagators The inward recursive iterations of the articulated body method serve to encode the inertia! properties of the arm and to propagate bias forces from the tip to the base. The force propagator PFi can be defined as P* = 1 - (3-10) Si I{ Si Using the force propagator, the tip to base recursion of the articulated-body method can be rewritten as follows I? = u + x i + 1 / p F i + J i + l x l + h i , (3.11) Chapter 3. Solving the Forward Dynamics Problem 27 and • " A -Pi = A ( P F i + 1 p i + 1 + — t . a ). (3.12) The outward recursive iteration of the articulated body method propagates accelera-tions from the base to the tip. The unitless acceleration propagator operator P^ can be defined as 3 PA, = (1 - ^ - ) . (3.13) Using the acceleration propagator, the tip to base recursion can be rewritten as follows o-i = P.AiXi,i-iOi_i -I . (3.14) S i S i It is interesting to note that this operator is in fact a weighted least squares projection operator. In other words, the acceleration propagator could be interpreted as an operator that projects motion vectors along the direction that moves the least mass. This would be in agreement with Gauss' principle of least constraint [34, 57]. 3.4 Summary and comparisons In an attempt to better compare the various forward dynamics algorithms, Jain [32] used the generalized spatial notation, to obtain a closed-form expression for the JSIM based on generalized spatial operators. Jain showed that the interpretation of the generalized operators in this expression leads directly to the A B M algorithm. In [5], using a unified formulation of multibody dynamics and basic linear algebra, it is shown that the C R B M and A B M are in fact simply different elimination methods for solving the same system of algebraic equations of motion. At a first glance, a good way to compare the performance of the C R B M and the A B M is to look at the number of operations required. The following table summarizes 3 This type of spatial transformation can be called a spatial articulated transformation as opposed to a composite-rigid body transformation [32]. Chapter 3. Solving the Forward Dynamics Problem 28 the operation counts for the two methods. These results were obtained by Featherstone [20] and account for the fact that each method has been optimized so as to remove redundant calculations. Method Number of multiplications Number of additions A B M 300n - 267 279n - 259 C R B M (Cholesky) n 3 /6 + 11.5n2+ (160 + l /3 )n-109 n 3 /6 + 7n 2+ (138 + 5/6)n - 102 Based on this table, in terms of the number of operations, the A B M is more efficient than the C R B M for n > 9. However, as will be demonstrated in the following chapters, the performance of a simulation program does not only depend on the efficiency of the forward dynamics algorithm but also on its numerical performance. Most of the literature on forward dynamics does not discuss stability and other numerical issues related to the implementation of these methods [18, 55]. Chapter 4 Numerical Properties of Forward Dynamics Algorithms In this chapter, the focus will be placed on the analysis of the numerical properties of two of the most popular methods for performing forward dynamics: • the Composite Rigid-Body Method (CRBM) • the Articulated-Body Method (ABM) In addition to examining the numerical properties of these forward dynamics methods, we will examine the effects on the whole simulation process (i.e. including motion inte-gration) . 4.1 Review of issues in numerical computation Roundoff errors Computers store floating-point numbers with finite precision. The machine accuracy em is defined as the smallest floating-point number that when added to 1.0 gives a result different from 1.0. A computer operation between two floating-point numbers adds a fractional error of at least em to the results. This error is called the roundoff error. For example, subtraction of two nearly equal numbers can result in relatively large roundoff errors. A good numerical algorithm should attempt to avoid such inaccuracies. One of the important issues in numerical computation is determining how roundoff errors accumulate in successive/iterative computations. An algorithm will be numerically unstable if roundoff error is propagated in such a way that the errors are magnified. 29 Chapter 4. Numerical Properties of Forward Dynamics Algorithms 30 Some of the current methods for estimating the numerical stability of algorithms include interval analysis and statistical analysis. A very good reference on the topic of roundoff errors in numerical computation is [60]. Matrix computations The condition number n(A) = ||J4||||./4. - 1|| is a popular way of measuring the sensitivity of a linear system Ax — b. The two-norm condition number of a matrix corresponds to the ratio of the largest singular value to the smallest singular value [28]. A large condition number indicates a poorly conditioned system. The condition number is only a worst case estimate of the potential sensitivity to errors, it can only provide a bound on the relative errors. If K(A) is the condition number of A and 6b is a perturbation in b then Sx, the perturbation in x, will be bound by: Truncation errors Approximate numerical algorithms result in truncation error. In the case of robot sim-ulation, truncation error is introduced only in the numerical integration process since forward dynamics is an "exact" algorithm. Adaptive step size integration schemes pro-vide a good means by which truncation error can be controlled. 4.2 Numerical properties of the C R B M As previously discussed, an 0{nz) forward dynamics algorithm typically calculates the joint space inertia matrix (JSIM) completely, and then solves an n-dimensional linear system. Forming the JSIM can contribute to a loss of precision due to round-off errors (round-off errors occur when relatively "large" quantities are added to smaller ones during Chapter 4. Numerical Properties of Forward Dynamics Algorithms 31 the computation of the JSIM). Using different linear solution algorithms does not help since once the JSIM is formed, the precision is lost and cannot be easily recovered. Furthermore, the JSIM is symmetric positive definite (SPD) and therefore should not require any pivoting to better solve the linear system unless roundoff errors are such that the matrix is no longer SPD. However, the JSIM could be poorly conditioned. In robot simulation, ill-conditioning of the JSIM is more common than generally appreciated. Popular robot designs have intersecting joint axes at the base of the manipulator and at its wrist, leading to one or more links of small mass and link length. In addition, there are typically large dis-crepencies in link length and mass between links of the regional structure and those of the orientation structure. Some ill-conditioning may be caused by the process of model-ing, when the masses of some links are neglected or poorly approximated. Finally, the condition number is also a function of the configuration of the manipulator. Ill-conditioning of the JSIM is an important indication of difficulty. This follows from the well known properties of the sensitivity of linear systems solutions [28]. Suppose the final step of the C R B M is subjected to 0(e) perturbations {H + eG )q (e) = r -b + eg , (4.2) where G € &NxN, g 6 XN and K{H_ ) = ||/f || \\H_ is the condition number, and the relative error in quantity x is written as ,(*) = | | x(6 ) - s | | / | | * | | (4.3) then the relative error in g is bounded above as follows p(q)<K(K)[p(H) + p(T -b)]+0(e2) (4.4) We observe that this bound is met in practice, as can be seen in the experimental results of chapter 5. Chapter 4. Numerical Properties of Forward Dynamics Algorithms 32 4.3 Numerical properties of the A B M We now show that the 0(n) articulated-body method naturally possesses good numerical properties for ill-conditioned problems. A precise analysis of roundoff error propagation in this case is difficult to perform but the results of this section and the numerical experiments of the following chapter clearly demonstrate the phenomenon. In this section we will consider an example which will illustrate the numerical dif-ferences between the forward dynamics methods and will provide an explanation for the ABM's observed robustness to ill-conditioning. A more detailed analysis of this phe-nomenon in terms of underlying matrix factorizations can be found in [5]. In related work, [18] showed that the A B algorithm is equivalent to performing Gaussian elimina-tion on a symmetric positive definite matrix, and therefore stable as long as the positive definiteness is maintained. However, they did not observe any numerical superiority of the A B M over other methods such as C R B M . Consider the example of a two link manipulator. To compute Q\ using the A B M (see previous chapter), we solve where c 2 is a computed bias force and lx is the articulated-body inertia of link 1 (4.5) C2 = (TI - bi) - sf X 2 , i M2s2(s2I2 s2) \T2 - b2) (4.6) 11 = Mi+X2,\ (M2 - M2s2s2TM2(s2TM2s2)~l)X2ti. (4.7) Writing out the expression for 9\ gives: (n - bi) - sY X 2 , i Af 2s2(s212 s2) (*2 - h) (4.8) a f ( M i + X2,i ( M 2 - M2s2s2TM2(s2TM2s2)-l)X2A)si Chapter 4. Numerical Properties of Forward Dynamics Algorithms 33 Now consider the C R B M . For the two link robot manipulator, the joint space inertia matrix can be written as (see, e.g. [36]) s r ( M 1 + i 2 , 1 T M 2 i 2 . i ) 5 i SJX2/M2S2 ^ & s?M2X2jiSi s2M2s2 J (4.9) To solve for q\ using the 2 x 2 matrix H_ , we eliminate the upper right element by a block-row operation (multiply the second row by sjx2,i M 2s2{sli M 2s2)~l and subtract it from the first row). This gives the following expression for 0\. n (jl ~ 6l) - s\X2^ M2S2{slJA S2)~KT2 - 62) , . "I,CRBM — ; ; — — 1 ;—~T~~2 1 1 — I • v*-*^/ i [ ( M i + i 2 , i M2X2,x)si - sjx2,i M2s2(sTM2s2)-ls2M2X2,ih This expression is of course equivalent to equation 4.8. The numerator is computed in exactly the same order in both cases. The difference between the A B M and C R B M lies in the order in which the denominator (s, Ix si) is computed. Suppose the distal inertia (M 2 ) is much larger than the proximal inertia (Mi) (ill-conditioned chain). In the C R B M we must, while forming H_ , add the small M\ to the large X 2 , i M2X2ji. Then the other large term s 1 X 2 , i M 2 s 2 ( s 2 M 2 s 2 j s2 M2X2,iSi is subtracted during the solution process. This may cause cancellation error which appears as a random function of time t, even though all the quantities are smooth functions of t. This error does not appear in the A B M algorithm. However, it is important to take note that the A B M can be implemented so as to perform poorly in the presence of ill-conditioning. For example, if the articulated body inertia is computed in the following manner: -A , ~ T - - - TM2S2S2TM2 -Ix =(Af 1 + A- 2 ,i M2X2A)-X2A — ^ s2 M2s2 then cancellation error may also appear. Chapter 4. Numerical Properties of Forward Dynamics Algorithms 34 4.4 Impact on motion integration: formulation stiffness The numerical simulation problem is usually treated as two separate problems: the for-ward dynamics problem for computing system accelerations, and the numerical inte-gration problem for advancing the state in time. However, when the noise caused by cancellation error in the C R B M becomes relatively large, the forward dynamics presents a nonsmooth solution profile. This profile may slow down an unsuspecting adaptive step size motion integration routine. In such cases, the interaction between the forward dy-namics and the numerical integration is an important factor that must be considered when evaluating the overall efficiency of the robot simulation process. Define the formulation stiffness in multibody simulation as any operational stiffness encountered by the motion integration routine that is due to a poor numerical formulation of the simulation algorithm (i.e. any stiffness that can be avoided by a reformulation of the solution algorithm). Based on this definition, we may state the following: the C R B M will exhibit formulation stiffness when used in the simulation of ill-conditioned multibody chains. The A B M should perform better then the C R B M in ill-conditioned situations. This will be experimentally verified in the following chapter. Chapter 5 Implementation and Experimental Verification In order to observe the numerical properties discussed in the previous chapter, three computer programs were developed: a program for testing only the forward dynamics component, an off-fine robot simulation program, and an animation program used to vi-sualize the simulation results. The first part of this chapter describes the implementation of these programs and the following sections present some of the results obtained with this program. In section 5.2, we first focus on the numerical performance of the forward dynamics component before broadening the analysis to incorporate the whole simulation process (forward dynamics and motion integration). Section 5.3 presents some other interesting operational differences between the C R B M and A B M algorithms. 5.1 P R R R o b o t simulation software Robot simulation software was developed for comparing the properties of the C R B M and A B M algorithms. The software consists of a specialized program to test and compare forward dynamics (called prrrfdt), an off-line robot simulation program (called prrrsim) and an animation program (called prrranim) used to visualize the simulation results. Both the A B M and C R B M algorithms were implemented for the basic case of a general n-link planar, revolute jointed, rectangular linked robot. 5.1.1 Assumptions and definitions In general, a robot simulation program manages the following types of data: 35 Chapter 5. Implementation and Experimental Verification 36 • Robot parameters: characteristics of the robot that remain constant throughout the simulation (e.g. link lengths, weights, widths, joint friction coefficients and link moments of inertias). • Robot variables (state): characteristics of the robot that may vary at each time step (e.g. joint positions, velocities and accelerations). • Environment parameters: characteristics of the environment that remain constant throughout the simulation (e.g. gravity, environment friction coefficients). • Environment variables (state): characteristics of the robot that may vary at each time step (for simulation in dynamic environments, e.g. positions of moving obsta-cles). 1 • Simulation parameters: e.g. numerical integration accuracy, minimal step size. • Simulation statistics and results: e.g. history of joint positions, number of forward dynamics function evaluations, execution time, energy measures. The implementations consider a special class of robot manipulators. These robots have the following attributes: • robot joints are single axis revolute joints and have a full 360 degree freedom of motion (inter-link collisions are not considered) the joints are also considered to be free of friction, • motion of the robot motion is constrained to a two dimensional plane (i.e. all joint axis are parallel), lAt this stage, robot interactions with the environment are not considered. Therefore, environment variables are non-existent and environment parameters are greatly simplified. See chapter 6 for a dis-cussion of the issues involved in simulating environment interactions. Chapter 5. Implementation and Experimental Verification 37 LinkN A : > Link 2 Figure 5.1: N-link PRRRobot manipulator. • robot links are flat rectangles having a uniform mass distribution and variable lengths and widths. We will refer to this basic type of robot manipulator as a PRRRobot. Planar, Revolute jointed, Rectangular linked Robots (see figure 5.1). 5.1.2 Forward dynamics test program The PRRRobot forward dynamics test program systematically compares the different for-ward dynamics solutions for a two link PRRRobot by performing a static sweep through the possible joint configurations. For example, the program may be configured to com-pute the forward dynamics for discrete values of 92 € [0,27r] while maintaining all other input variables constant. For each configuration of the chain, the program computes the joint accelerations using the articulated-body method (ABM) and the composite rigid body method (CRBM). The relative difference between the two solutions is then Chapter 5. Implementation and Experimental Verification 38 computed as follows: Ar9i = \{9ITABM — &i,CRBM)/9i,ABM\ where: • ®i,ABM is the computed joint i acceleration using the A B M , • 9i,CRBM is the computed joint i acceleration using the C R B M , • A r # t is the relative difference between the joint i accelerations computed using the A B M and the C R B M . If we consider O^ABM to be the "correct" solution, then Ar9i represents the relative error in O^CRBM-In this chapter, the prrrfdt program will be used to illustrate the effect of ill-conditioning on the C R B M and A B M forward dynamics algorithms. 5.1.3 OfF-line simulation program The off-line simulation program (prrrsim) reads the robot description from an input file and writes the simulation results to an output animation file and to various statistics files. The user may choose from two distinct forward dynamics algorithms (ABM or CRBM), various linear solution methods (e.g. Cholesky or LU decomposition) and vari-ous numerical integration methods (e.g. Runge-Kutta with adaptive or fixed step size). Figure 5.2 illustrates some of the run-time options provided by the program. Performance measures The PRRRobot simulation program generates statistics which allow monitoring and com-parison of the numerical performance of each forward dynamics method over the whole simulation process. The statistics are generated at a user-specified frequency and are compiled into several output files. These statistics include: Chapter 5. Implementation and Experimental Verification PRRRobot Simulator Test i d e n t i f i e r (max. 3 characters): T01 Name of robot data f i l e : prrr3.dat Forward dynamics method (1=ABM 2=CRBM): 2 Linear s o l u t i o n method: (1=CD 2=LU 3=GJE 4=SVD): 2 Integration method (1=RKA 2=STIFF 3=STIFBS): 1 Star t and stop times ( s ) : 0 10 Use d e f a u l t i n t e g r a t i o n parameters (y/n)? n F i r s t i n t e g r a t i o n step siz e ( s ) : 0.01 Minimal i n t e g r a t i o n step size (s): 0.0 Integration accuracy: 0.0001 Save i n t e r v a l ( s ) : 0.01 Maximum number of steps saved: 1000 Use d e f a u l t i n i t i a l robot st a t e (y/n)? n I n i t i a l j o i n t p o s i t i o n s (rads) and v e l o c i t i e s (rads/s): q[l) = 1.5 0 q[2] = 0 0 q[3] = 0 0 Use other d e f a u l t parameters (y/n)? n Gravity : 9.8 C o e f f i c i e n t of j o i n t viscous f r i c t i o n : 0.0 J o i n t actuator torques: t a u [ l ] =1.0 tau[2] =1.0 tau[3] =1.0 Working ... Done. Execution s t a t i s t i c s : Execution time: 1.923077 Average energy: 1313.936411 58 su c c e s f u l steps 38 bad steps . 801 function evaluations 87 steps saved Figure 5.2: PRRRobot off-line simulator interface and options. Chapter 5. Implementation and Experimental Verification 40 • The evolution of joint positions 9 ((), velocities 9 (t) and accelerations 9 (t). • The evolution of the total energy E(t) (kinetic and potential) of the robot. This can be used to confirm the conservation of energy throughout a simulation. • The number of forward dynamics function calls needed to advance the solution by A i (a user-specified update frequency). This number is a good measure of the evolution of the computational effort required to simulate a particular motion. • The evolution of the condition number K ( / £ . (t)) of the joint space inertia matrix of the robot. The two-norm condition number of a matrix is equal to the ratio of its largest singular value to its smallest singular value. In this implementation, the singular values are found using the singular value decomposition Numerical Recipes function [46] (not a particularly efficient technique, but it does the trick). • The evolution of the relative difference between joint accelerations (A r #i(£)) calcu-lated using the CRBM vs. joint accelerations calculated using the A B M based on the current C R B M solution (see equation 5.1.2). 2 5.1.4 Animation program The animation program (prrranim) was developed to enable simple visualization of the simulation results. This graphical display program was implemented using V O G L , a portable graphics library based on SGFs G L and available from the University of Mel-bourne via ftp. 2 To compute the relative difference at time t, both the A B M and C R B M accelerations are computed using the joint positions and velocities of the C R B M solution at time t. Chapter 5. Implementation and Experimental Verification 41 5.1.5 Implementation issues This section presents some of the details related to the implementation of the forward dynamics and simulation programs. The validity of the implementations was confirmed by comparing the obtained results to those of a commercially available simulation package (SD/FAST 3 ) . 4 Computation of inverse dynamics As was mentioned in chapter 3, the calculation of the inverse dynamics is used to compute the bias vector in both the A B M and C R B M forward dynamics algorithms. An efficient inverse dynamics procedure for a general n link PRRRobot was implemented. The inverse dynamics procedure is based on a recursive Newton-Euler formulation of the dynamics developed by Luh, Walker and Paul [39]. For the purposes of this project, the equations have been simplified by taking advantage of the attributes of the special class of robots being considered (PRRRobots). In particular, the following simplifications were possible: • For 2D motion, the rotational quantities: inertia matrices, torques, angular veloci-ties and accelerations can be represented by scalars instead of matrices and vectors (these scalars implicitly represent the component of rotation about the z axis). In addition, these values are invariant to changes in coordinate systems since all coordinate systems have the same z axis. 5 • In 2D, forces, linear velocities and accelerations can be represented by 2D vectors instead of 3D vectors. 3 S D / F A S T is a trademark of Symbolic Dynamics, Inc., 561 Bush Street, Mountain View, C A 94041 USA. 4 A l l programs were implemented in C and are highly portable. Versions of these programs were developed for P C compatibles, for S U N type workstations and for SGI workstations. For more details on the software implementations, the reader is encouraged to examine the documented source code. 5Inertia matrices are transformed based on the parallel axis theorem. Chapter 5. Implementation and Experimental Verification 42 • In 2D, the gyroscopic torque (u x Iu>) is zero. • By taking advantage of the rectangular structure of the links and the uniform mass distribution assumption, the calculation of the positions of the link center of mass and the calculation of the link inertias is also greatly simplified. Composite rigid-body forward dynamics Walker and Orin's composite-rigid body algorithm for computing the joint-space iner-tia matrix was implemented. The implementation takes full advantage of the simple structure of the special class of planar robots being considered. In addition to the pos-sible simplifications which were stated in the previous section, the calculation of the composite-rigid-body inertias simplify to simple scalar operations. Once the JSIM is computed, and once the bias forces are calculated using the inverse dynamics algorithm, a linear system of equations must be solved in order to obtain the joint accelerations. For this step, the simulation program offers a choice of several different linear system solvers ([46]): • Cholesky decomposition (CD) having 0(l /6n 3 ) running time (CD can be used since the JSIM is positive definite), • LU decomposition with partial pivoting having 0(l/3n 3 ) running time, • Gauss Jordan elimination (GJE) with full pivoting, • Singular value decomposition (SVD). Although theoretically pivoting is not required, given that the JSIM is symmetric positive definite, it is possible that the numerical condition of the JSIM could degrade to a point where it becomes no longer positive definite. It is only in this case that algorithms such Chapter 5. Implementation and Experimental Verification 43 as LU decomposition with partial pivoting would be preferred over Cholesky decomposi-tion. The Numerical Recipes [46] implementations of these algorithms were used in our simulation program. Articulated-body forward dynamics Featherstone's articulated-body algorithm [19] was implemented. Due to the simplified class of robots being considered, the following simplifications were possible: • In 3D, spatial articulated-body inertias and spatial coordinate transformations are represented by 6 x 6 matrices (36 components). However, for the planar robots being ~ 4 ~A ~A ~A considered, only 9 components need actually be computed (It x 3, I{ x 4, It x 5, I{ 2 3, -A -A ~A -.4 ~A Ii 2,4' 1% 2,5' 'i 6 ,3' *i 6,4 a n " M 6,5)-• In 2D, spatial forces and motions can be represented using 3-vectors instead of 6-vectors (e.g. spatial forces can be represented by three components: the linear forces in x and y and a torque about z). Motion integration The main O D E integration method used is a fifth-order Runge-Kutta method with adap-tive step size control. Other numerical integration methods are also provided, including two integrators for stiff ODE: an adaptive step size Fourth-order Rosenbrock method and an adaptive step size semi-implicit method (Burlirsch-Stoer). These methods will not be discussed in this study. They were found to be much less efficient and difficult to adjust. The simulation program uses the Numerical Recipes [46] routines to perform all numerical integration. Chapter 5. Implementation and Experimental Verification 44 5.2 Results and observations: simulation of ill-conditioned systems 5.2.1 Performance of the forward dynamics In this section, simple two link chains are used to illustrate the effect of a poorly con-ditioned JSIM on individual computations of the forward dynamics. Ill-conditioning is obtained by increasing the length and mass of the second link while keeping the first link's dimensions constant. To illustrate the effect of joint configuration on the condition number, the forward dynamics test program prrrfdt was used. The forward dynamics were computed under the following conditions: gravity of 9.8 m/s2, no friction, no applied actuator torques, zero initial joint positions and velocities, unit fink widths, first link of unit length and mass, and second (distal) link length and mass were varied from 1.0 to 1.0e6. Cholesky decomposition was used to solve the linear system in the C R B M (note that similar results were obtained using LU decomposition with pivoting). For all computations, a double precision floating point representation (15 digits of precision) was used. Thus giving rise to relative errors (p) in the representation on the order of 10 - 1 6 . Figure 5.3 shows, for one particular two link chain (l2 = m2 = 102), the effect of the configuration on the condition number and computed accelerations. First, observe that although the configuration has a direct effect on the condition number (cnum = K(H_ )), the magnitude of the condition number remains on the order of 104 (such high condition numbers are actually not uncommon, for example the Stanford Arm is known to exhibit condition numbers as high as 11934 [3]). Next, observe that the size of the relative error between the C R B M and the A B M solution (Ar#i) follows the condition number and is on the order of 10 - 1 2 . As expected, this indicates that each forward dynamics implementation deals with ill-conditioning in a different way. It also illustrates the relationship between the condition number and the Chapter 5. Implementation and Experimental Verification 45 Figure 5.3: Effect of configuration on the condition number of 2 link chain. Note that the condition number is largest at the singularities in this case. The relative difference for joint 1 accelerations (Ar#i) is also shown; results for joint 2 accelerations are similar. Chapter 5. Implementation and Experimental Verification 46 Figure 5.4: Effect of link length relative error described in the previous chapter AR6>i « i f r 1 2 < i o 4 x i r r 1 6 = K{K )P. (5.1) Figure 5.4 summarizes the results obtained for several different lengths (and masses) of the distal link. Once again, this illustrates the direct relationship between the condition number and the relative error in the C R B M (note the logarithmic scales). Figure 5.5 depicts the computed accelerations for the A B M and C R B M for a particu-larly ill-conditioned linkage. This figure shows that the actual solution given by the A B M is smooth and that the solution given by the C R B M is not. In effect, this illustrates the poor behaviour of the C R B M in ill-conditioned situations and the exceptional behaviour of the A B M in these situations. In summary, these results show that the A B M does in fact deal better with ill-conditioned systems than the C R B M . Specifically, the size of the relative error in the C R B M solution is a function of the condition number of the JSIM and the floating point Chapter 5. Implementation and Experimental Verification 47 a) Results for link 2 length and mass = 4E4 -9.8 -9.80002 - • joint 1 -9.80004 - -(rad/s^) -9-80006---9.80008 -9.8001 J-joint 2 angle (rad) b) Results for link 2 length and mass = 5E4 -9.80001 -9.80002 * t joint 1 -980003 - • accel. -9.80004 • • (rad/s*2) -9.80005 - --9.80006 -9.80007 * joint 2 angle (rad) c) Results for link 2 length and mass = 1E5 joint 1 accel. (rad/s*2) joint 2 angle (rad) Figure 5.5: Comparison of calculated joint 1 accelerations (with l[2]=m[2]=le5). Note the smoothness of the A B M solution. Chapter 5. Implementation and Experimental Verification 48 representation being used. 5.2.2 Performance of the motion integration Now that we have illustrated the numerical properties of the forward dynamics, we con-sider how these properties carry over into the numerical motion integration. Results for the simulation of 10 seconds of motion of various two link planar arms are summarized in Figure 5.6. The forward dynamics were computed under similar conditions to those described in the previous section, with the exception that, in this case, each chain has the same total mass and length, only the ratio of the link lengths and masses is varied (this ensures that each system has the same total energy). Figure 5.6 displays the link length ratio (h/(h + h)), the order of magnitude of the maximum condition number (KMAX), the order of magnitude of the maximum relative error in individual evaluations of the C R B M forward dynamics (ArOi^^), the number of forward dynamics function evaluations for each method (TIABM, KCRBM), and the ratio rcRBM — ncRBM/n-ABM- For these simulations, an integration step size accuracy of le-6 was used. Notice that when the relative errors in the computed accelerations become greater than the desired step accuracy, the number of function evaluations in the CRB based simulation increases dramatically (in the last case, the integration for the C R B M failed to converge within a specified maximum number of iterations). This increase in the number of forward dynamics function evaluations is due to the adaptive step size integrator perceiving and reacting to the fluctuations in the computed accelerations (see figure 5.5). This effectively illustrates the phenomenon we refer to as "formulation stiffness". Although the results of figure 5.6 demonstrate an extreme case, similar problems can arise at much lower condition numbers. For example, in time critical simulations, single precision floating point representations (7 digits of precision) might be used in order to Chapter 5. Implementation and Experimental Verification 49 nABM nCRBM rCRBM le-1 le2 le-14 6526 6553 1.00 le-3 le6 le-10 31424 31469 1.00 le-5 le8 le-8 23787 23756 1.00 le-6 le9 le-7 23825 23888 1.00 le-7 lelO le-6 23797 23779 1.00 le-8 l e l l le-5 23812 24395 1.02 4e-9 l e l l le-5 23828 25392 1.07 2e-9 l e l l le-5 23792 31015 1.30 le-9 lel2 le-4 23792 64065 2.69 le-10 lel3 le-3 23822 — — Figure 5.6: Number of forward dynamics evaluations required to simulate 10 seconds of motion for various two link chains of increasing condition number. accelerate the mathematical computations. In such cases, the "formulation stiffness" would occur at much lower condition numbers. The situation would also be worse if a specific application required a very high level of step accuracy. Furthermore, the examples in this section dealt only with two link chains, however, preUminary studies have shown that these effects are made worse by the addition of more links to the chain. Finally, we note that we have performed similar experiments with a commercial dy-namics simulator (SD/FAST) where we compared simulations using Kane's formulation and using SD/FAST's 0(n) formulation which is similar to the A B M . While the specifics of the algorithms used in this code are not known to us, and the situation is complicated by symbolic simplification, we have observed the same difference in formulation stiffness between the two methods. In summary, these results have illustrated that the performance of a simulation al-gorithm does not depend only on the individual performances of the forward dynamics algorithm and the integration algorithm, but on the combined performance of these al-gorithms. Hence, the true measure of performance is given by the following expression Chapter 5. Implementation and Experimental Verification 50 for the total number of operations (n t o t aj): ntatal(ts, £) = nsingle x n F D ( t a , e) where: • ripr) is the total number of forward dynamics function evaluations required to ad-vance the solution by ts seconds with a step accuracy of e, • nsingle is the total number of operations required for a single forward dynamics function evaluation. 5.3 More simulation results 5.3.1 "Bad" implementations of the A B M As was discussed in the previous chapter, it is possible to implement the A B M algorithm so that its performance degrades in the presence of ill-conditioning. Figure 5.7 compares the well implemented A B M with an implementation of the A B M that computes the articulated body inertia using Ii = ( M i + ^M+l,t MiXi+ij) - X T - Xi+iti Si+i Mi+isi+i Accelerations were computed under conditions similar to those of figure 5.5. This result serves to confirm the observations of the previous chapter. 5.3.2 Conservation of energy In this section, using a simple example, we show that both the C R B M and A B M based implementations respect the conservation of energy principle. Figure 5.8 shows the evo-lution of the total energy for the simulation of a well-conditioned four link chain. This Chapter 5. Implementation and Experimental Verification 51 Results for link 2 length and mass - 1e5 -9.80000 • --9.80001 - • joint 1 -9.80002 - • accel. (rad/s"2) -9.80003 • • k -9.80004 • --9.80005 • • joint 2 angle (rad) Figure 5.7: Joint accelerations for "good" and "bad" implementations of the A B M (the smooth curve is the "good" implementation of the ABM). Chapter 5. Implementation and Experimental Verification 52 3304.85 -I . . . 3304.80 -total energy 3304.75 -^ 3304.70 < 3304.65 • a) ABM (step accuracy=1E-5) t 2 V 4 6 8 ^ 14 16 TO 20 t(s) 3304.85 i _ , 3304.80 -total energy 3304.75 • ( J ) 3304.70 < 3304.65 • b) CRBM (step accuracy=1E-5) —, h——; n - A _ r \ i i M f i k 2 V 4 6 8 10 ^ 1 2 14 16 18 20 t(s) Figure 5.8: Evolution of the total energy. figure shows that the total energy is conserved within a 5 decimal margin for the duration of the simulation for each forward dynamics method (this error is in keeping with the desired step accuracy). However, as the simulation advances, it seems to be more difficult to conserve energy. This behaviour is typical of any finite precision numerical simulation and is due to the propagation of the truncation error from step to step. For all simulations presented in this chapter energy was conserved in a similar fashion, regardless of ill-conditioning. The condition of the JSLM did not significantly affect the correctness of the results. Only the computational efficiency was degraded. Chapter 5. Implementation and Experimental Verification 53 5.3.3 Simulation of long chains This section presents results and observations related to the simulation of long chains. The competing forward dynamics methods (ABM and CRBM) are compared for regular chains varying in length from 2 to 14 links (each individual link is of equal length and mass). Observations are made regarding the execution time, the number of forward dynamics evaluations and the conservation of energy. Simulations were executed for 60 seconds of motion, with an integration step accuracy of 10 - 4 , for unactuated chains swinging from a horizontal position (initial joint positions and velocities are zero) in a constant gravity field (9.8m/s2). For each chain, a verification of the evolution of the joint angles in time confirmed that each forward dynamics method did indeed obtain the same joint angle trajectories to within the desired accuracy. Figure 5.9 presents the results obtained in each simulation, where: • n is the number of links, • Einit is the total energy (potential and kinetic) of the robot at the initial configu-ration, • teXec is the total execution time (timing results were obtained on a SUN 4/75 work-station) , • npo is the number of forward dynamics function evaluations, • Edev is the deviation (in percentage) of the total energy from Einu calculated using the average total energy Eavg as follows: Edev = 100 x (Eavg - Einit)/Einit. Figure 5.10 summarizes the results by showing the evolution of the number of forward dynamics function evaluations and the evolution of the total execution time as a function of the number of links. Chapter 5. Implementation and Experimental Verification 54 A B M C R B M n Einit{J) totalnpr) Edev (%) totalnpr) Edev (%) 2 3920 0.68 2377 0.011 0.67 2324 0.003 4 15680 2.18 3897 0.017 2.67 4362 0.036 6 35280 4.85 5818 0.015 5.92 5912 0.021 8 62720 8.75 7668 0.028 10.62 6916 0.018 10 98000 11.95 8576 0.019 19.47 8892 0.030 12 141120 15.18 9234 0.025 28.85 9746 0.019 14 192080 20.78 10497 0.016 39.25 10039 0.015 Figure 5.9: Summary of results for the simulation of long regular chains. Observations As the measure of energy deviation indicates, the simulations seem to conserve energy quite well regardless of the number of links. For these problems, the implementation of the A B M seems to require less execution time than the C R B M even for small values of n. This would indicate the need for further optimization of the C R B M code. However, it is also possible that the A B M is better suited to this particular type of planar problem than the C R B M . Further study is indicated. As expected, the total execution time also increases with the number of links. Fur-thermore, in almost all the test cases, real-time performance is obtained. The number of calls to the forward dynamics function seems to increase with the number of links in a logarithmic fashion (see figure 5.10). This would indicate that as the number of links increases, the set of differential equations becomes more difficult to solve to the desired accuracy (smaller steps must be taken). This measure suggests that the operational stiff-ness of the differential equations increases with the number of links. This stiffness could simply be attributed to the increasing complexity of the resulting physical motion and should not necessarily be attributed to the forward dynamics formulation. Chapter 5. Implementation and Experimental Verification 55 Figure 5.10: Simulation of long chains. Chapter 5. Implementation and Experimental Verification 56 In summary, the increase in running time stems from two sources of additional work: 1. as n increases, each call to the forward dynamics function is longer to evaluate (complexity: 0(n) for A B M and 0(n 3) for C R B M ) , 2. as n increases, more calls to the forward dynamics function are required (complex-ity: 0(/(n)), where f(n) appears to be logarithmic). This second source of work explains the fact that the running time for the A B M simulation does not seem to progress linearly with the number of links. Based on these observations, the effective running time for the A B M simulation can be expressed as 0(n x f(n)) and 0(n3 x f(n)) for the C R B M simulation. Chapter 6 Extensions: Contact Simulation One of the basic functions of a robot manipulator is to act on and modify its environ-ment. In general, to accomplish these goals, the robot must come in contact with its environment. When the end-effector of a robot comes in contact with its environment, constraints are imposed on the robot's motion. The simulation of interactions between the robot and its environment is a natural extension to a robot simulator and is currently one of the most challenging problems of rigid body dynamic simulation. This chapter presents an overview of the principal problems involved in contact simulation. The contact simulation problem can be divided into the following subproblems: con-tact detection and contact resolution. Contact resolution can be further divided into impact simulation and continuous contact simulation. In the following sections each of these issues will be considered individually. 6.1 Notation This chapter will first consider the general problem of resolving the interactions between two single rigid bodies interacting at a finite number of contact points. Extending the proposed solutions to a robot manipulator interacting with a rigid body is straight for-ward and will also be considered. However, extending solutions to a non-finite number of contact points is an open problem. The following notation is introduced to deal with contact simulation between two rigid bodies B\ and B2 colliding at a single contact point V (see figure 6.1). 57 Chapter 6. Extensions: Contact Simulation 58 Figure 6.1: Contact simulation notation. • n is the unit contact normal vector directed from B2 to B\, and let the right subscript N denote the normal component of a vector with respect to h (i.e. vpj = nTv) and T denote the tangential components of a vector with respect to n. • Mi and M2 are the spatial inertias of Bi and B2. • vi~ and v2~ are the known spatial velocities of B\ and B2 before the impact, and d\+ and v2+ are the unknown spatial velocities of Bi and B2 after the impact. • u~ = v2~ — v{~ is the relative velocity of the contact point before impact, and v+ — v2+ — vy+ is the relative velocity of the contact point after impact. • Adi = Vi+ — Vi~ is the absolute velocity change of body Bi due to the impact, and Av — v+ — v~ is the change in relative velocity due to the impact. Chapter 6. Extensions: Contact Simulation 59 6.2 Contact detection From a geometrical perspective the contact detection problem can be defined as the problem of determining all points of intersection between a set of objects at a time t given the current position and orientation of the objects. In the context of contact simulation, the specific time at which two objects collide tc is also of interest. For a relatively complex environment, three-dimensional collision detection can be a very difficult and time consuming process. Resolution of this problem is regarded by many as the performance bottleneck of the contact simulation process (see simulation results given by Hahn [29]). Conventional methods result in 0(n2) complexity, where n is the number of bodies [43]. However, more efficient algorithms based on efficient distance computation methods have recently been proposed [27, 37]. Cremer [16] uses a simple iterative approach in order to determine the time of col-lision. By assuming that the integration time step is small enough so that no collisions are missed, the time of collision is obtained by performing a binary search of the time step. The iterative search can be improved by making use of the current acceleration of the contact point to better divide the search space [7]. However, assuming that the integration time step is small enough so that no collisions are missed may impose the choice of an unnecesseraly small integration time step. More precise methods perform a search for intersections between the four-dimensional space-time swept volumes of the objects in order to determine the exact time and point of collision [26, 13, 42]. For the contact simulation process, it is common to assume that all objects can be modelled as convex polyhedra (non-convex objects can be decomposed into unions of convex polyhedra). Reducing the number of objects tested to a minimum is perhaps the most basic way of reducing the cost of collision detection. A spatial segmentation based on bounding volumes can be used to determine objects that are near to another. Chapter 6. Extensions: Contact Simulation 60 Also, knowledge of velocity and acceleration bounds can be used to reduce the number of collision tests. Using such bounds, it is possible to construct swept bounding volumes which can be tested for intersection in a preprocessing phase [42]. Once contact points have been determined, they must be organized into an efficient data structure for the simulation process. Mirtich and Canny [42] propose an interesting approach based on constructing a "collision heap" which gives a prioritization of the contact events in the simulation. This prioritization of contact events can facilitate the choice of appropriate integration time steps. 6.3 Friction modelling When two bodies are in contact, friction creates a force that opposes relative motion. According to Coulomb's law, the friction force is proportional to the normal force acting at the contact point. The Coulomb friction model is expressed as follows: The first condition is referred to as dynamic friction and the second is referred to as static friction, u-r is the tangential velocity of the point of contact, fT is the friction force, is the component of the total force along the normal, us is the static friction coefficient and \id is the dynamic friction coefficient. The coulomb friction model defines an exact friction force only in the dynamic case. In the case of static friction, all that is known is an upper bound on the friction force (cf. friction cone). In the following, we will assume Although the Coulomb friction model is approximate, it is the most widely used. As will be discussed in the following sections, modelling friction in contact simulation is a very difficult task. (6.1) (6.2) Chapter 6. Extensions: Contact Simulation 61 6.4 Impact simulation Impact simulation is defined as the simidation of the motion interactions of two colliding bodies for which the point of contact has a non-zero relative normal velocity Au^r- This section reviews some of the popular methods for impact simulation. 6.4.1 Penalty methods Penalty methods 1 are approximate methods which handle collisions by applying a force at the contact point which penalizes interpenetration. In general, the penalty force is obtained using a spring model. A spring is inserted in the normal direction between the colliding bodies at the point of contact. If d is the interpenetration depth then the resulting penalty force is f = kd (linear spring model). To be effective, penalty methods must necessarily permit a certain degree of interpenetration.2 The spring constant k is chosen ad hoc so as to be large enough to prevent too much interpenetration. However, the larger the spring constant, the stiffer the resulting differential equations of motion and the smaller the required integration step size. Due to the approximate nature of the numerical simulation process, a spring contact model cannot guarantee the conservation of energy. In order to obtain a more stable contact model and in order to simulate frictional effects during impact, it is common to add a damping component to the spring model. In [1], Anderson uses such a spring-damper system to simulate force feedback in a telerobotics application. The elasticity of an impact can be specified by using two different spring constants in order to simulate the loss of energy during a collision. One for the compression phase (Aujv < 0) and the other for the relaxation phase (Au# > 0), e.g. fcre/ax = ekcompress 'This approach is based on penalty function methods used in nonlinear optimization. 2 Certain methods minimize interpenetration by extending the boundaries of the bodies so as to begin applying a penalty force when two bodies come within a distance e of each other. However, this does not necessarily result in a more realistic simulation. Chapter 6. Extensions: Contact Simulation 62 where 0 < e < 1 is the coefficient of "restitution" [43]. Another variation of the spring model consists of using a non-linear spring constant [21]. Furthermore, by adding tangential spring-damping components, one can also simulate the effects of frictional forces. An in depth analysis of penalty methods can be found in Baraff [7]. 6.4.2 Analytical methods Analytical methods generally assume that the impact between two rigid bodies occurs instantaneously. The collision results in a discontinuity in the velocity of the bodies. When an impact occurs, the motion integration routine must be stopped and reset with new post-impact initial velocity conditions. This discontinuity in velocity is a result of an impulsive constraint force that acts normal to the contact surface Av1N = -If^ (6.3) (n Iin' A u 2 i V = - - # T ) - ( 6 - 4 ) (n I2n' The impulsive constraint force is computed using the coefficient of restitution. The coefficient of restitution (or elasticity) is a function of the materials of the two bodies in contact. There are two principal interpretations of this coefficient: Newton's and Poisson's. This section, presents rigid body simulation systems based on each of these interpretations. Good references on the topic of impact modelling include [33, 10, 59]. Analytical methods based on Newton's law of restitution Newton's law of restitution states that the coefficient of restitution relates the relative normal velocity after impact to the relative normal velocity before impact: VN = ~evN- (6-5) Chapter 6. Extensions: Contact Simulation 63 Using this relationship and equations 6.3 and 6.4 relating velocity variations to applied impulses, a linear set of equations is obtained. For the frictionless case, this set of equations can be directly resolved to obtain the change in velocity [16, 29]. In [29] the equations include a Coulomb friction model for tangential impulses. Tan-gential impulses are related to normal or constraint impulses in a fashion similar to equations 6.1 and 6.2. The solution algorithm first assumes that the bodies stick, and if the impulse is larger than the static friction bound then it is assumed that the bodies are sliding and the frictional impulses are determined. This method assumes the relative tangential velocity does not change during a collision (the problem of reversed sliding [59] is not considered). This method is similar to that used by Moore and Wilhelms [43]. Baraff [7] uses Newton's interpretation of the coefficient of restitution to obtain the following constraints: (vN+ + evN-)>0, (6.6) PN > 0, (6.7) PN{vN+ + evN~) = 0. (6.8) These constraints can be used to obtain a L C P formulation similar to that used for continuous contact simulation and can be solved by quadratic programming techniques. Analytical methods based on Poisson's hypothesis of restitution Poisson's hypothesis of restitution distinguishes between a compression and a relaxation phase of impact. The hypothesis states that the coefficient of restitution relates the total normal impulse to the normal impulse at maximum compression PN = (1 + e)PNC (6-9) where pNc is the total normal impulse up to the point of maximum compression and pN is the total normal impulse due to the whole impact. Chapter 6. Extensions: Contact Simulation 64 As with the previous analytical methods, the velocity variations are related to the total impulse by equations 6.3 and 6.4. In the frictionless case, Poisson's hypothesis can be used to obtain a closed form expression for the total normal impulse [20]. In this case, the resulting equations are the same as for Newton's interpretation of restitution. Mirtich and Canny use a detailed impact model based on Poisson's hypothesis. This model incorporates a complete friction model that deals with the transitions between sticking and sliding during impact. The method consists of numerically integrating the velocity during impact using the cumulative normal impulse as the integration parameter (instead of time). When the velocity reaches zero, the maximum compression has been reached and the total impulse at compression is obtained. Then, by multiplying by (1 + e) the total impulse at the end of the impact is obtained. This represents the value for the "time" of end of collision, integration continues to this point and the variation in velocity is obtained. This approach is based in part on the work of Keller [33]. Incorporating friction Dealing with friction in the impact model is a very difficult problem. Wang and Mason present a treatment of impact with friction for planar collisions involving single contact points. One of their main findings is that Poisson's hypothesis is better than Newton's law of restitution for modelling impact with friction (they show that Newton's law of restitution does not always conserve energy). This result in based on the fact that Pois-son's hypothesis is expressed as a dynamic law and is consistent with energy conservation. Brach [10] states that the fundamental problem in modelling impacts with friction is the interpretation of the coefficient of restitution and friciton. Brach uses Newton's law of restitution and resolves the increase in energy by using a slightly different interpretation of the coefficient of friction and restitution. Baraff discusses an approximate method for simulating impact with friction [7]. Chapter 6. Extensions: Contact Simulation 65 Simultaneous impacts It is very difficult to model multiple-contact point collisions. In general, most models assume that no two collisions occur at the same time and treat simultaneous collisions serially. Cremer [16] presents some of the advantages and disadvantages of using a simul-taneous model of collisions. In general, the problem of resolving simultaneous collisions with friction is currently unsolved. Mult ibody systems Extending impact simulation to include collisions with multibody systems, generally involves propagating the velocity discontinuities through each joint of the multibody system. These issues are discussed in [6] and [43]. Summary Even though resolving impacts using analytical methods may seem like a lot of work, it is often more efficient than a spring model since it only has to be computed once for each impact [43] (spring models involve several evaluations over many small steps of integration to give realistic results). 6.4.3 Finite element methods Finite element simulation methods (FEM) consist of simulating contacts at a very de-tailed (elementary) level. These methods can be used to study the stresses and strains oc-curring during the compression and restitution phase of a collision, however such methods are too computationally expensive for real-time simulation. Furthermore, it is generally believed that such detail is uneccessary and that simple empirical rules can sufficiently describe the collision behaviour [42]. Chapter 6. Extensions: Contact Simulation 66 6.5 Simulation of continuous contact Baraff defines the continuous contact simulation problem as follows [8]: "Given a sys-tem of rigid bodies, contacting at n points and nowhere colliding, with known spatial and velocity variables, calculate the contact forces that should arise between the bodies to prevent interpenetration (...) the contact forces must also obey constraints due to friction." In this section, several approaches to continuous contact simulation will be reviewed. 6.5.1 Penalty methods Penalty methods similar to those used for impact modeling can be used in continuous contact simulation. The main advantage of using such models is that they are very simple to implement. The use of such a model during impact and continuous contact may also ensure smoother transitions between the impact and continuous contact motions. For certain applications, these types of models may generate results which give an acceptable level of realism. However, many applications demand a greater level of realism; hence the interest in analytical methods. 6.5.2 Analytical methods Analytical methods or constraint-based methods impose kinematic constraints on the configuration of contacting bodies such that interpenetration is strictly prohibited. Forces at a contact point are separated into normal forces and tangential forces. Normal or constraint forces are workless forces which prevent interpenetration, while tangential forces or friction forces oppose motion along the contact surface. For two bodies in continuous contact with no friction, the following local motion Chapter 6. Extensions: Contact Simulation 67 constraints can be written for each contact point: d{vN)/dt > 0, (6.10) IN > 0 (6.11) fNd(vN)/dt = 0. (6.12) These constraints define a convex linear complementary problem (LCP) which may be solved by quadratic programming (QP) techniques. Lotstedt [38] was one of the first to formulate the contact simulation problem in such a manner. Recently, BarafF [7] has used such a formulation to simulate problems with and without friction. BarafF studied the inconsistencies related to frictional contact simulation [47, 59] and showed that detennining whether a configuration is consistent is NP-hard. Recently, Lee et al. [35] showed how Baraff's contact simulation methods can be ex-tended to the context of a robot simulation system. An L C P formulation is also used by Featherstone [20] to solve the problem of constrained robot dynamic simulation. Feath-erstone combines the dynamic equations of the robot manipulator with n simultaneous quadratic equations in the unknown contact forces. For an equality constrained problem, constraint forces can be computed by solving a linear system of equations. Such a formulation is used by Lilly [36] to simulate the motion of a robot whose end-effector is constrained by continuous contact. Lilly uses an operational space formulation of the robot's inertia matrix to obtain an expression for contact constraints which results in a linear system of algebraic equations. This approach is similar to the contact simulation methods used by Witkin et al. [61] (this work also considers multiple contact points). For a robot manipulator constrained by a contact with the environment, the con-straining forces can be much larger than the joint forces resulting in a system of stiff Chapter 6. Extensions: Contact Simulation 68 ordinary differential equations. Therefore, equality constrained contact problems are very well suited to D A E formulations. D A E formulations of continuous contact problems can be found in [40, 21]. Baiardi et al. [6] present the differential algebraic equations of motion for a kinematic chain in contact with a hard motionless surface. Three formulations are obtained for constraints under the following conditions: friction free, with static friction and with dynamic friction. The transitions between each of these conditions are also explored. One of the drawbacks with analytical approaches is in dealing with changing contact profiles. Constraints must be monitored in order to keep track of the changing geometry of the contact surface. Another disadvantage is that impacts must be dealt with separately. 6.5.3 Impulse based methods Mirtich and Canny [42] propose an interesting approach to contact simulation called impulse-based simulation which resolves all types of interactions as impacts. Under this approach, continuous contact is modelled as a series of numerous impacts called micro-collisions. This approach is both conceptually simple and computationally efficient. The underlying assumption is that collision forces in nature are based on local prop-erties like contact velocity, and not on the global state configuration of bodies in contact. Mirtich and Canny use an impulse simulation algorithm based on [33]. One of the advantages of this approach to contact simulation is that it unifies all types of contacts (bouncing, rolling, sliding, etc..) under a single model. No artificial boundary is created between the different modes of motion. It is also claimed that this approach results in physically accurate behaviour. Chapter 7 Conclusion In this thesis, we have explored several key issues involved in dynamic simulation of robot manipulators. We have showed how a good notation can play an important role in the formulation and understanding of the equations of motion of robot dynamics. We have identified an important phenomenon which we call "formulation stiffness". Different formulations of the forward dynamics problem were shown to deal differently with ill-conditioned joint-space inertias. In particular, we have shown that the articulated-body method can be less sensitive to ill-conditioning, while the composite rigid body method produces errors that are proportional to the condition number of the joint space inertia matrix. This is important because the errors in the C R B M solution can incapacitate explicit adaptive-step size integration methods by requiring very small integration step sizes. In this light, we demonstrated with experimental results that the articulated-body method is superior to the composite rigid body method for ill-conditioned problems. We thus highlight the importance of measuring the performance of simulation algorithms by taking into account both the computational complexity of the forward dynamics and its numerical properties. Awareness of such issues is an essential step in properly identifying the performance and limitations of computer simulation. Of particular interest to the author is the possible contribution of robot dynamic simulation to the development of advanced telerobotics operator interfaces. Recently, many have showed how graphical simulations can be used to improve operator interfaces to telerobotics control systems. 69 Chapter 7. Conclusion 70 Graphical robot simulations can be used to create what are known as predictive displays [53, 9]. Predictive displays are used to bridge the time delay gap in situations where the remote robot is situated far from the operator (e.g. space telerobotics can involve time delays of the order of 10 seconds). Experimental results indicate that the use of predictive displays improves the performance of the human operator [53, 9, 14, 31] (in comparison with simple "do and wait" methods). A graphical simulation of the robot and its environment may also serve to complement available camera views of the worksite [12]. Recently, Funda and Paul have shown how simulation can be used to create a preview display for teleprogramming [24, 45]. In addition, a robot simulation interface can include elements such as virtual tools and fixtures which can further aid robot operators [51]. Most of the current systems used to enhance telerobotics interfaces incorporate only a kinematic simulation of the robot manipulator. Therefore, it is of particular interest to see how dynamic simulation, including simulation of contacts with the environment, could be used to aid telerobotics. Bibliography [I] R.J. Anderson. Teleoperation with virtual force feedback. In SPIE Int. Symposium on Optical Tools for Manufacturing and Advanced Automation, pages 32-39, Boston, M A , 1993. [2] T . Andrzejewski, H.G. Bock, E . Eich, and R. von Schwerin. Recent advances in the numerical integration of multibody systems. In W. Schiehlen, editor, Advanced Multibody System Dynamics, pages 127-151, Kluwer Academic Publishers, 1993. [3] J. Angeles and 0. Ma. Dynamic simulation of n-axis serial robotic manipulators using a natural orthogonal complement. International Journal of Robotics Research, 7(5):32-47, 1988. [4] W.W. Armstrong. Recursive solution to the equations of motion of an n-link ma-nipulator. In Proceedings of the 5th World Congress on Theory of Machines and Mechanisms, pages 1343-1346, Montreal, 1979. [5] U.M. Ascher, D.K. Pai, and B.P. Cloutier. Forward Dynamics, Elimination Methods, and Formulation Stiffness in Robot Simulation. Technical Report, University of British Columbia, 1995. [6] P. Baiardi, G. Cannata, G. Casalino, and P. Pagano. Modelling contact phenomena within the dynamic simulation of advanced robotic structures. In Proceedings of the IEEE International Conference on Robotics and Automation, pages 196-203, 1993. [7] D. BarafF. Dynamic Simulation of Non-Penetrating Rigid Bodies. PhD thesis, Cor-nell University, 1992. [8] D. BarafF. Issues in computing contact forces for non-penetrating rigid bodies. Al-gorithmica, 10:292-352, 1993. [9] A . K . Bejczy, W.S. Kim, and S.C. Venema. The phantom robot: predictive dis-plays for teleoperation with time delay. In Proceedings of the IEEE International Conference on Robotics and Automation, pages 546-551, 1990. [10] R . M . Brach. Rigid body collisions. Journal of Applied Mechanics, 56:133-138, 1989. [II] H. Brandl, R. Johanni, and M . Otter. A very efficient algorithm for the simula-tion of robots and similar multibody systems without inversion of the mass matrix. 71 Bibliography 72 In Proceedings of IFAC/IFIP/IMACS International Symposium on the Theory of Robots, Vienna, Austria, December 1986. [12] R.A. Browse and S.A. Little. The effectiveness of real-time graphic simulation in telerobotics. In Proceedings of the 1991 IEEE International Conference on Robotics and Automation, pages 895-898, 1991. [13] J. Canny. Collision detection for moving polyhedra. IEEE Transactions on Pattern Analysis and Machine Intelligence, 200-209, 1986. [14] L. Conway, R. Volz, and M. Walker. Tele-autonomous systems: methods and archi-tectures for intermingling autonomous and telerobotic technology. In Proceedings of the 1981 IEEE International Conference on Robotics and Automation, pages 1121-1130, 1987. [15] J. Craig. Introduction to Robotics: Mechanics and Control. MIT Press, Cambridge, Mass., 1989. [16] J.F. Cremer. An Architecture for General Purpose Physical System Simulation. PhD thesis, Cornell University, 1989. [17] P.E. Dupont. The effect of friction on the forward dynamics problem. International Journal of Robotics Research, 12(2):164-179, 1993. [18] R.E. Ellis, O .M. Ismaeil, and I.H. Carmichael. Numerical stability of forward-dynamics algorithms. In Proceedings of the IEEE Conference on Robotics and Au-tomation, Nice, France, 1992. [19] R. Featherstone. The calculation of robot dynamics using articulated-body inertias. International Journal of Robotics Research, 2(l):13-30, 1983. [20] R. Featherstone. Robot Dynamics Algorithms. Kluwer Academic Publishers, Nor-well, MA, 1987. [21] G. Ferretti, C . Maffezzoni, and G. Magnani. Dynamic simulation of robots interact-ing with stiff contact surfaces. Transactions of The Society for Computer Simulation, 9(l):l-24, 1992. [22] A. Fijany and A . K . Bejczy. Techniques for parallel computation of mechanical ma-nipulator dynamics, part ii: forward dynamics. Control and Dynamic Systems, 40:357^110, 1991. [23] A. Fijany and R.E . Scheid. Efficient conjugate gradient algorithms for computation of the manipulator forward dynamics. In Proceedings of NASA Conference on Space Telerobotics, Pasadena, C A , January 1989. Bibliography 73 J. Flmda and R.P. Paul. Efficient control of a robotic system for time-delayed environments. In Proceedings of the IEEE Conference on Robotics and Automation, pages 219-224, 1991. W . G . Gear. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, 1971. E . G . Gilbert and S.M. Hong. A new algorithm for detecting the collision of moving objects. In Proceedings of the 1989 IEEE International Conference on Robotics and Automation, pages 8-14, 1989. E . G . Gilbert, D.W. Johnson, and S.S. Keerthi. A fast procedure for computing the distance between complex objects in three-dimensional space. IEEE Journal of Robotics and Automation, 4(2):193-203, 1988. G.H. Golub and C.F. Van Loan. Matrix Computations. John Hopkins University Press, Baltimore, 1989. J.K. Hahn. Realistic animation of rigid bodies. In Computer Graphics (ACM), pages 299-308, 1988. E. Hairer and G. Wanner. Solving Ordinary Differential Equations II - Stiff and differential Algebraic Problems. Springer-Verlag, 1991. G. Hirzinger, K. Landzettel, and J. Heindl. Rotex - space telerobotic flight exper-iment. In Telemanipulator Technology and Space Telerobotics (SPIE Vol. 2057), pages 51-72, 1993. A. Jain. Unified formulation of dynamics for serial rigid multibody systems. Journal of Guidance, Control, and Dynamics, 14(3):531-542, 1991. J.B. Keller. Impact with friction. Journal of Applied Mechanics, 53:1-4, 1986. C. Lanczos. The Variational Principles of Mechanics. University of Toronto Press, 1970. P.U. Lee, D.G. Ruspini, and 0. Khatib. Dynamic simulation of interactive robotic environment. In IEEE International Conference on Robotics and Automation, pages 1147-1152, 1994. K.W. Lilly. Efficient Dynamic Simulation of Robotic Mechanisms. Kluwer Academic Publishers, Norwell, MA, 1993. M.C. Lin. Efficient Collision Detection for Animation and Robotics. PhD thesis, University of California, Berkeley, 1993. Bibliography 74 P. Lodstedt. Numerical simulation of time-dependent contact and friction problems in rigid body mechanics. SI AM J. Sci. Stat. Comput., 5(2):370-393, 1984. J.Y.S. Luh, M.W. Walker, and R.P.C. Paul. On-line computational scheme for mechanical manipulators. Journal of Dynamic Systems, Measurement and Control 102:69-76, 1980. N.H. McClamroch. Singular systems of differential equations as dynamic models for constrained robot systems. In IEEE International Conference on Robotics and Automation, pages 21-28, 1986. J.P. Meijaard. A comparison of numerical integration methods with a view to fast simulation of mechanical systems. In E . J . Haug and R.C. Deyo, editors, Real-Time Integration Methods for Mechanical Systems Simulation, Springer-Verlag, 1990. B. Mirtich and J. Canny. Impulse-based, real time dynamic simulation. In Euro-graphics Proceedings, pages 1-13, 1994. M . Moore and J. Wilhelms. Collision detection and response for computer anima-tion. In Computer Graphics (ACM), pages 289-298, 1988. J.J. Murray and C P . Neuman. Organizing customized robot dynamics algorithms for efficient numerical evaluation. IEEE Transactions on Systems, Man, and Cyber-netics, 18(1):115-125, 1988. R. Paul, T. Lindsay, and C. Sayers. Time delay insensitive teleoperation. In Pro-ceedings of the 1992 IE EE/RSJ International Conference on Intelligent Robots and Systems, pages 247-254, 1992. W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical Recipes in C: the Art of Scientific Computing. Cambridge University Press, second edition, 1992. V . T . Rajan, R. Burridge, and J.T. Schwartz. Dynamics of a rigid body in fric-tional contact with rigid walls. In IEEE International Conference on Robotics and Automation, pages 671-677, 1987. R.E. Roberson and R.F. Schwertassek. Introduction to the Dynamics of Multibody Systems. Springer-Verlag, 1987. G Rodriguez. Kalman filtering, smoothing and recursive robot arm forward and inverse dynamics. IEEE Journal of Robotics and Automation, 3(6):624-639, 1987. Bibliography 75 [50] G. Rodriguez, A. Jain, and K. Kreutz-Delgado. A spatial operator algebra for manipulator modeling and control. The International Journal of Robotics Research, 10(4):371-381, 1991. [51] L.B. Rosenberg. Virtual fixtures as tools to enhance operator performance in telep-resence environments. SPIE, 2057:10-21, 1993. [52] R. Schwertassek and W. Rulka. Aspects of efficient and reliable multibody system simulation. In E.J. Haug and R.C. Deyo, editors, Real-Time Integration Methods for Mechanical System Simulation, pages 55-96, Springer-Verlag, Berlin, 1990. [53] T .B . Sheridan. Telerobotics. Automatica, 25(4):487-507, 1989. [54] W . M . Silver. On the equivalence of lagrangian and newton-euler dynamics for ma-nipulators. International Journal of Robotics Research, l(2):60-70, 1982. [55] M . et al. Sofer. Ode formulations for multibody dynamics: numerical aspects. In W. Schiehlen, editor, Advanced Multibody System Dynamics, pages 397-402, Kluwer Academic Publishers, 1993. [56] G. Strang. Linear Algebra and its Applications. Harcourt Brace Jovanovich College Publishers, Orlando, Florida, 1988. [57] A . F . Vereshchagin. Computer simulation of the dynamics of complicated mecha-nisms of robot-manipulators. Engineering Cybernetics, 6:65-70, 1974. [58] M.W. Walker and D.E. Orin. Efficient dynamic computer simulation of robotic mechanisms. Journal of Dynamic Systems, Measurement, and Control, 104:205-211, 1982. [59] Y . Wang and M.T. Mason. Two-dimensional rigid-body collisions with friction. Journal of Applied Mechanics, 59:635-642, 1992. [60] J .H. Wilkinson. Rounding Errors in Algebraic Processes. Prentice-Hall, 1963. [61] A. et al. Witkin. Interactive dynamics. In Symposium on INteractive 3D Graphics (ACM), pages 11-21, 1990. 

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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0051428/manifest

Comment

Related Items