An Architecture for General Purpose
Physical System Simulation--
Integrating Geometry, Dynamics, and
Control
James F Cremer
Ph. D. Thesis
TR 89-987
April 1989
Department of Computer Science
Cornell University
Ithaca, NY 14853-7501
AN ARCHITECTURE FOR GENERAL PURPOSE
PHYSICAL SYSTEM SIMULATION
INTEGRATING GEOMETRY, DYNAMICS, AND
CONTROL
A Dissertation
Presented to the Faculty of the Graduate School
of Cornell University
in Partial Fulfillment of the Requirements for the Degree of
Doctor of Philosophy
by
James F. Cremer
May 1989
Oc James F. Cremer 1989
ALL RIGHTS RESERVED
AN ARCHITECTURE FOR GENERAL PURPOSE PHYSICAL SYSTEM
SIMULATION --H INTEGRATING GEOMETRY, DYNAMICS, AND
CONTROL
James F. Cremer, Ph.D.
Cornell University 1989
Simulation of physical systems has long been of interest to scientists and en-
gineers, and significant efforts have been directed toward the development of
general purpose computer aided design and analysis systems. To date, however,
success has been largely limited to the production of tools suited only for particu-
lar aspects of design: computer aided design systems have primarily emphasized
specification of geometry; simulation systems from the mechanical engineering
field have concentrated mainly on formulation and integration of an unchang-
ing set of equations describing object behavior; and work by computer graphics
and animation researchers has been aimed at producing good-looking animations
without much regard for whether the generated motions were physically realistic.
Flexible systems integrating all aspects of design and analysis have not yet been
built.
This thesis addresses the issues involved in developing fundamentally more
powerful simulation systems. A system architecture for general purpose physical
system simulation is proposed, and a prototype implementation, the Newton
system, is described. The architecture is based on a rich model-based object
representation and provides a level of automatic analysis that encourages the
kind of experimentation necessary for successful design. In particular, it is argued
that by using geometric modeling techniques to include a full description of object
geometry, previously difficult-to-incorporate simulation system features, such as
collision detection and resolution, can be handled much more naturally.
The Newton architecture incorporates a uniform exceptzonal event handling
mechanism that allows the system to respond to a variety of simulation events
that necessitate modification of object behavior descriptions. Combined with a
general method for the formulation of object motion equations, the event han-
dling mechanism supports automatic handling of collisions, changing contact
relationships, control program state changes, and other events that can cause
discontinnities in object motions.
Biographical Sketch
Jim Cremer was born in San Francisco, California, on September 12, 19(30. After
living in thirteen interesting places in seventeen years, he came to Cornell in the
fall of 1978. He earned a B.S. in Computer Science in 1982. While pursuing
his graduate studies at Cornell, he has been employed by the Computer Science
Department's modeling and simulation research group. He received a M.S. degree
in Computer Science in August, 1988.
in
To Anne
lv
Acknowledgements
I am very grateful to my advisor, John Hopcroft, for providing persistent en-
couragement, constant support, and valuable insight. He is responsible for the
stimulating research environment that I have enjoyed for the past three years.
Special thanks are also due to Alberto Segre, Peter Kahn, and Robert Con-
nelly, for serving on my thesis committee
Many members of the modeling and simulation research group have supplied
expertise and friendship. In particular, James Stewart, Anne Neirynck, Dinesh
Pal, Daniela Rus, and Michael Karasick have made important contributions to
this work. Christoph Hoffmann, Bill Bouma, and Joe Kearney provided addi-
tional ideas and encouragement.
My graduate career would not have been successful without the strong early
support of Joseph Bates and Robert Constable.
Finally, I thank my family for their love and for their patience.
v
Table of Contents
2
3
4
Introduction
1.1 Organization of the thesis
1.2 Review of related research
1.2.1 Mechanical engineering work on physical system simulation
1.2.2 Work in the computer graphics and animation fields
1.2.3 Other related work
Summary
1.3
General Purpose Simulation and Analysis Systems: Architectural
Issues
2.1			Simulation architectures
2.2			The Newton system
2.3			An example simulation			. .
2.4			Newton architecture overview
2.4.1			Definition and representation module
2.4.2			Analysis module
2.4.3			Report module
Representation and Definition
3.1 Objects and hinges
3.2 Models
3.3 Quantities
3.4 Equations . .
Newton's Analysis Module: Dynamics
4.1 Introduction
4.2 Basic dynamic analysis
4.3 Motion equation formulation
4.4 Equation storage, collection, and modification
4.5 Integration methods
vi
5
7
8
11
14
15
17
17
20
21
26
26
29
30
32
32
37
41
43
46
46
47
49
57
59
5
6
7
4.5.1 Stability
4.6			Dependence
4.7			Other formulation schemes
Newton's Analysis Module: Handling of Exceptional Events
5.1			Handling of exceptional events
5.1.1 Details of the event handling mechanism
5.2			Impact.
5.3			Contact
5.4			Friction			.
Motion Control
6.1			Introduction
6.2			Control quantities in Newton
6.3			Basic control of motion . . .			. .
6.3.1 The direct control approach
6.3.2 The control equation approach
6.3.3 Comparison of the two control mechanisms
6.3.4 Underdetermined equations
6.4			Control and exceptional events
6.5			Control programming
6.6			Control examples
6.6.1 Simple low-level controllers
6.6.2 Control of biped walking
6.6.3 Control of redundant manipulators
Conclusions and Continuing Research
7.1 Contributions
7.2 Continuing research . .
Bibliography
vii
63
65
69
71
71
74
79
85
94
98
98
99
101
102
103
104
105
106
106
107
107
108
113
120
120
123
129
List of Figures
1.1 Different hand configurations
2.1
2.2
2.3
2.4
2.5
2.6
3.1
3.2
3.3
3.4
3.5
4.1
4.2
4.3
4.4
4.5
4?6
4.7
4.8
.3
An example simulation
Object class definitions for the example simulation
Scene specification for the example simulation
Basic Newton architecture diagram
Newton `s top-level loop			.
Redundant robot arm			simulation and graphs
Basic structure of object and hinge classes
Composition graph for simple arm model
Structure graphs for simple arm model
Quantities
Equation schemata and equations
22
23
24
27
29
31
34
36
36
42
45
48
53
54
55
56
60
62
Newton's top level
Pin (revolute) hinge
Prismatic hinge . .
Printed form of automatically generated equations for example of
Section 2.3.
Two closed-loop linkages .
Simple integration method
Fourth-order Runge-Kutta integration formula.
Fourth-order Adams-Moulton predictor-corrector integration for-
mula62
4.9			Four-bar closed loop linkage with pin hinges			and bottom object
fixed			66
4.10			Four-bar linkage in bifurcation position and			two possible later
configurations. 			67
4.11			A critical form			68
5.1			Newton's event handler interface			75
viii
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
5.10
5.11
6.1
6.2
6.3
6.4
6.5
6.6
6.7
6.8
6.9
6.10
Impact resolution scheme
Different collision interpretations. . .
Collision requiring large sequence of impact equation solutions
Collision for which the simplest model formulates inconsistent
impact equations
Changing kinematic relationships .
Stack of bricks
Stack of bricks simulation
Contact regions
Simple simulation involving friction.			. .
Graphs for the simple friction simulation
PD controller used for positioning
Low-level controllers
Simulated biped model
Swing phase constraints. .
Walking cycle
Newton statistical output for walking simulation
Definitions for the walking algorithm
Walking algorithm initialization procedure
Abbreviated walking algorithm
Movement of a redundant robot arm
82
83
84
85
89
91
91
93
97
97
109
109
111
112
113
114
115
115
116
118
lx
Chapter 1
Introduction
Right from the time of their invention, computers have been used extensively for
the simulation and analysis of physical systems. The aerospace and automotive
industries, for example, have extended great efforts toward building simulation
and analysis tools to aid in their design and testing processes. The capabilities
of simulation programs have grown in proportion to the growth in available com-
puting power, but the basic techniques used have not really changed that much.
This thesis claims that computer science techniques can be used to build funda-
mentally different simulation systems, with much more general purpose analysis
capabilities.
Past and current simulation systems have been very well suited for analysis
of particular problems. They often use a sophisticated, special purpose model
of a particular class of physical systems that enables them to do a good job of
simulating the aspects of behavior accounted for by that model. N\Thlle certainly
2
useful, such systems still leave much to be desired. In particular, if, when using
such a system, a designer wishes to account for aspects of behavior not already
included in the model, there is often almost no simple way to do so. It would
clearly be desirable to be able to easily incorporate new aspects of behavior into
the model used by a simulation system.
Current computer-aided design systems do not really support most aspects
of design. Many are little more than geometric modeling systems, allowing for
detailed specification of object geometry and possibly some material properties
but providing very little in the way of tools for design analysis. More powerful
systems will allow building and testing of "electronic prototypes," sophisticated
object models incorporating information about object behavior from a variety
of domains. An electronic prototyping system will allow easy specification of
object designs, support a wide range of analyses of such designs, and enable use
of design analysis information in any necessary design modifications.
Consider, for example, the design of an anthropomorphic multi-fingered robot
gripper. Attempts to integrate a control algorithm and a particular hand design
could expose flaws in either the control algorithm or in the mechanical design.
A powerful design system would allow immediate redesign and testing of both
of these components. A designer could construct an electronic prototype of the
hand and experiment with several different configurations before comrnitting to a
specific design. Figure 1.1 shows designs modeled after the Salisbury hand. The
designer might first choose a three-fingered model but be unable to find control
algorithms that achieve design objectives. ?Vith minimal effort the designer could
3
Figure 1.1: Different hand configurations
test control algoritlims on a four-fingered model to determine whether such a
model would better meet specifications. Ease of redesign facilitates discovery of
an optimal match between control algorithms and mechanical design.
By demonstrating the effectiveness of the following ideas, this thesis argues
that more powefful simulation and design tools can indeed be built:
o+ we can significantly extend traditional simulation system capabili?ties by
fully integrating geometric modeling techniques into such systems. There
are a number of advantages to cleanly interfacing a complete geometric
modeling system to the analysis component of a simulator. In particular,
collision detection and handling is made a much more natural and man-
ageable feature of the system; no longer will ad hoc handling of object
interference be necessary.
o+ we can provide significant amounts of automatic aid to analysis. Ease of use
is critical to the success of design and analysis systems. A primary mech-
4
anism for achieving such ease of use is the capability of incrementally and
automatically modifying object behavior descriptions during a simulation.
In the domain of dynamics, for example, we show how systems can de-
tect changes in kinematic relationships between objects and incrementally
modify their internal description of the objects' behavior.
o+ we can use a rich, structured, hierarchical object representation as the basis
for sets of design and analysis tools. The representations encapsulate ob-
ject information from a wide variety of domains of interest; the information
is not organized solely in accordance with the needs of a specific type of
simulator or analysis tool, but instead in a more natural manner that al-
lows for general access to and reasoning about the object. Thus, structured
descriptions of objects can be used by various simulators simply by provid-
ing means for extracting the relevant information in the particular format
required. For example, a kinematics simulator needs ouly to extract the
inter-object kinematic constraints in order to carry out a simulation. On
the other hand, the same object model can be used for dynamics simulation
if the simulation module is able to extract from the object representations
the appropriate masses, volumes, and inertias, in addition to the kinematic
constraints. Path planning algorithms and assembly planners could be run
and their results then used by other modules. For instance, a path planner
could synthesize a rough plan, and then a control module of a dynamics
simulator could effect this plan in simulation. If problems with the plan
were discovered during simulation, the control module could reinvoke the
5
path planner to request a new, more detailed plan before continuing. By
sharing a common general object representation, communication among
various analysis modules is made simpler and more productive.
The thesis supports these claims in the context of the description of a general
purpose dynamics simulator. The system, called Newton1, is based on a hierar-
chical, m9del-based object representation scheme. It automatically formulates
and integrates object motion equations, uses geometric models to detect and
resolve impacts, uses a general-purpose event-handling mechanism to support
incremental modification of object behavior descriptions as object relationships
change (due, for instance, to temporary contacts being made or broken), and pro-
vides facilities for specification of control algorithms. Newton was implemented
during the development of the ideas outlined above and, while not incorporating
them completely, demonstrates their effectiveness. The system has been used in
various asp?cts of physical system modeling and simulation research at Cornell
University, Purdue University, and the University of Iowa.
1.1 Organization of the thesis
Chapter 1 outlines the goals and ideas presented in the thesis. Relevant related
work on physical system simulation is reviewed to clarify the contributions of the
thesis. The ideas and contributions of the thesis are then reviewed and contrasted
with the previous research.
1"Newton" is actually the name for a larger set of design, simulation, and geometric modeling
tools developed at Cornell and Purdue Universities. In this thesis, however, "Newton" will he
used primarily to refer to the dynamics simulation system.
6
Chapter 2 presents ideas concerning the structure of general purpose design
and simulation systems. A general simulation system architecture is proposed.
The architecture of the Newtonsystem is summanzed and an example simulation
presented.
Chapter 3 examines the issues of object representation and definition in a
general simulation system. The basic representation scheme, in which object
descriptions are organized as a collection of models from different ddmains, is
described. For the Newton system, descriptions of the kinds of objects that can
be modeled and the particular representations used are given. The structures of
the modeling systems underlying the object representation are also presented.
In Chapter 4, the analysis component of Newton is described. In particular,
the problems of automatic motion equation formulation and constraint main-
tenance are investigated. As simulations progress constraints on the relative
motion of objects change, due, for example, to temporary contacts. At all times
a consistent (usually independent) set of constraints must b? iiiaintained in ordei-
to integrate the system of dynamics equations.
Chapter 5 discusses the handling of exceptional events in simulation systems,
and argues that incorporation of general-purpose event-handling mechanisms can
provide critically needed automatic aid to experimenters. In particular, the New-
ton system's facilities for modeling impact, contact, and friction are presented.
Those aspects of dynamic analysis that still present difficult problems are also
detalled: tractable resolution of multiple contacts, suitable impact resolution
models, constraint maintenance, and so on.
7
Chapter 6 presents a detailed examination of Newton's control module. The
control module is becoming an increasingly important component of the system.
It is used to provide ways for users to influence the evolution of motion in a sim-
ulation through the specification of control programs. Some such programs are
simple low-level controllers for producing joint torques while more interesting in-
vestigations involve methods for high-level specification of complex control tasks
and development of algorithms for realizing such specifications. The problems of
integrating a control module into a simulation system are also discussed.
Chapter 7 reviews the goals and contributions of this work, lists a number
of interesting and difficult remaining problems, and discusses the focus of our
current and future work on physical system simulation.
1.2 Review of related research
This section examines related previous and current research on physical system
simulation. Its focus is on revealing the common structure of previous simulation
systems so that comparisons can be made with the contributions of this thesis.
An important early influence on our work was the observation that the large
amount of work directed toward development of computer aided design systems
was producing programs that were good for specification of geometry and for
making pretty pictures, but that productive use of these systems was limited by
the failure to incorporate sufficient design analysis tools. During the development
of the general purpose dynamics simulator, Newton, the primary related work
came from the mechanical engineering domain. More recently, the computer
8
graphics and animation community's interest in what they call physically based
modeling has influenced work on Newton.
1.2.1
Mechanical engineering work on physical system
simulation
Extensive mechanical engineering research has led to many developments in phys-
ical system simulation and to many implementations of simulators. The ADAMS
[15j and DADS [31? systems are examples of large state-of-the-art systems from
the mechanical engineering domain. In many ways such systems are very sophis-
ticated: efficient formulations of mechanism dynamics are supported, intricate
numerical techniques for solving equation systems are used, object flexibility and
elasticity are often handled, etc. However, from a computer science perspective
some things are lacking. Richer object representations are needed. Typically,
systems have airnost ignored geometric considerations and represented objects
simply as point masses with associated inertias and coordinate systems. Geo-
metric modeling techniques have matured enough to allow object representations
used by dynamic simulations to include a complete geometric description usable
by a sell-contained geometry processing module. Furthermore, impact, contact,
and friction are typically handled by current systems in an adhocor rudimentary
manner, if at all.
Surveys of mechanical engineering simulation systems can be found in [32,
56].
9
ADAMS
ADAMS[14,15,52jwas developed by N. Orlandea,M. Chace, and others at the
University of Michigan and related work continues commercially at Mechanical
Dynamics, Inc. ADAMS has been a successful, widely used system. It includes
sophisticated automatic formulation and integration of motion equations and the
use of efficient numerical analysis techniques. As is typical of most mechanical
engineering systems, object geometry is not fully represented. Shapes of rigid
components of the mechanisms may be described by wire frame models at the
user interface level, but internally the bodies are simple point masses having a
local coordinate system and user-supplied mass and inertia properties. Impact is
handled in a limited way using markers. A marker is a distingnished point in an
object's coordinate system. To model impacts, a user associates a pair of markers
from the two colliding objects, and associates with those markers a function
providing a repelling force. To model impact well, the function would have
negligible effect except at very small inter-marker distances. Such a scheme for
impact handling is not very satisfying, but is typical of many systems. Modeled
this way all possible collisions must be explicitly foreseen by the experimenter
before the simulation. In contrast to the way that Newton handles impact and
contact, the system's equations of motion don't change over time. As will be
discussed in detail later, it is advantageous for a system to be able to modify
the equations of motions as new kinematic relationships are established in the
course of a simulation. Especially in such tasks as robotics parts assembly, the
constantly changing contact relationships make a scheme like the use of markers
10
unappealing.
DADS
DADS is a similar large three-dimensional dynamic analysis system. It was
developed by E. 11aug and others at the University of Iowa's Center for Computer
Aided Design and also commercially by Computer Aided Design Software, Inc.
As with ADAMS, it has been used to model many sophisticated mechanisms. It
too, however, has a weak object representation scheme and its modeling of impact
is, for the most part, similar to that of ADAMS. A few more recent papers [35,36,
37] do discuss automatic reformulation of constraints as kinematic relationships
change, but the general event-handling issue is not examined in detail and it
appears that not much of this has yet been incorporated into DADS. DADS has
a notably sophisticated and integrated control module [64], though it does not
seem to address issues of high-level specification of control algorithms. Haug's
group has published extensively about their work, detailing current techniques in
many areas of simulation of mechanics [33,34] ([34] is a comprehensive University
of Iowa technical report containing many of Haug's group's papers for the period
1980--H1986).
Other systems
The system MAP, developed at General Electric, Inc., is notable for its inte-
gration of flexibility effects with rigid body dynamics analysis. Collisions and
contact are also handled in a relatively sophisticated manner, though geometric
modeling techniques again play a minimal role. It is a proprietary system, and
11
thus difficult to assess in detail.
J. Wittenburg's MESA VERDE system[74] is another quite capable mechani-
cal analysis system. It too efficiently formulates motion equations for both closed-
and open-loop systems, handles impacts and friction, and provides facilities for
mechanism control. And like the other systems, it lacks full representation of
object geometry, high level specification facilities and extensibility of coverage.
The products of SILMA, Inc., of Los Altos, California, and Deneb Robotics,
Inc., of Troy,Michigan, are systems representing the state of the art in commercial
CAD systems for robotics applications. The analysis capabilities of these systems
are, for the most part, kinematic, though SILMA has just recently released an
early version of a dynamics component for their system. Some of the features of
SILMA's system, CimStation, are described in a paper by Craig[17].
1.2.2 Work in the computer graphics and animation
fields
Computer graphics and animation researchers are becoming increasingly inter-
ested in what they term phys?caily-based modeling. Their desire to create increas-
ingly complex animations has made traditional keyframing and in-betweening
techniques less successful. In order to produce more realistic animation these
researchers have become very interested in modeling and simulating object dy-
namics. While simulation systems from mechanical engineering are generally
much more sophisticated in their dynamic analysis capabilities, a number of as-
pects of the animation-oriented systems are of interest. In particular, the focus
of many engineering systems is on simulation of highly constrained mechanical
12
systems (e.g. four-bar linkages). Animation-onented systems, on the other hand,
have concentrated more often on simulation and control of objects having many
degrees of freedom, particularly models of humans and animals.
Spacetime Constraints
Some of the most recent work of interest is that by Witkin and Kass [72]. Their
concept of spacet?me constraints has been demonstrated by impressive amma-
tions of a Luxo larrip. The idea of spacetime constraints is to combine an overall
physical constraint that objects obey Newton's laws with other, user-defined,
constraints that specify goal positions, velocities, etc. at different times. The
method works by discretizing the (time-dependent) constraints into related con-
straints for each instant of time. Then a constrained optimization technique is
used to solve the constraints and produce, all at once, values of the state van-
ables for each instant. The work emphasizes the desirability of being able to
express constraints in a high-level level manner. In the Luxo examples, it is
demonstrated that constraints such as "clear the hurdle," "don't waste energy,"
and "land softly" are simply and clearly expressible.
The spacetime constraints scheme might be called an off-line simulation method
since state variable values for all time instants are solved for concurrently; each
iteration of the optimization technique produces a complete new approximation
of the animation sequence that would optimize the constraints. For production
of certain kinds of realistic-looking animations the method is clearly useful and
powerful. However, for other animations, use of general purpose simulation sys-
13
tems would probably be more effective. For instance, the spacetime constraints
paradigm is not likely to easily incorporate certain impact models or the handling
of kinematic relationship changes for which the precise times of the changes are
not known beforehand.
Other graphics and animation work
Barzel and Barr's notion of dynamic constraints [10] is another approach being
used to produce better animation through modeling of dynamics. A general
approach to specification and formulation of constraints on object motion is pre-
sented. The formulation scheme is related to methods of constraint stabilization
and also PD controller specification from traditional mechanical simulation. A
benefit of this idea is that dynamic constraints solve "for free" the problems of
getting objects into desired initial positions; if intended kinematic and position-
ing constraints are not initially met, the stabilization properties of the derived
dynamic constraints cause the objects to move toward configurations where the
constraints are met. Handling of impact, however, is achieved via schemes like
the use of "force-fields" mentioned earlier, and the problem of user-directed con-
trol of high-degree-of-freedom mechanisms is not addressed in depth. Platt, from
the same group, applies related techniques to animate flexible bodies[58]
Numerous other simulation-related papers have appeared recently in the graph-
ics and animation literature. Some of the recent work has begun to address
issues involved in providing general-purpose simulation capabilities (e.g. deal-
ing with impact and contact problems) [30,49]. The most sophisticated work is
14
actually from the traditional CAD/simulation domain, but addresses the graph-
ics community[18J. Various other works [3,4,13,29,42,43,68,69,7o,71,75,76,77]
concentrate on particular aspects of dynamic and/or kinematic simulation and
related motion-control problems, while handling of flexible bodies is emphasized
by Terzopoulos[60,61,62,63].
As with the mechanical engineering systems, the animation-oriented systems
have been used to do very intersting simulations of complex objects. However,
overall, these systems also exhibit many of the characteristics that make the
engineering systems inadequate as truly general purpose computer aided design
and analysis tools. They lack powerful representational facilities; some produce
pretty pictures using complex geometric models but these models are used very
little in the actual analysis. They are often restricted in coverage; some deal
only with open-loop systems or are meant to deal just with human movement
simulation. Eriction, impact, and contact are handled, if at all, in simple or ad
hoc ways.
1.2.3 Other related work
Research from a number of other related areas has also influenced our work.
As mentioned, geometric modeling research is relevant since a large burden of
analysis is put on the Newton system's geometric modeling component. The
bibliography contains a number of geometric modeling references[40,44,50,55].
Work on positioning of articulated figures is another area of interest, since before
simulations can take place objects must be put into their initial positions. This
15
involves meeting both the structural constraints between the components of an
articulated figure and any user-chosen initial positioning constraints [6,9,45,63].
Other related works come from a number of areas, including modeling, constraint
systems, and simulation[1,2,8,12,38,39,54,65].
1.3 Summary
There is great interest in the development of general purpose computer aided
design systems. To date, however, success has been largely limited to the pro-
duction of tools that are suited only for particular aspects of design: computer
aided design systems have primarily emphasized specification of geometry; sim-
ulation systems from the mechanical engineering field have concentrated mainly
on formulation and integration of an unchanging set of equations describing ob-
ject behavior; and work by computer graphics and animation researchers has
often concentrated on producing good-looking animations without true regard
for whether the generated motions were physically realistic. Flexible systems
that integrate many aspects of design and analysis have not yet been built.
What we would like are systems supporting electronic protot?ping. These sys-
tems will both support specification of object design and provide a full range of
tools for analyzing and testing of designs. Using electronic prototyping systems,
the design process can be made more efficient. Costly physical prototypes will
no longer need to be built. Since after design specification the "building" of
electronic prototypes is essentially free, designers do not have to worry about
breaking or damaging expensive prototypes and will be able to carry out a wider
16
variety of analyses. It should also be feasible to experiment with a wider range
of designs since minor design variations should be easily specified and tested. In
other situations it is not even feasible to build a physical prototype because the
environment in which the product will be used (e.g. space) cannot be satisfac-
torily recreated in the lab.
This thesis, then, argues that simulation and analysis systems can be made
significantly more powerful. In particular, by using geometric modeling tech-
niques to include a full description of object geometry, previously difficult-to-
incorporate simulation system features, such as collision detection and resolu-
tion, can be handled much more naturally. In addition, design can be made
easier by providing more automatic aids to analysis. One primary candidate for
such automation is the incremental modification of object behavior descriptions
during simulations, and, in particular, the thesis proposes and demonstrates the
use of a general purpose event handling mechanism to stipport automatic han-
dling of collisions, changing contact relationships, and other events that can cause
discontinnities in object motions. Finally, it is argued that structured object rep-
resentations that model object information from a variety of domains can provide
the basis for more flexible simulation environments. These environments allow
addition and modification of analysis tools and, by supporting communication
between the different tools, increase their overall usefulness.
Chapter 2
General Purpose Simulation and
Analysis Systems: Architectural
Issues
In this chapter, we elaborate on the ideas outlined in Chapter 1 to propose basic
architectural features for general purpose simulation and analysis systems. Then,
the overall structure of the Newton system is described, revealing the partial
realization of these ideas in a particular implementation.
2.1 Simulation architectures
The shortcomings of systems intended to be used for computer aided design and
analysis have been that they generally are suited only for particular aspects of
design. Using some systems, geometric specifications can be made, after which
17
18
a limited set of analyses may be run on the resulting representations: integral
property calculations, such as volumes, surface areas, and inertias; finite element
mesh generations for stress analysis, swept volume computations for use in an-
alyzing tool cutting paths; and so on. Simulation systems, on the other hand,
can formulate and integrate a complex set of equations describing the behavior
of a set of objects. However, use of geometric information for collision detection
and analysis of contact relationships, and automatic modification of object be-
havior descriptions (as encoded in the set of dynamics equations) as inter-object
relationships change are not usually supported.
These problems can be alleviated by the use of richer object representations.
Object representations should contain more thorough coverage of the many fea-
tures aifecting object behavior and more complete descriptions of the features
included. Of course, geometry should be a primary component of such a repre-
sentation. Along with geometry, however, should be information about how an
object is related to other objects (e.g. what kinematic relationships directly in-
volve this object or what component objects are contained by this object), what
its material properties are, how it should be drawn on a screen, etc.
By recording the object information in an organized manner, it will be possi-
ble to build a variety of simulation tools that are based on and that can commu-
nicate through this common representation. By providing appropriate interfaces
between the representation and analysis tools, the analysis tools can be made
relatively independent of the specific details of the particular implementation of
the object models.
19
Thus, object-oriented programming techniques seem well suited to serve as
the basis of object representation and as the basis for communication between
various analysis modules. Object information is then divided by domain into a
set of models. Each model type has an associated modeler which is responsible
for handling requests for object information for its particular domain. Again,
interfaces are used between analysis systems and modelers to separate the actual
implementation of the model information from the rest of the system. This allows
for easy modification of model organization or modeler algorithms.
The use of a rich object representation will only be worthwhile if it can indeed
support a variety of simulation and analysis tools. The clean separation of knowl-
edge and explicit representation of object relationships makes reasoning about
object behavior feasible. Kinematic or dynamic analysis tools can examine the
interconnections between objects in order to formulate equations describing ob-
ject motion. Explicit inclusion of object geometry allows for interference checks.
Material property information allows reasoning about friction effects. The 5up-
port for reasoning about various aspects of object behavior makes creation of a
variety of general purpose simulation tools possible. Cooperation among these
tools is possible because they can communicate through the common representa-
tion, e.g. by recording results of various analyses with the object for other tools
to use.
20
2.2 The Newton system
Using Newton, a designer can define complex physical objects and mechanisms
and can represent object characteristics from a wide range of domains. An ob-
ject is made up of a number of "models," each responsible for the organization of
obje?t characteristics from a particular domain. In most simulations the basic do-
mains of geometry, dynamics, and controlled behavior are modeled. A dynamic
modeling system, for example, is responsible for maintaining an object's posi
tion, velocity, and acceleration, and for automatically formulating the object's
dynamics equations of motion. A geometric modeling system is responsible for
information about an object's shape, distingnished features on the object, and
computation of geometric integral properties such as volume and moments of in-
ertia. It also detects and analyzes object interpenetrations so that an interference
modeling system can deal with collisions between objects. With this kind of flexi-
bility, mechanism design and analysis is made simpler; a number of simulations of
a physical system might be carried out, with different sets of modeled properties
being accounted for each time. New modeled aspects might be added to increase
the overall accuracy of the simulation, or certain domains might be abstracted
or ignored to allow the experimenter to focus on the contribution of other do-
mains to the observed behavior. In addition, since the details of each modeling
system are hidden from the rest of the system by an interface, modeling systems
may be replaced easily just by rewriting the modeler interface. Thus, a simple
but fast geometric modeler might be used for rough early experiments; a more
21
sophisticated modeler could then be used in later stages of design refinement.
To help clarify some of the issues that will be discussed in this chapter, a
simple simulation example is presented in the next section. The later sections
then present the actual description of the system architecture.
2.3 An example simulation
In Figure 2.1 fifteen frames from a simple simulation are shown. The simulation
exhibits several 0 the novel analysis features supported by Newton, including
use of the geometric modeling subsystem and the event handling w?chanism to
detect and handle collisions and to incrementally modify object descriptions in
accordance with changes in kinematic relationships.
Object and Scene Definitions
Figure 2.2 shows the three class definitions needed to define the objects in the
scene. Two primitive object types are defined: a cuboid-shaped object t?pe
named link and a sphere-shaped object type named ball. A composite type,
pendulum, is defined, having four links, a ball, and four ball-and-socket hinges
(spherical joints). The hinge attachment points were placed off-center in order
to produce spinning and jangling motions in the links. The description of the
initial scene, a configuration structure called test, is shown in Figure 2.3. It
specifies the objects that make up the scene along with their initial positions and
velocities.1
`Note the explicit specification of initial object locations using three numbers for each object's
position and four more numbers (Euler parameters) for each object's orientation. It is clearly
undesirable to require users to numerically specify the exact initial positions and velocities of all
22
-1 flu I a
Th
4/
Figure 2.1: An example simulation
23
primitive link (thickness, height, rho)
begin
properties: (density: rho)
geometry: cuboid (thickness, height, thickness) where begin
topback: (0, height/2, -thickness/2 );
botback: (0, -height/2,-thickness/2 );
botoffmid: (0, -(height/2) - (.75 * thickness), 0)
end
dynamics:;
interference:;
control:;
end
primitive ball (radius, rho)
begin
properties: (density: rho)
geometry: sphere (radius);
dynamics:;
interference:;
end
composite pendulum (thickness, linkheight, rho)
components
10 : link (thickness, thickness, rho);
11,12,13 : link (thickness, linkheight, rho);
ball: ball (.75 * thickness, rho);
jI,j2,j3,j4 : ball-and-socket
structure
join 10
join 11
join 12
join 13
to
to
to
to
begin
dynamics:;
control:;
end
11 with ji matching (botback topback);
12 with j2 matching (botback topback);
13 with j3 matching (botback topback);
ball with j4 matching (botoffmid center)
Figure 2.2: Object class definitions for the example simulation
24
configuration test
components
block : link (1.5, 7. 1);
ground : link (20, .5, 1);
pendulum : pendulum (1 3, 1)
state
pendulum.10:
pendulum.11:
pendulum.12:
pendulum.13:
pendulum.ball:
block:
ground:
((			0, 15, 0)			(1, 0, 0, 0));
((			1.5,			14.5,			0) (.70710677, 0, 0,
((			4.5,			14.5,			0) (.70710677, 0, 0,
((			7.5,			14.5,			0) (.70710677, 0, 0,
((			9.75			, 14.5, 0) (.70710677, 0, 0,
((--H2, 2.75, 0) (1.0, 0, 0, 0));
((0, --H1, 0) (1, 0, 0, 0))
.70710677));
.70710677));
.70710677));
.70710677));
Figure 2.3: Scene specification for the example simulation
The simulation
The simulation begins with the objects in the initial positions shown in the first
frame in Figure 2.1. Initial transiational and rotational velocities of all objects
are zero and two objects, the ground and the small top cuboid of the pendulum,
pendulum. 10, are controlled to remain fixed.
Before entering the main simulation loop, Newton notices that the tall block
and the ground are in contact but have no relative velocity in the direction of
the ground's surface normal. So, as will be discussed in Chapter 5, the system
puts a temporary hinge between the ground and the block to model this contact
relationship. A new composite object is also created to model the pair of objects
objects. Some small initial positioning packages have been written tor Newlon but the problem,
in general, is a difficult one. Chapter 7 briefly addresses this problem.
25
and their constrained relationship. As simulation then progresses, the ground
and block remain stationary while the pendulum swings down. Impact between
the pendulum's ball and the block is detected by the geometry system after
about 1.5 seconds (sixth frame of figure).2 The system formulates and solves a
set of equations to deal with this collision, and in this case uses a coefficient of
restitution of .7 between the ball and the block.
After the collision the contact relationship between the block and the ground
changes, since the impact causes the block to gain translational and rotational
velocity components. The contact between block and ground changes from face-
face to edge-face contact and the temporary hinge representing the contact re-
lationship is modified accordingly. As the block continues tipping the contact
relationship continues to change. Since the block is also rotating around its ver-
tical axis, the contact changes from edge-face to vertex-face and then back to a
different edge-face relationship. Later, when the top corner of the block hits the
ground (block-ground collisions were treated as totally inelastic in this simula-
tion), the contact becomes edge-face again. Eventually, after a few more changes,
the block slides off the ground altogether, causing removal of the temporary hinge
between the objects and deletion of the composite object containing the hinge
and the objects.
2In this example, the system was directed not to detect collisions between components of the
pendulum. For example, adjoining links in the pendulum were allowed to pass through each
other.
26
2.4 Newton architecture overview
Newton has three main components: the definition and representation module,
the analysis module, and the report system. The definition module is responsible
for analyzing high-level language descriptions of Newton entities, for organizing
information in the appropriate data structures, and for providing the interfaces
to the various modeling subsystems. The analysis component implements the
top-level control loop of simulations and the event handling mechanism, and
coordinates the working of various analysis subsystems. The report system han-
dles generation of graphical and statistical feedback to users during simulations
as well as recording of relevant information for later regeneration of animations.
The basic organization of Newton is depicted in Figure 2.4.
2.4.1 Definition and representation module
Simulated objects and simulation scenarios are described in a high-level definition
language. The language is then compiled into appropriate data structures, after
w?? sim?ations can ?e run.
Newton's definition language is used to describe a variety of simulation en-
tities, including objects, hinges, constraints, models, equations, and quantities.
Simulated objects are represented as instances of the object class. This class is
further divided into two subclasses: primitives, corresponding to single indivisible
bodies, and composites, representing collections of objects related by constraints.
A primitive object is single indivisible body. The constrained relationships usu-
ally correspond to material hinges such as ball-and-socket or pin joints and are
27
Definition Module			Analysis Module			Report Module
defluition
language
abstract
modeler
interference
modeler
handler
Impact
II
contact
control
others
other domains
simulation
top level
equation
formulation
equation
integration
equation
solving
clock
statistics
equation
manipulation,
simplification
Figure 2.4: Basic Newton architecture diagram
28
represented by instances of the hinge class. The components of composite objects
can be either composites or primitives. Thus, description of complex mechanisms
is made simpler by breaking the description down into natural part-component
relationships.
Quantities are data structures used to represent complex values relevant to
simulations. Such values include positions, velocities accelerations, masses, vol-
umes, coefficients of restitution and friction, forces and torques, and so on. Mod-
els are data structures used to partition information about an object according
to particular domains. Thus, the model class is subdivided into subclasses of ge-
ometric models, dynamic models, interference models, control models, etc. Each
model type is responsible for representing object information for its particular
domain. Chapter 3 presents the details of these representations.
One advantage of the hierarchical object representation scheme is that it
facilitates automatic incremental reformulation of an object's motion equations.
During simulation, Newton's analysis module makes requests for an object's set of
motion equations. The set is constructed recursively, by requesting the equation
sets for each of the object's components and for each of its hinges. At every
level of an object's composition hierarchy, the set of equations for that level,
once derived, is stored in the appropriate dynamic model. Then, when events
occur that dictate a change in the equations for a component of an object, only
the equations for that component need be rederived. The other components'
equation sets are still available from their dynamic models. This issue will be
discussed in more detail in Chapters 3 and 4.
29
(until (stop-simulation-p)
(save-current-state)
(integrate-from-current-state)
(handle-exceptional-events)
(report))
Figure 2.5: Newton's top-level loop
2.4.2 Analysis module
Newton's analysis module is responsible for overall coordination of a simulation.
After a simulation scenario has been defined using the definition module, the
main simulation loop executes according to the code in Figure 2.5.
At the beginning of each iteration, the current state of all objects is saved in
case an unacceptable attempt at stepping forward in time by some At occurs and
necessitates trying again with a smaller time increment. Next, the integration
subsystem is invoked, with the purpose of producing new positions and velocities
for all objects from their old positions, velocities, and accelerations. Since object
accelerations are needed during the integration process, the system collects and
processes the set of motion equations representing the current dynamic state and
then sets up and solves a linear system to produce the accelerations and hinge
reaction forces and torques.3 Once tentative new positions and velocities have
been computed, the call to handle-exceptional-events passes control to the
event handling subsystem. The event handler is responsible for detecting and
handling collisions, contact formations and breakages, control algorithm state
3Actually, control forces and torques can also be solved for, since New?on supports a mixture
of forward and inverse dynamics. This is described in the chapter on control.
30
changes, and other events that yield discontinuities in object velocities or accel-
erations. These events can invalidate the newly proposed state and necessitate
restoration of the previous state, integration using a different time increment,
resolution of collisions or other events, and so on. The event handler is described
in detail in Chapter 5.
2.4.3 Report module
As stated above, the report system is responsible for generating output that can
be of use to the experimenter in analyzing simulations. This includes display
of the simulation scene as execution takes place, and display of values of any
quantities that the user wishes to monitor. Typically, users choose to display the
evolution of values of a number of variables using a set of graphs. In Figure 2.6,
four frames from a simulation of a redundant robot arm ar? shown. The graphs
exhibit values of the position and acceleration of the arm's end effector, and joint
angles and accelerations for the distal three joints of the arm. The report package
is also responsible for the recording of information that allows later redisplay of
the simulation as a real-time `?
31
Y14M???i?_? _
___
. .-
-.
:.:.t?? :":M
. .-
1:___
.-
fflf 1 Hi ??
Figure 2.6: Redundant robot arm			simulation and graphs
Chapter 3
Representation and Definition
The following sections describe the definition component of the Newton system,
concentrating on the object representation scheme and the way in which it sup-
ports novel aspects of mechanical analysis. Using Newton's definition language
and interfaces, users can describe and manipulate a number of simulation enti-
ties. Foremost among these are objects, hinges, models, quantities, constraints,
and equations.
3.1 Objects and hinges
Instances of the object class are used to represent physical objects. In this thesis,
objects consist primarily of collections of rigid bodies related by constraints.1
Objects can have other objects as components, allowing complex mechanisms to
be described naturally and abstractly. Thus, the object class is divided into two
`There has been some preliminary work on incorporating analysis of nonrigid bodies into
Nett,ton.
32
33
subclasses, primitive and composite.
Definition 1 A primitive object is an indivisible physical body, not composed of
other objects.
Definition 2 A composite object is a set of objects, called the components of the
composite, interrelated by a set of constraints. The components of a composite may
be either primitive objects or other composite objects. The graph induced by the set
of components and constraints must be connected. It may be cyclic.
By these definitions, primitives correspond to the basic material bodies of
simulations, while composites provide for abstraction based on constrained rela-
tionships between components.
Composite objects usually correspond to logical units of interest to a user
of the system. The hierarchical organization allows easy definition of initial
simulation models, then refinement of components later on. For instance, an
ARM composite is often defined as being composed of two primitive objects, an
UPPER-ARM and a LOWER-ARM, and another object, a HAND. For initial
testing of general arm movement algorithms, the hand might be defined as a
simple primitive object. Subsequent experimentation with control algorithms
for fine manipulation tasks would substitute a fully developed composite hand
model. The hierarchical organization also plays a role in supporting incremental
modification of the equation sets representing object behavior.
Composites do not always correspond to user defined units. As contacts are
made and broken during simulations, the system creates and destroys composite
34
Object class
name
id
models
containing-composite
Primitive class
rncident-hingelist
Composite class
components
hinges
component-hingelinked-list
user-deflned-p
Hinge class
name
id
models
containing-composite
object-feature-pairs
position-constraints
velocity-constraints
acceleration-constraints
Figure 3.1: Basic structure of object and hinge classes
objects. Primitives, on the other hand, persist throughout a simulation.
The hierarchy induced by the part-component relationship must be strict.
That is, an object can be an immediate component of only one other object.
Hinges impose further restrictions on the structure of object hierarchies.
Definition 3 A hinge represents a constrained relationship between two objects. A
hinge constraint may be kinematic or nonkinematic, holonomic or nonholonomic.
Commonly defined hinges often represent standard material hinges such as
bali-and-socket, pin, and universal joints. However, hinges are also used to model
nonholonomic relationships such as surfaces rolling or sliding on each other, and
nonkinematic relationships such as gravitational attraction between bodies.
The structures of the object and hinge classes are shown in Figure 3.1.
35
We associate two kinds of graphs with composite objects. The first, called the
composition graph, represents the user-onented view of part-component structure
of an object. This graph is actually a tree. The root of the tree represents a
composite object viewed at the highest level of abstraction. The children of any
node are just the component objects of that node. Thus, all interior nodes of the
tree represent composite objects and all leaves represent primitives. The leaves,
clearly, represent all of the primitive components of the root object. Figure 3.2
shows the composition graph for a simple model of an arm.
The second kind of graph is called a strnctnre graph. For any composite
object, there is an associated set of structure graphs. The number of graphs
in the set is one greater than the depth of the object's composition graph. An
object's depth 0 structure graph simply consists of a single node for that object.
The depth 1 structure graph has the object's components as nodes and its hinges
as edges. A depth n+1 structure graph is then constructed from a depth n graph
by replacing all nodes in the depth n graph by the depth 1 structure graphs of
those nodes We distinguish an object's largest structure graph, containing all
of the primitive objects and all of the connecting hinges, and call it the pmmitive
st'ncture graph. Figure 3.3 shows the structure graphs for the simple arm model.
As mentioned earlier, one advantage of the hierarchical object representation
scheme is that it provides support for incremental reformulation of object motion
equations. During simulation, Newton's analysis module makes requests for an
object's motion equations. If the object is a composite, the set is constructed
recursively, by requesting the equation sets for each of the object's components
36
arm
base
thumb			fi			f2			f3
Figure 3.2: Compositiou graph for sin?pIe arm mo?Ie?I
upper arm			lower arm			hand
0
arm
thumb
fi
f2
upper arm			lower arm			base
upper arm
thumb.ji			thumb.j2			thumb.j3
f4.ji			f4.j2			f4.j3
Figure 3.3: Structure graphs for simple arm model
37
and for each of its hinges. At every interior node of the composition graph, the
complete equation set for the corresponding composite object, once derived, is
stored in that object's dynamic modd. Then, when events occur that dictate
a change in the equations for a component of an object, only the equations for
that component will need to be rederived. The other components' equation sets
remaan available from their dynamic models.
3.2 Models
Every object contains a set of models. A model for an object organizes informa-
tion about the object for a particular domain. The primary domains currently
are geometry, dynamics, interference, control, and abstract.
Associated with each modeling domain is a modeter or modeting system. A
modeler is responsible for performing calculations about its particular domain.
Requests for object information are handled by interfaces that direct the requests,
along with the relevant object models, to particular modeling systems. In ?his
way the specific details of the model and modeler implementation are hidden
and can be modified easily. For instance, we can use any geometric modeling
system to handle the geometric calculations in simulations so long as it supports
a certain set of queries. An important issue involves determining what operations
a model or modeler must support.
38
Abstract model
All object's abstract model is responsible for abstract descriptive information
about the object. In particular, the abstract model contains the information
necessary for constructing the object's composition and structure graphs. Thus,
it records the object's component objects and the hinges that relate the compo-
nents. In addition, other object properties are maintained in abstract models.
Included are things such as the name of the material the object is made of, the ob-
ject's color, its density, etc. The distribution of object properties among model
types is not rigid. Some properties that are currently maintained in abstract
models belong naturally in object models not yet covered by the system. For
instance, if future applications depended more heavily on factors related to the
material composition of an object, a material model could be added to record the
material properties of interest. Through an associated material modeling sub-
system, this could include computation of coefficients of friction and restitution
from the material information, and so forth.
Geometric model
This thesis has stressed the importance of integrating geometric modeling tech-
niques with physical simulation techniques. Thus, simulated objects contain
geometric models that are manipulated by Newton's geometric modeling system
to handle any geometric problems that arise.
In order to carry out simulations, Newton's dynamic analysis module needs
access to a variety of geometric information. An important aspect of the design
39
of a geometric modeling system is deciding what operations the modeler needs
to support. For simulation of dynamics these operations include:
o+ calculation of volume and volumetric moments of inertia.
Calculations of other integral properties, such as surface area, are needed
for sophisticated handling of contact and friction problems and should be
supported as well.
o+ intersection of geometric models.
Testing of object interpenetration is necessary for handling of impact and
contact. Thus, efficient geometric model intersection tests must be sup-
ported. For dealing with contact problems, it is, of course, not enough to
simply report that the objects in contact intersect. The modeling system
must be able to analyze tile nature of the contact (e.g. surface-surface,
face-edge, vertex-edge, etc.) and provide a useful description of the region
of contact. The inclusion of space-time swept-volume operations would aid
in the efficient calculation of impact occurrence times.2
o+ identification of geometric entities (e.g. points, line segments, surface re-
gions) as distinguished features.
It is useful to allow users to reason about the objects at an abstract level.
Thus, hinge specifications in the definition language relate names for ge-
2Full support for boolean operations on geometric models is not normally required. Tests
for interpenetration are sufficient for collision detection. For contact analysis, the continuous
contact regions are at most tw??dimensional, so general three-dimensional boolean operations
are not needed. however, since they provide a mechanism for defining complex primitive object
geometries, they should be included in geometric modelers used by general purpose simulation
systems.
40
ometric features of objects rather than explicit mathematical descriptions
of the features.
Different geometric modeling systems using different geometric representation
schemes can be used by Newton as long as they support the necessary operations.
The simulation examples shown in the thesis use a simple geometric modeler
supporting only a few geometric shapes cuboids (parallelpipeds), spheres,
cylinders, cones, and pipes (cylinders with cylindrical holes along the main axis)
__ and the bare minimum of operations. To replace this modeler by another,
more sophisticated modeling system, it is only necessary to modify the small set
of interface functions that translates simulator requests into appropriate modeler
requests and interprets modeler replies for the simulator. One of the reasons for
not currently using more sophisticated modelers in Newton is that geometric
modeling systems supporting the required operations are not readily available;
most systems are either too slow or are not sufficiently robust.
Dynamic model
An object's dynamic model contains the bulk of the information crucial to carry-
ing out the numerical integration component of simulation. Symbolic representa-
tions of the object's equations of motion are stored in this model. Niore precisely,
motion equation schemata, the actual motion equations, and any compiled mo-
tion equations are included. Equation schemata and equations are discussed
in Section 3.4 below. Dynamic models are also responsible for maintaining an
object's state vector. The state vector consists of an object's position, orienta-
41
tion (in Euler parameters), velocity, angular velocity, acceleration, and angular
acceleration, all in global coordinates.
Interference model
The purpose of the interference model is to aid in the detection and re?solution of
impacts and contacts. Toward this, information is maintained about how object
intersection tests should be carried out. A hierarchy of simple geometric models
enveloping the object can be used to make the intersection tests more efficient,
for example.
In addition, equation schemata, and corresponding equations, modeling the
response of the object to collisions are formulated and maintained in the inter-
ference model.
Control model
Control models are used to maintain control force and torque quantities that
influence the motion of objects, and also to record equations used by control
algorithms. Control is discussed in detail in Chapter 6.
3.3 Quantities
Quantities are an important class of data structures representing primitive typed
values relevant to objects and simulations. Values that should be referenceable
through their type and name information and that need to be accessed, modified,
or computed in such a manner are represented using quantity structures. Quan-
tities are organized by type into classes that share behavior. Example classes
42
Quantity instance
value
owner
timestamp
Quantity class
name
valuetype
valuestateproperties
ancestors
descendants
owner-type
get-method
set-method
computemethod
Some Built-in Newton Quantity Types
position
orientation
velocity			(translational and
_____________________ acceleration angular)
mass
density
volume
moment-of-inertia
constraint force
control force			(and torque)
Figure 3.4: Quantities
include mass, volume, moment of inertia, as well as various transformation matri-
ces, forces and torques, and so forth. Thus, the quantity class hinge-control-torque
includes all quantities representing torques applied by controlled hinges. Names
are used to distinguish between instances belonging to one particular hinge, since
a hinge can apply as many torques as it has controllable axes. Figure 3.4 shows
the structure of quantity classes and instances and lists some of Newton's built-in
quantity types.
Three generic access functions, get-prop, set-prop, and compute-prop, are
used to read, modify, and compute quantity values. Information about the spe-
cific implementation of the operations is contained in the quantity's class. As-
sociated with the class are three methods: a get-method, a set-method, and a
compute-method. Thus, when a quantity is accessed through one of the generic
functions, the referenced quantity's class information is examined and the ap-
43
propriate access method invoked. Quantity values are computed on demand and
timestamped so that they can be recomputed automatically after they become
invalid (due to changing of values of quantities they depend on).
This organization frees system code and users from knowing the details of the
implementation of quantities. Object mass, for instance, is currently stored in the
object's abstract model. The mass quantity's access methods, invoked through
use of the generic access functions, direct requests to the abstract modeler inter-
face. The abstract modeler interface functions, in turn, may make new quantity
requests through the generic functions in order to process the original ones. For
instance, if the value of the mass was requested and the abstract modeler in-
terface found that no value was immediately available, it would make requests
for the values of the object's volume and density in order to compute the mass.
The volume and density requests are again made through the generic functions
so that the interface doesn't have to know their particular implementations. In
the implementation, the density request is itself handled by the abstract modeler
while the volume request gets directed to the geometric modeler interface. If,
as suggested earlier, a material model was added to the object representation,
object density would be maintained in the material model. Ouly the quantity
class methods for density would need to be modified to reflect this change.
3.4 Equations
Equations, clearly, are an important component of a simulation system since they
are used to describe mathematically the behavior of the simulated objects. In
44
Newton we commonly refer to three classes of equations: motion equations, con-
straint equations, and control equations. The distinctions are not strong. Gener-
ally, motion equations are those equations corresponding to Newton-Euler equa-
tions of motion relating object accelerations to forces and torques. Constraint
equations encode constrained relationships between objects and usually model
material hinges or other kinematic relationships between the objects. Control
equations represent any additional user-imposed constraints on object motions
(e.g. corresponding to inverse dynamic control as discussed in Chapter 6). Equa-
tions of all three types are represented in Newton using the equation schema and
equation data structures.
Equation schemata are structures indicating the algebraic relationship be-
tween various types of physical quantities. Equations, on the other hand, are
data structures that are derived from the equation schemata and reflect a cer-
tain view that the simulator takes of the relationships; equations correspond to
schemata whose basic parts have been bound to particular quantity instances
rather than quantity classes (see Figure 3.5). In addition, certain of the terms in
equations are viewed as knowns, having as value, for instance, the current value
of the corresponding quantity. Other terms are considered to be unknowns and
will have their values computed when the dynamic analysis module solves the
equations. In some circumstances, a third kind of equation data structure, corn-
piled equation, is employed. Compiled equations encode efficiently manipulable
versions of equations.
45
equation schema-			equation -
F			=			MA
external-forcequantityclass			external-force-i
external-force-? ??
acceleration-i
mass quantity class _
accelerationquantityclass			?			?mass-2
Object i
Figure 3.5: Equation schemata and equations
Chapter 4
Newton's Analysis Module:
Dynamics
The details of the basic elements of Newton's analysis module are presented in
this chapter. In particular, the issues of automatic motion equation formulation
and integration are discussed, leaving the problems of equation modification and
event handling for Chapter 5.
4.1 Introduction
As mentioned in Chapter 1, there has been a great deal of work done on methods
of dynamic analysis, especially in the mechanical engineering community. In the
past, formulation of equations of motion for a mechanical system was often done
by hand and encoded directly into a program. This is still done today when
extreme efficiency is required. Automatic techniques (e.g. as used in DADS and
46
47
ADAMS) have become common. Some of the techniques are quite sophisticated,
allowing very efficient simulation of certain classes of mechanisms. For the most
part, however, systems have not supported general automatic formulations that
allow efficient reformulation of equations as object relationships change during
simulations. Newton's object and equation representation scheme and the sys-
tem's modular design allow it to support a variety of dynamics formulations
without requiring the user to be concerned with its details. In addition, the ba-
sic formulation scheme used by Newton is powerful, while intuitively simple. It is
general purpose, applicable to mechanisms with both open and closed kinematic
chains and able to handle a large variety of inter-object constraints, and it sup-
ports incremental modification of the set of equations as relationships between
objects change.
4.2 Basic dynamic analysis
During a simulation, Newton's analysis module handles the top-level flow of con-
trol. It is responsible primarily for formulating and integrating the set of equa-
tions that represents the dynamic behavior of the simulated objects' and for co-
ordinating the detection and handling of events (called exceptional events) that
cause discontinuities in object velocities and accelerations.
After a simulation scenario has been defined using the definition module, the
simulation proceeds as shown in the code in Figure 4.1.
Simulation scenario specifications generally specify only initial positions and
1These equations will be called, a bit loosely, "motion equations." More precisely the equations
consist of differential equations of motion and algebraic equations of constraint.
48
(defun simulation-top-level ()
(establish-complete-consistent-initial-state)
(until (stop-simulation-p) do
(save-current-state)
(integrate-from-current-state)
(handle-exceptional-events)
(report)))
Figure 4.1: Newton'?s top level
velocities for objects and do not include descriptions of any initial inter-object
contact relationships. Thus, estabiish-complete-consistent-initial-State
is responsible for establishing a complete and consistent initial state before en-
tering the main simulation loop. It causes an initial formulation and solution of
the motion equations in order to compute object accelerations. In doing so, it
uses contact information provided by a call to the event handler to ensure that
the proper hinge constraints are derived.
The main simulation loop is then executed, each iteration beginning with
a complete state in which all object positions, velocities and accelerations are
known. First, the current state of all objects is saved in case an unacceptable
attempt at stepping forward in time by some occurs and necessitates trying
again with a smaller time increment. Next, the integration subsystem is invoked
in order to produce new positions and velocities for all objects based on their
old positions, velocities, and accelerations. Since object accelerations are needed
during the integration process, the system collects and processes the set of motion
49
equations representing the current dynamic state and then sets up and solves a
linear system to produce the accelerations and hinge reaction forces and torques.2
The integration routine is "atomic" in the sense that it assumes that object mo-
tion is continuous over the given interval. Once tentative new positions and
velocities have been computed, the call to handle-exceptional-events passes
control to the event handling subsystem, which determines whether or not the
proposed state is satisfactory. The event handler is responsible for detecting and
handling collisions, contact formations and breakages, control algorithm state
changes, and other events that yield discontinuities in object velocities or ac-
celerations. These events can invalidate the newly proposed state and cause the
event handler to restore the previous state, modify the equation system, reinvoke
the integration subsystem using a different time increment, resolve collisions or
other events, and so on. The event handler is described in detail in Chapter 5.
4.3 Motion equation formulation
The basic scheme for formulation of motion equations is simple. Its primary
advantage lies in its support for efficient reformulation of equations as inter-
object relationships change during simulation.
We begin by associating Newton-Euler equations of motion with each prim-
itive object. These equations are stored and maintained in each primitive's dy-
2Actually, other user-defined unknowns, and control forces and torques, can also be solve for.
See, for example, the discussion of mixed forward and inverse dynamic control in Chapter 6.
50
namic model. At the time an object is created the equations are of the form
mr =0
Jw+w x Jw --H 0
where m is the mass, i: is the acceleration of the center of mass, J is the 3-by-
3 inertia matrix, and w and w are the angular velocity and acceleration (with
respect to a coordinate system located at the center of mass but oriented like the
global coordinate system) of the object.
Hinges between objects impose constraints on object motion. Equations rep-
resenting these constraints are stored with the hinge. ?Vhen a hinge is created
between two objects, the constraint is analyzed to assess its effect on the objects.
Generally, for holonomic constraints, the second derivative of the constraint will
be used in conjunction with the basic motion equations to solve for object ac-
celerations and constraint forces. For example, a specification that two objects
are to be connected with a spherical hinge is realized by the inclusion of one
new vectorial constraint equation and the addition of some terms, representing
constraint forces, to the motion equations of the constrained objects. Thus, the
constraint representing the 3-degree-of-freedom spherical hinge is
ri t ci = r2 + C2,
where Cj is the vector from object i's center of mass to the location of the hinge.
Taking the second derivative of this constraint, we get
r4 + wi x c? + Wi X (wi x ci) = r2 + W2 X C2 + W2 X (w2 x c2).
51
The original motion equations for objects 1 and 2 represent 12 scalar equations
in 12 unknowns. The addition of the acceleration equation for the constraint,
however, results in 15 equations in 12 unknowns. In order to realize the effect of
the constraint, some terms need to be added to the original motion equations.
These terms represent the reaction force that will be produced in the hinge
to ensure that the conditions of the constraint are met. Thus, adding equal
and opposite constraint force terms, fhinge, to the original motion equations, a
determinate set of 15 equations in 15 unknowns results,
miri = fhinge
J1w1 + Wi X J1w1 = ci x fhinge
rn2r2 = --Hfhinge
J?w2 + W2 X J2w2 = c? X fhinye
?i + Wi X C1 + Wi X (wi xci) = r2 + W2 x c2 + W2 X (w2 x
Given positions and velocities, this set of equations can be solved as a linear
system at each simulation time step to produce object accelerations and hinge
constraint forces. As will be discussed later on, the constraint forces could be
eliminated from the equation formulation. However, in many cases the constraint
force values are of direct interest to experimenters; programs can monitor these
forces and model limited strength hinges that break when specified force or torque
values are exceeded.
Figures 4?2 and 4.3 show the constraint equations, their second derivatives
and the corresponding constraint term variables for the one-degree?ofreedom
52
revolute (pin) and two-degree-of-freedom prismatic hinges. Note that the pris-
matic joint's second degree of freedom allows rotation about the axis along which
translation is allowed.
Nonkinematic constraints can be handled by the same method. For instance,
a hinge representing a simple spring imposes the constraint
fhinge = k(c1 --H c2).
In addition, the object motion equations are modified to become
mir?i = fhinge
J1w1 + wi x J1w1 = cl X fhin9e
m2r?2 = fhinge
J2w2 + W2 X J2w2 = C2 X fhinge
If gravity is to be accounted for during the simulation, the system will auto-
matically add gravitational force terms to every objects' transiational motion
equations. Figure 4.4 shows the printed form of the equations Newton formu-
lated for the pendulum of the simulation example of Section 2.3. The terms
labeled CSjj are stabilization terms like those discussed in Section 4.5.
In contrast with equation formulation methods used by many other systems,
ours does not require special case analysis in order to deal with closed kinematic
chains. It will be seen later, however, that closed loop mechanisms sometimes lead
to constraint force values becoming indeterminate due to redundancies between
constraints. Figure 4.5 shows two four-link closed-loop mechanisms. The first,
53
Let hiz = Ci2 --H Cj1
and choose hi: and hiy to define coordinate system axes with hiz as the z axis.
position constraints:
ri + cli = r2 + C2l
hiz h2z = 0
hi: - h2y = 0
acceleration constraints:
t:1 +W1 X cii +wi x (wi xcii) = ?2 +W2 X C21 +W2 x (w2 x c2i)
--H W2 --H (wi x w2)) h2y = 0
(wi --Hw2 --H (wi x w2)) h2z = 0
constraint term variables: fhinge (3-vector),
Ty (scalars)
Figure 4.2: Pin (revolute) hinge: 1 rotational degree-of-freedom
54
Let hiz = Ci2 --H Cj1
and choose hi: and hj? to define coordinate system axes with hiz as the z axis
and let hp = (r? + c2l) --H (ri + cii).
position constraints:
acceleration constraints:
where h& =
ftp=
h?in
hi?
h2z hi: = 0
h2z			hiy = 0
hp			h1: = 0
hp			hi? = 0
(w2 --HWi --H(w2 xwi))
hp h?i: + 2hphi: + hp
(w'2 --H wi --H (w2 x wi)) hi? = 0
hi: = 0
hi: = 0
hp hiy + 2hphi? + h?p hiy = 0
(r'? +W2 X C2i --H ri --HWi xcii),
(i:2 +W2 X C2i +W2 x (W2 x c2i) --Hr4 --HWi xcii --HWi x (wi xcii)),
= Wi X hiu, and
= wi x hitt + Wi X (wi x
constraint term variables: f:, f?, T:, Ty (scalars)
Figure 4.3: Prismatic hinge: 1 transiational and 1 rotational degree of freedom
55
(?rint.object.equations `pendu1u?)
PENOULUfl
COMPONENTS OF PENflULUfl
BflLL
L3
L2
LI
L9
-t1 C + -?CFi? + M R.+. =
HCFY4 ? Ce? + J ti. w
-t1 C + HCFr? + -HCFn +
Cl! x ?CF?? + ?CFn ?< CL!
-t1 C + HCtJ! + -HcF? +
cl! x HCFY! + HCF7! x CL!
-t1 0 + HCFT2 + -HCF? +
Cli x ?CF72 + RCFri x CLI
-t1 0 + t1 R. =
J w. + W x (J W) =
HINCES OF PENDULUM
Hinge J4
Ninge J3
Hinge J2
Hinge Ji
(J W)			= 0
t1 R.. =			0
+ J ?			+ W x (J			W)			= 0
t1 ... =			0
+ J II			+ W x (J			W)			= 0
t1 ... =			0
+ J ?			+ W (J			w)			= 0
.L!			+			W.I!			x CL! + WI!			?			(WI!			?			CI! )			+
?LL			+ C?a?? ?< ?. 1111 +			(			WB?l			x			Ca??? )			x			WB?l			+			CSr? = 0
.L2			+			t?.I!			? CL! + WI!			?			(WI!			x			CI! )			+
.1! + CI!			>s W.I! +			(WI!			?< CI!			)			x WI!			+			CSY!			=			0
LI			+			?.II			x CLi + WI!			?			(WIt			?<			CII )			+
I! + CI!			x ?. 1! +			(WI!			? CI!			)			x WI!			+			CST!			=			0
+			?.Is			?< CL. + WI.			?			(WI.			x			CI. )			+
--Ha. .11 + Cti			x ?.Ii +			(WlI			? CI!			)			x WI!			+			CS;t			=			0
NIL
Dynamic Lisp Listener I
___________"? " "\			?-`-"-?-
Figure 4.4: Printed form of automatically generated equations for example of
Section 2.3
56
Figure 4.5: Two closed-loop linkages
made with four spherical joints, has determinate reaction forces. The second,
composed of pin hinges, contains a redundant constraint and so has some inde-
terminate reaction forces. This is often not a serious problem since the motion
of the object will still be uniquely determined. The issue of equation dependence
is discussed further in Section 4.6.
Nonholonomic kinematic constraints can also be handled by the system. Some
one-sided constraints (also called unilateral or inequality constraints), such as a
block sliding, are handled as special cases of holonomic constraints; since they are
piecewise holonomic, if we are able to detect the discontinuities we can model the
constraints as holonomic in the intervals between discontinuities. The more diffi-
cult cases of nonholonomic relationships, such as a disk or ball rolling on a plane,
correspond to constraints that are not expressible as equalities or inequalities be-
tween position variables. Our formulation method again handles such constraints
in a manner much like the one it uses to handle holonomic constraints. If the
constraint is given in terms of velocities, say, we differentiate once to get the ac-
57
celeration constraints. Thus, we are again able to solve for object accelerations
and can use our standard integration routines. The difficulties with this kind
of nonholonornic constraint, however, arise during the integration phase. With
holonomic constraints it will be seen that the position and velocity constraints
are used to ensure that the values produced by the integration routines satisfy the
constraints. However, much more care is needed to accurately integrate a system
containing general nonholonomic constraints. Neimark[51] presents a detailed
investigation of nonholonomic systems.
4.4 Equation storage, collection, and
modification
The basic motion equations are stored in the dynamic models of the primitives;
constraint equations are stored in the corresponding hinges. During simulation,
the integration routines need to collect the set of equations in order to calculate
accelerations and produce new positions and velocities. A distinguished set of
objects, each called a top-level object, is maintained by the system. Each top-
level object represents a maximal connected structure graph. Each primitive in
a top-level composite object is related to some other primitive in the composite
and is not related to any objects outside the top-level object. For purposes of
integration and equation solving (though not for event handling), then, the top-
level objects can each be treated individually, making for more efficient matrix
computations.
Construction of the set of equations for a top-level object is done recursively
58
on demand. A request for the equations of an object is turned into requests for
the equations of the component objects and hinges. At every level of an object's
composition hierarchy, the set of equations for that level, once derived, is stored
in the appropriate dynamic model. The system keeps track of the constraints
responsible for the various terms in the motion equations. ?Vhen, for example,
events occur that cause a hinge to be removed, the constraint equations and the
corresponding terms in the hinged objects' motion equations can be removed
without necessitating complete rederivation of the entire system of equations for
the containing top-level composite object. The equation sets of the ancestor
dynamic models need to be collected again, but the process will be inexpensive
since the other equation sets are still stored in dynamic models of interior nodes
of the composition graph.
Once the equation set has been constructed, object accelerations, hinge re-
action forces and torques, or other variables that the user has designated as
unknowns of the equations can be computed. In order set up a standard Ax = b
matrix problem, the system evaiuates the equations. Each equation is examined
in turn. Since the equation representation allows both scalar and vector equa-
tions, the system first determines how many rows of the matrix to allocate to
the given equation. Then, by examining each of the equation's terms, matrix
columns are assigned to unknown quantities and coefficients of the unknowns
computed and added into the appropriate matrix slots, while values for known
terms are computed and added into the bvector. Recall that the equations do not
contain specific numerical values, but instead refer to quantity data structures.
59
Thus, the equations do not change between discontinuities in a simulation, but
evaluation of the equations at successive time instants produces different matrix
entry values corresponding to the changes in quantity values.
Various tricks can be used to ensure that the equation evaluation proce-
dure does not spend a lot of time doing symbolic analysis and manipulation of
terms. Newton employs a simple equat?on cornpi1at?on scheme that creates corn-
piied equattons from the original equations. The terms of the compiled equations
are compiled Lisp closures that, when executed during the equation evaluation
process, efficiently perform the same matrix entry modifications that evaluation
of the noncompiled equations would. Equation compilation is generally only
worthwhile when exceptional events are rare.
Once equations have been evaluated, standard linear system solving methods
can be used to compute values for the unknowns. The values of the corresponding
quantity data structures are then set using these values.
4.5 Integration methods
In order to begin to examine the issues involved in integrating the motion equa-
tions, we first assume that all constraints are holonomic. In this case, the full set
of equations consists of second-order differential equations (the basic object mo-
tion equations) plus algebraic equations (the constraint equations). Integration
of such systems of equations is quite difficult. The difficulties and approaches
to overcoming them are discussed in numerous papers[24,31,34]. The simplest
method of integration, one we often employ, makes use only of the basic mo-
60
v a
p,v ?p
v, a			v
Figure 4.6: Simple integration method. Here, p, v, and a are the current po-
sitions, velocities, and accelerations, while p' and v' are the newly calculated
positions and velocities.
tion equations and the acceleration constraints, the second derivatives of the
constraints. This method is attractive since all we need to do is solve a lin-
ear system to find accelerations and then apply a standard ordinary differential
equation integration scheme to produce new positions and velocities. Thus, as
shown in Figure 4.6, given the object positions and velocities at a time t, the
motion equations and acceleration constraints are first used to calculate the cur-
rent object accelerations. Then, the positions and velocities are used to produce
new positions for time t + ?t, and the velocities and accelerations are used to
produce new velocities.
In Newton a number of integration methods are made available. These include
the first-order Euler method, the fourth order Runge-Kutta method and fourth
order Adams-Moulton predictor-corrector method, a hybrid method using both
the Runge-Kutta and the Adams-Moulton methods, and an adaptive-time-step
method based on the Runge-Kutta-Merson and Adams-NIoj?ton methods. The
61
hybrid and adaptive-time-step methods are used most often.
The basic fourth order Runge-Kutta formula is shown in Figure 4.7. The
fourth-order Runge-Kutta methods work well but require four evaluations of the
derivative at each step. In our situation, this entails setting up and solving the
large linear system involving the motion equations and constraints. Since the
solving of this system is often the most costly step in simulations, integration
methods that require fewer derivative evaluations make simulations mo?e effi-
cient. The Adams-Moulton method is, like the Runge-Kutta, a fourth order
method, but it requires only two derivative evaluations per step. It cannot nor-
mally be used alone, however, since it requires a history. For this reason, the
hybrid scheme uses Runge-Kutta until enough state history is saved and then
it switches to Adams-Moulton. In addition, when discontinuities in state occur
due to exceptional events, the system will need to suspend its use of Adams-
Moulton and revert to Runge-Kutta until a new history is established. The
Adams-Moulton formula is shown in Figure 4.8.
The adaptive-time-step method is similar to the hybrid method but, in ad-
dition, it mouitors the local error produced during integration steps and auto-
matically modifies the size of the simulation time step accordingly. Thus, when
conditions permit, the step size is increased to permit more efficient simulation.
On the other hand, the step size is decreased when the error bound is exceeded
(e.g. when nearing mechauism singularities).
62
= qt + --H(k1 + 2k2 + 2k3 + k4)
6
k1 = At *
k2 = At *
k3 = At *
k4 = At * q't+At?q+k3
Figure 4.7: Fourth-order Runge-Kutta integration formula.
At
qt+At = qt + --H24(55q't --H 59q'?--H?? + 37q t--H2At --H
At
qt+At = qt + (97t+At + l9q'? --H 5?'t--H?t + q't--H??t)
24
predictor:
corrector:
Figure 4.8: Fourth-order Adams-Moulton predictor-corrector integration for-
mula. The ?`t--HnAt, for i --H 1 9 3, are the saved derivatives.
63
4.5.1 Stability
The standard integration schemes outlined above are not always satisfactory.
It is well known that integration of mixed systems of differential and algebraic
equations is difficult [25]. The matrices used in solving for accelerations can have
very high condition numbers, and, in general, the equations being integrated are
stiff.
When, as described above, only the acceleration constraints are used to drive
the integration, numencal error can produce (unbounded) constraint violation
over time[73,11]. There are a number of methods of dealing with these problems,
The Gear stiff integration algorithm [24] can be used to handle some sets of
differential-algebraic equations and we are investigating its applicability to the
equation sets we encounter.
Constraint Stabilization method
The constraint stabilization method [11] is simple to implement and works well
in many situations. Use of the method only requires adding two terms to the
acceleration constraints. The method transforms an unstable constraint, ? = 0,
into a stable constraint that has the same solution, ? + kv? + k?p? = 0, where
kv and kp are user-chosen coefficients (various people [31,73] claim that values
between 1 and 50 usually work well). Intuitively, the added terms act like a
spring and damper that pull back any values that have strayed from precisely
meeting the constraints.
64
Coordinate Partitioning method
Wehage and 11aug [67] introduced another scheme, called Generalized Coordinate
Partitioning, and claim that it has a number of advantages over Gear's algorithm
and the constraint stabilization method. This method works by examining the set
of constraint equations and numerically determining a set of independent coordi-
nates. Determination of accelerations takes place as usual but integration to get
new positions and velocities is done differently. When using the simpler integra-
tion methods, all of the current accelerations are used to produce new velocities,
and the current velocities are used to produce new positions. As indicated, if
the original constraint equations are ignored, this method does not ensure that
the constraints will remain satisfied. The coordinate partitioning method uses
only the accelerations, v?, of the chosen independent position variables, v, to get
new velocities, v. Similarly, the current velocities of the independent variables
are integrated to produce new positions. The positions, u, and velocities, u, of
the dependent variables can then be determined from v, 1), and the constraints
by using iterative procedures such as Newton's method. Thus, the coordinate
partitioning method can ensure that the constraints remain satisfied to within a
user-specified tolerance.
There are some complications with the method. The choice of independent
vanable set affects the accuracy and convergence properties of the method. The
set must be monitored during simulation to determine when a new partitioning
of coordinates is required. This kind of task can be handled by Newton's event
handling mechanism.
65
11aug et al. [34] have also developed a hybrid method combining constraint
stabilization and coordinate partitioning, since both have advantages in particu-
lar situations.
4.6 Dependence
As mentioned in Section 4.3, the use of our method of formulation will sometimes
produce a dependent set of object motion equations and acceleration constraints.
A number of formulas have been developed in the mechanics domain to aid
in the determination of the number of degrees of freedom in a mechanism. They
generally are grouped into a class known as Greubler formulas[57,59]. For ex-
ample, for 3-space mechanisms containing ouly one-degree-of-freedom joints, the
formula for the total number of degrees of freedom in the mechanism is 6n--H5h--H6,
where n is the number of rigid bodies in the object and h is the number of hinges.3
?Vhile such formulas are useful in many cases, they must be used with care since
there are also many exceptional cases for which the formulas do not give the
correct result. Consider again the kinematic loop of Figure 4.9, this time with
the bottom object fixed in space. The mechanism is a standard four-bar linkage.
The Greubler formula predicts 6 * 4 --H 5 * 4 --H 6 = --H2 degrees of freedom for the
object. Clearly, however, the mechanism has one degree of freedom. The dis-
crepancy is due to the fact that the twenty constraint equations contributed by
the four pin hinges are dependent. To see this, imagine the loop being broken at
hinge 4. The other three pin hinges are enough to restrict the attachment points
3This assumes that exactly one of the n bodies is fixed in space.
66
fixed
Figure 4.9: Four-bar closed loop linkage with pin hinges and bottom object fixed.
of hinge 4 to the proper plane and hinge 4's axis to the proper orientation. This
makes three of the five pin hinge constraints redundant (1 for the restriction to
plane, 2 for the axis orientation). Thus, we get the right number of degrees of
freedom with 6 * 4 --H 3 * 5 --H 1 * 2 --H 6 = 1.
Depending on the experimenter's focus, this type of dependency might or
might not be problematic. When only the motion of the mechanism is of interest,
the redundancy is not a problem. The accelerations in many cases are still
ui:iiquely determined; they are determined when there is no friction, but not
always in situations involving friction. However, redundancies cause trouble when
hinge reaction forces are of interest, since while certain of the constraint forces
may still be uniquely determined, the forces corresponding to the redundant
equations are indeterminate. Thus, in the example, the values of constraint forces
are individually indeterminate. The total value of the these forces is known, of
course, but it can be spread among the various hinges in an infinite number of
ways.
67
Th			____
-1
Figure 4.10: Four-bar linkage in bifurcation position and two possible later con-
figurations.
Further dependencies can arise in other situations. For instance, when the
four-bar linkage (counting the ground as the fourth bar) is in the position shown
at the top of Figure 4.10, only sixteen of the twenty acceleration equations derived
from the hinge constraints remain independent, rather than the usual seventeen.
Again, this doesn't pose any particular problem when the dynamics of the prob-
lem is being calculated and the individual constraint forces are not needed. If a
kinematic analysis is to be done however, this configuration represents a singu-
larity where the mechanism has, instantaneously, two possible motions (toward
the bottom two configurations shown in Figure 4.10).
Finally, there are times when singular configurations and redundancy, can
lead to trouble if the wrong integration scheme is chosen. Examine the situation
68
Figure 4.11: A critical form
in Figure 4.11, with the leftmost and rightmost dark objects being fixed in space,
and all hinges being spherical hinges. The acceleration equations for the hinge
constraints are redundant but solvable. Thus, given certain forces on the middle
two rods, the motion equations can be solved to produce nouzero accelerations for
the rods. Since the configuration should obviously be rigid, something is wrong.
If an integration routine is chosen that uses only the acceleration constraints and
makes no use of the position or velocity constraints, the computed acceleration
values will be used to produce new nonzero velocities and positions for the rods,
contradicting the necessary rigidity. (Usually, the linear system constructed from
the dynamics equations will then be inconsistently singular on the next attempted
simulation step.) The coordinate partitioning method, by making use of the
position and velocity constraints, can generally do the right thing, since it uses
those constraints to determine an independent set of values to be integrated.
Since the size of this set corresponds to a numerically computed value for the
system degrees-of-freedom, in this case 0, the method will be able to "detect"
and ensure the rigidity of the mechanism. This example is an instance of a kind
of singularity called a cmticai form. Some discussion of the characterization of
69
critical forms and other singular and dependent configurations can be found in
Paul[57].
4.7 Other formulation schemes
Throughout this thesis, discussion of the formulation of object motion equa-
tions has concentrated on the use of the Newton-Euler formulation method. In
addition, the general method described in this chapter makes no use of special
object structure to produce computationally efficient formulations. In particular,
many other systems take advantage of the fact that many mechanisms contain
no closed kinematic loops. For such mechanisms, other methods, like the recur-
sive Newton-Euler formulation, can be used and result in smaller equation sets
(having, for instance, a number of equations equal to the number of degrees of
freedom in the mechanism). Featherstone examines a number of these methods
in [23]
Our method was chosen because of its simplicity and generality and because
it supports the incremental modifications to the equation set required when con-
strained relationships between objects change. Hybrid methods could be used
that combine, say, a method for tree-structured objects and a "loop-closure" tech-
nique. ?Vith such a method, certain hinges in a closed loop are "cut" until the
object becomes tree-structured. The motion equations are then formulated for
the new tree-structured object Ultimately the loop-closure technique accounts
for the cut hinges by modifying the equations produced from the tree-structure
analysis in a manner similar to Newton's basic formulation technique. In essence,
70
Newton's formulation method is one large loop-closure scheme. All the hinges
are cut initially, then accounted for by the loop-closure scheme.
Along with various Newton-Euler-based methods of motion equation formu-
lation, Lagrangian methods of equation formulation must be considered. There
are a number of advantages and also some disadvantages to the use of Lagrangian
methods of motion equation formulation. One advantage of the Lagrangian
method is that is a scalar rather than vectorial method, making the mathematics
somewhat cleaner and simpler. On the other hand, if, as is often the case, it is
inconvenient to choose an independent set of coordinates (corresponding to the
number of degrees of freedom in a mechanism), or if hinge reaction forces are of
interest to the experimenter, then the standard Lagrangian formulation method
must be modified to include constraint equations and Lagrange multiplier terms.
This Lagrange multiplier form of the Lagrangian formulation method is similar to
our formulation --H the terms we add to object motion equations due to hinge con-
straints correspond to the Lagrange multiplier terms of the Lagrangian method
and are derived in an analogous manner. Most of these equation formulation
issues are beyond the scope of this thesis, but are being considered in further
work on architectures for simulation systems.
Chapter 5
Newton's Analysis Module:
Handling of Exceptional Events
General purpose simulation systems need the ability to respond to a variety of
exceptional events. The first section of this chapter presents a uniform mecha-
nism for coordinating the handling of such events. Subsequent sections describe
Newton's basic collision detection and impact resolution schemes, its handling of
continuous contacts between objects, and its handling of frictional effects.
5.1 Handling of exceptional events
During the course of interesting simulations, a variety of events can occur that re-
quire special processing. Collisions need to be detected and resolved, constraints
modeling contacts between objects need to be added and deleted, friction forces
need to be monitored to determine when an object changes from sticking to
71
72
sliding on a surface, and so on. In general, exceptional events can cause dis-
continuities in object velocities or accelerations and necessitate corresponding
modification of the internal models of object behavior. It is crucial to the suc-
cess of general-purpose simulation systems that they be able to deal with such
events.
Previous simulators have done little to provide an event handling mecha-
nism. Once the set of equations used to model object behavior is formulated,
whether automatically or not, it is generally fixed. Furthermore, these equations
are assumed to represent a continuous function of the object state variables. To
some extent, handling of impact and control can be incorporated into a scheme
requiring continuity. Some systems deal with impact by including "force-field"
constraints in the equation set. These constraints, while continuous, have neg-
ligible effect when the objects are not in contact but exert large repulsive force
on the objects when they are very near to contact. Control problems are dealt
with in a similar manner by some researchers. Equations modeling aspects of a
control strategy that should be active only during a particular time interval of
the simulation are rewritten to be continuous but to have negligible effect outside
the interval of interest. Work on this kind of scheme for the handling of events,
where discontinuities are approximated by continuous functions, is presented in
[20,53]. Explicit detection and resolution of specific kinds of discontinuities is
examined in [66,35,23].
A more natural approach to dealing with such events involves the inclusion
of a general purpose event handling mechanism. The event handler serves as an
73
interface between the top-level loop of the simulation, responsible for equation
formulation and integration, and a set of subsystems responsible for detecting
and then handling exceptional events.
Newton has a general purpose event handler that is currently responsible for
coordinating:
0 detection and resolution of impacts,
o+ analysis of continuous contacts between objects and corresponding main-
tenance of tempora? hinges, special kinds of hinges that model one sided
constraints between objects in contact, and
o+ detection of control events (object state dependent conditions usually sig-
nifying state changes in user-supplied control programs) and corresponding
modification of the control equation set.
Since it deals with discontinuous changes in system behavior, the event handler
is also responsible for things such as restarting parts of the integration subsys-
tem. For instance, one of the normally used integration schemes starts by using
a Runge-Kutta method for the first few simulation steps. ?Vhen enough previ-
ous states have been saved, the system switches to using an adaptive-timestep
Adams-Moulton predictor-corrector method. However, impacts, formation or
breakage of temporary hinges, control-program directed modification of the set
of motion equations, or other events, generally cause discontinuities in object
velocities and accelerations, invalidating further use of the saved states. Thus,
when the event handler detects such an event, it makes sure that the integration
74
module will restart itself the next time that it is invoked, probably by reverting
to use of the Runge-Kutta method.
5.1.1 Details of the event handling mechanism
Figure 5.1 shows an outline of the event handling interface. The specific event
handling modules are called in order to determine what events, if any, have
occurred. Descriptions of these events, indicating the particular subsystem re-
sponsible for the event, the type of event, and whether or not it is necessary to
compute a closer approximation of the time of the event occurrence, are collected
by the event handler. We currently distinguish between two types of events, those
that cause velocity discontinuities and those that cause acceleration discontinu-
ities. Events, such as impacts, that cause velocity discontinuities are normally
handled by deriving a set of impulse-momentum equations from the dynamic
equations of motion and solving them to ensure momentum conservation while
causing instantaneous velocity changes in the object states. Acceleration discon-
tinuities generally involve only removal or addition of acceleration constraints on
the objects. These events are handled more simply, mainly requiring reinitializa-
tion of the integration module.
In many cases, the time at which the occurrence of an event is detected is not
a close-enough approximation of the actual time of its occurrence. For instance,
the integration step might cause two previously noncontacting objects to over-
lap. The event detection routines can specify tolerances for preciseness of event
occurrence. Given the complete set of events collected from the detection rou-
75
;;; Assume consistent previous state exists, say for time t,
;;; and that a new state has just been proposed for the current time, t + At.
;;; Upon leaving the function, a consistent new state will have been
;;; established for some time in the interval (t,t + At].
(defun handle-exceptional-events ()
;; Collect the various types of events.
(collect-collision-events)
(collect-contact-events)
(collect-control-events)
;; Examine the events.
;; Search interval between previous state and current
proposed state to determine time of earliest event.
;; Use interpolation scheme to compute state for that time.
(??n?-time-o?rst-event?
;; If there are velocity discontinuity events, process them
;; first. Events will need to be recomputed after velocity
;; modification.
(handle-velocity-discontinuities)
;; Given the new velocities, calculate corresponding
;; new accelerations. Then determine and handle any
;; events causing acceleration discontinuities.
(calculate-new-accelerations)
(recollect accderation-discontinm.ty-events)
(handle-acceleration-discontinuities))
Figure 5.1: Newton's event handler interface
76
tines, if any event tolerances are not met, the event handler will use information
supplied in the event descriptions to search the interval between the previous
time instant and the current one for a time satisfying the tolerances.
Typically, a binary search of the time interval is done. The search ends at the
earliest time where at least one of the detected events does occur and satisfies its
tolerance. (Clearly, such a time exists; if we go back far enough in the interval,
none of the events occur since they didn't occur at the previous time instant.)
The system state must be computed at each time instant checked by the search.
Instead of using the normal integration routines, some efficiency is gained by a
method that interpolates between the previous state and the state at which the
events were originally detected. In addition, during the search, the set of events
that need to be checked grows smaller since some of the originally detected events
drop out as time is backed up. This simple method of determining the time of
event occurrence has worked sufficiently for our purposes. Possible improvements
are discussed in the sections on impact and contact and in Chapter 7.
Once the event occurrence time is isolated, the appropriate action can be
taken. If both velocity discontinuity and acceleration discontinuity events oc-
cur, the velocity events are handled first. Acceleration events generally need to
be recomputed following the state modification caused by the velocity events.
Again, handling of the velocity events usually involves formulation and solution
of impulse-momentum equations and will be discussed in the section 6n impact
and in the section on control. After resolving the velocity events, a complete,
though not yet necessarily event-free, state is again established by computing new
77
accelerations in accordance with the positions and new velocities. If this state
yields events causing acceleration discontinuities, these are handled to produce
a final consistent state.
A number of issues complicate this scheme for handling events. Depending on
how the particular methods of handling velocity or acceleration discontinuities
work, it is possible for new velocity or acceleration events to occur at the current
time instant due to the instantaneous velocity increments caused by handling the
velocity events. Thus, there is the temptation to iterate the handling processes
at an instant. That is, the event handler could iteratively detect velocity events
and then resolve them, terminating ouly when no further velocity events occur.
Following resolution of the velocity events, a similar strategy could be followed
for acceleration events; events would be detected and then resolved until no more
events were found. There are a number of disadvantages to using such a strategy.
It is counterintuitive, allowing a large set of states to exist and then be modified
in what is supposed to be a single time instant. Contrived control strategies
could benefit from such a scheme, but control should probably not be so closely
linked with the implementation of the event handler. Overall, such a scheme
would provide extreme generality with few clear benefits.
This is not to say that iteration schemes cannot be a part of exceptional event
handling, only that our top-level event handling mechanism does not iterate at
a time instant. The actual event handling functions (handle-velocity-events
and handle-acceleration-events) can handle the events in any manner, in-
cluding use of iterative schemes. In fact, Newton's standard impact resolution
78
method, detailed in the next section, is iterative. Once the impact resolution
procedure has returned control to the top-level event handler interface, however,
no further resolution of impacts can be done for that particular time instant.
There are still a number of difficulties with our event handling scheme. The
noniterative scheme can leave events unresolved. If such events are discovered
during the next simulation step, the search procedure for isolating the event
occurrence time can lead back to the current time. This presents a problem if the
event occurrence tolerances are not then met; since the current time represents
a discontinuity in the simulation, the search can go back no further. Of course,
if the tolerances are met, the events can be handled. In such cases, however,
it might be more appropriate to signal an error or to ignore the event since the
situation is likely to have ansen due to an error on the part of a control algorithm
programmer.
Another difficulty with the event handling mechanism lies in the potential for
missing events. If particular events are to be detected by, for instance, changes
in the sign of certain functions, then the assumption must be made that small-
enough time steps will be taken that no sign changes will be missed. In order
to deal with this problem, the individual event detection modules should be able
to generate events that, when handled, influence the choice of integration time
step size. Thus, as collision detection routines notice that collisions might occur
in the near future they could generate events to decrease the time step size.
79
5.2 Impact
For the purposes of this section and that on contact, we restrict the following dis-
cussions to impacts and contacts between polyhedral objects, though the Newton
system is not restricted in this way. Then, the various contacts may be described
as surface-surface, edge-surface, point-surface, edge-edge, and so on. Impacts are
distinguished from other contacts as those contacts where the velocity of a con-
tact point on one object relative to the corresponding point on the second object
is directed into the second object's interior. For a contact between a point, Pi,
of an object Ol and a corresponding point, P2, on the surface of another object
02, the condition is stated more precisely as
(Pi --H P2) ?contaa < 0,
where flcontact is the normal to the contact surface (directed toward 02's exterior).
When the normal-direction relative velocity is greater than zero, the contact is
in the process of breaking. When the velocity is zero, the contact will remain
and may result in the creation of a temporary hinge.
When the event handler begins its impact analysis, the integration module
has just proposed a set of positions and velocities for time, t + ?t. The sys-
tem then uses its geometric modeling subsystem to determine whether any im-
pacts occurred in the time interval. While there are many difficulties in properly
computing the intersection between two geometric representations [44,401, the
problem of determining the precise time of any impacts makes matters still more
complex. To do things correctly the four-dimensional space-time swept volumes
80
of two objects must be intersected. In the current implementation, however, we
count on time steps being sufficiently small that we don't miss collisions between
steps and, when it is determined that an impact does occur between times t and
t + ?t, the moment of impact is found by binary search of the time interval.
We repeatedly halve the time increment, reintegrate, and analyze contacts for
the new time, until any object interpenetrations are within a user-controllable
tolerance.
After determining the time of any impact occurrences, the geometry system
is used to analyze the nature of all inter-object contacts for that time. For the
moment, assume that all such contacts are indeed impacts. To resolve impacts,
Newton formulates impulse-momentum equations for each object, and contact-
point velocity equations for each impact (using coefficients of restitution based on
object properties), and then solves this equation system to compute appropriate
instantaneous changes for the object velocities.
The equations are automatically derived in a fashion analogous to the formu-
lation of motion equations described earlier. For point-on-surface impacts, the
process involves formulating equations of the following form:
m1Arj			=			fimpactncontact
=			ci x (fimpactncontact)
m2Ar2			=			fimpactncontact
J2Aw2			=			C2 X (f?mpactncontac?)
A			A
Pi --HP2 =?e(pB1?p2B),
81
where Cj is the vector from object i's center of mass to the location of the hinge,
fhinge is the (scalar) impact impulse, flcontact is the surface normal at the point
of impact, Ar? is the difference between object i's center of mass velocity before
and after impact, p1B. and pAj are the velocities of the impact point on object i
before and after impact, respectively, and e is a coefficient of restitution that
depends on the material properties of the colliding objects.
When composite objects are involved, impulses due to impact are transmitted
through hinges by formulating impact constraint equations for the hinges and
adding appropriate impulse terms to equations for the hinged objects.
Thus, if object 2, from above, and a third object are related by a spherical
hinge, the hinge equation is
and the object equations are
A + Ax C2 = A + Ax C3
r2			w2			r3			w3
m2A?2			=			fimpactncontact + fhing?
J2?w2			=			C2 X (fimpactncontact) + c2 x fhinge
?3Ar3			=			?fhinye
J3?w3			=			C3 X --Hfhin9e,
where fhinye is a global coordinate system vector representing the impulse trans-
mitted through the hinge.
Our current model of impact is extremely simple and clearly not satisfactory
in some situations. When all of the contacts are impacts, use of this collision
82
,,, The state has been saved for time t.
New positions and velocities have been proposed for time t + ?t
1. Compute the earliest impact time ti ill [t, t + ?t].
2. Integrate from time t to produce state for time ti.
3. Determine and analyze all contacts.
4.
WHILE there are still impacts in the contact set DO
4.1 Formulate equations for resolving the impacts.
4.2 Solve the equations to obtain new velocities.
Figure 5.2: Impact resolution scheme
resolution scheme yields instantaneous velocity changes that do not imply any
further impacts. That is, after the impact resolution procedure there are no
longer any contacts that meet the conditions for being impacts. In many situ-
ations, however, only a subset of the set of contacts represent true impacts. In
dealing with such situations, we do not produce impulses for the non-impact con-
tacts. Thus, the impact resolution scheme treats these contacts as if there were
an infinitesimal separation between the objects. Resolution of a set of impacts
then often produces post-impact velocities for which some new subset of the con-
tacts represent impacts. To handle this, we currently iterate the procedure until
the set of contacts contains no impacts. Newton treats this entire process of solv-
ing a sequence of instantaneous impacts problems as occurring instantaneously.
The impact handling scheme is summarized in Figure 5.2.
We use this impact model at present because it was relatively simple to imple-
ment and produces expected behavior in many cases. For instance, in Figure 5.3,
using our model to resolve the collision between spheres 2 and 3 yields the nor-
83
A			?
before impact
A
after impact: case 1
after impact: case 2
Figure 5.3: Different collision interpretations
mally expected behavior in which, after impact, spheres 2 and 3 are at rest and
sphere 1 moves off to the left. Other simple impact models can be used in Newton.
Featherstone [23], for example, details a different scheme for resolving impacts
in the presence of contacts. In Featherstone's scheme, impulses are transmitted
through the non-impact contacts, soit would appear to be a more realistic model.
However, under this model, the spheres of the example behave in a less expected
manner. After impact, sphere 3 moves back to the right and spheres 1 and 2
behave as if they were connected by a true hinge, moving off to the left. Still,
the model does produce better results than ours in other cases.
One case that our model does not handle satisfactorily involves a large mass
block falling onto a smaller mass block that is resting on a table top. See Fig-
ure 5.4. Though modeled as an instantaneous event, a large sequence of impact
equation solutions will be executed before the large block rebounds upward.
Even for the least troublesome case, in which the collisions are considered to be
perfectly elastic, the size of the sequence is proportional to the square root of
84
v
Figure 5.4: Collision requiring large sequence of impact equation soii?i ????is
the ratio of the two masses. Thus, for these situations, our model is both time
consuming and counterintuitive.
In addition, there are situations involving multiple concurrent impacts for
which the current implementation of our model is clearly unsatisfactory. ?Ve
model surface-surface impact as a set of point-surface impacts --H one set of
point-surface impact equations is set up for each corner of the convex hull of
the contact region (shown as a dotted rectangular region on the fixed slab in
Figure 5.5). Thus, impulses are, except in cases of redundancy, generated on
every one of the modeled impact points. In the situation shown in Figure 5.5,
where a falling slab is striking another, fixed, slab, this formulation can result
in an inconsistent set ofequations. The solution to this is the use ofthe kind of
technique employed by Featherstone's impact resolution method to determine a
subset ofthe impact points on which to generate impulses.' This modified ver-
sion ofour implemented impact resolution technique, then, produces consistent
1The techniquefor choosing the proper subset iscompletely analogous to those described in
the next section forchoosingcontact support sets.
85
K
k			fixed
fixed
Figure 5.5: C()ii1?1ull for which tliesimplest model formulat??s i?ic??nsistent impact
equations.
equation sets while retaining the desirable behavior of the current model. It is
still iterative and does not generate impulses on the non-impact contacts. So,
for example, the behavior for situations like that of the colliding balls discussed
above remains the same.
It is clear that neither our model nor Featherstone's model is sophisticated
enough to do realistic impact modeling. ?Ve are currently investigating more
complex models that can better account for the elasticity properties of objects.
5.3 Contact
Newton was designed to be able to handle continuous contacts between objects.
By continuous contacts, we mean contacts in which two objectsremain in contact
for a finite (rather than infinitesimal, as for impacts) amount of time. Such
86
contact relationships as in a block sliding or a ball rolling on a table top are
modeled in our system by extensions to hinges called temporary hinges.
When a continuous contact is made between two previously noncontacting
objects, a new kinematic relationship is formed. Constraints corresponding to
contact relationships are part of a subclass of nonholonomic constraints called
unilateral or one-sided constraints[57]. These constraints are modeled by inequal-
ities rather than equations. Recall that the most easily handled constraints are
the holonomic constraints, represented by equations relating position variables.
If the variables related by an inequality are position variables, then the inequality
falls into a special class of nonholonomic constraints called piecewise holonomic
constraints. For piecewise nonholonomic constraints there are simply different
holonomic constraints active during different time intervals.
Newton's temporary hinge data structure is a subclass of the hinge data struc-
ture used to model contacts and other relationships that can be represented by
inequality constraints. Temporary hinges are not restricted to the class of piece-
wise holonomic constraints, but that is the class that will be focused upon here.
For most simulations we want to enforce noninterpenetration constraints be-
tween objects. Rather than explicitly creating the large (approximately 0(n2))
set of inequalities modeling these constraints, Newton uses its event handler and
geometric modeling system to monitor and manage contact relationships between
objects. During the geometric analysis of contacts, normal-direction velocities
of contact points are monitored. As covered in the previous section, negative
contact velocities correspond to collisions and are handled accordingly. On the
87
other hand, when contact velocities are zero2 there is continuous contact and the
system creates a temporary hinge to model this relationship. During the course
of simulation, the system's event handling module continually monitors contact
velocities, removing temporary hinges when the contact constraints are no longer
met.
For instance, when a vertex of an object, 01, achieves continuous contact
with a plane of another object, 02, the kinematic relationship that should then
be maintained is
(pi --H pcontact) ?contact > 0
where pi is the position of the contact point of 01, Pcontad is a point on the
contact plane of 02, and fleontact is the normal to the contact surface (directed
toward 02's exterior). Thus, the contact point must remain either on or above
the contact plane on 02.
Because of the intended inequality constraints, determination of object accel-
erations is made complicated when temporary hinges exist. Constraint equations
for temporary hinges are formulated in the same manner as for other hinges, and
constraint force terms are again added to the motion equations of the hinged ob-
jects. For instance, for point-on-plane contact without friction, the instantaneous
acceleration constraint is
(Pi --H P2) flcontact + 2(?i --H p2) (w2 x flcontact) = 0.
where pi and P?i are the velocity and acceleration of the contact point on 01,
2Within some epsilon, of course. For purposes of this section we avoid the crucial issue of
numerical difficulties.
88
and P2 and P?2 are the velocity and acceleration of the contact point on the
contact plane of 02. However, using equality constraints instead of inequalities
to solve for object accelerations necessitates checking the results for consistency.
Since the equation solving procedure calculates values for hinge reaction forces
in addition to calculating object accelerations, the system is able to check that
the values of the reaction forces for any temporary hinges are consistent with
that hinge's intended inequality constraint. For the point-plane contact case the
system needs to check that the normal-direction component of the reaction force
is not tensile, since a contact hinge should only sustain compression.
For two convex polyhedral objects in contact, the region of contact will be
either a polygon, a line segment, or a point. For the case of a polygonal region,
it is sufficient to use ouly vertices of the polygon's convex hull in formulating the
temporary hinges constraints. The system first assumes that polygonal support
will be maintained and searches for a "support triangle" among the convex hull
vertices. If no support triangle produces accelerations and reaction forces consis-
tent with the contact conditions, the system successively searches for supporting
contacts having more degrees of freedom, i.e. it attempts to find a supporting
line segment and, failing that, a single support point. If no set of contact points
yields consistent solutions, the system will remove the hinge. In the simula-
tion of Figure 5.6, the kinematic relationship between the small and large blocks
changes twice. After sliding across the top of the large block and maintaining
plane-plane contact, the kinematic relationship changes to plane-edge contact
as the small block tips over the end of the larger block. Ultimately, contact is
89
In -ff1			___
broken altogether.
Figitre 5.6: Changing kinematic relationships
Practical contact support set determination
For a single pair of objects in contact, determination of a consistent set of support
points is simple, taking time (at worst) linear in the number of contact vertices.
However, for the case of multiple objects and contacts, a naive algorithm postu-
lates support sets for a contact between a pair of objects independently of sets
proposed for other contacts. The result is an exponential search of the space of
possible contact support sets. Figures 5.7 and 5.8 show a simulation of three
bricks stacked on a table in an initially stationary but unstable position. The
lower picture of Figure 5.7 indicates the contact region boundaries with dot-
ted lines. Since there are three rectangular regions, the naive contact support
90
set search algorithm might check a few thousand possible support sets.3 Since
simply checking the consistency of any particular candidate support set involves
solving the dynamics equations to determine accelerations and reaction forces,
this process can be prohibitively expensive even for small examples.
The exponential complexity exhibited by the naive contact set search algo-
rithm can be avoided in a number of ways. Heuristics and continuity arguments
can be used to guide the search, or a quadratic programming or linear comple-
mentarity problem can be solved to determine the appropriate support set.
Quadratic programming and linear complementarity techniques
As detailed by L5tstedt[46,47] and Featherstone[22,23], the problem of determin-
ing contact support sets can be formulated as a quadratic programming or linear
complementarity problem. The development proceeds by first observing that a
(nonzero) reaction force occurs at a contact point if and only if the contact is not
breaking, and then combining this with the standard dynamics equations. The
resulting problem is of the form
minimize xTAx + xTbI
subject to x > 0
Ax + b > 0.
3More precisely, the simplest algorithm searches through
((43)+(42)+(41)+(4o)f=2744
possible support sets. A slightly smarter, but still exponential algorithm would examine up to
(4 + 4 + 4 +			= 2197 support sets.
91
L
F
:?`ire 5.7: Stack of bricks. Dotted lines are contact re?i()n b??r?l??i
mm			yMmTh
Mm			mm			mm			?
?Th			?
Figure 5.8: Stack of bricks simulation
92
Here, x is a size n vector representing the magnitude of the reaction force at each
contact point, while A is a positive definite or positive semi-definite matrix de-
rived from the object masses and inertias, the noncontact kinematic relationships
between the objects, and the contact normals. Though the general quadratic pro-
gramming problem is NP-hard, the definiteness properties of A make the above
particular form of the problem solvable in polynomial time; practical algorithms
that are guaranteed to run in polynomial time, however, do not yet exist. The
problem can be solved efficiently in practice by a variety of algorithms that are
exponential in the worst case.
Heuristics
In many situations, heuristics can be used to guide the search for contact support
sets. Featherstone[22,23] presents a heuristic search algorithm that is efficient
in practice on a subset of the general problem. The subset consists of those
situations for which there exist unique consistent distributions of reaction forces
among the set of support point candidates (i.e. the vertices of the convex hull
of the contact region). In Figure 5.9, when the leftmost object is resting on
the table, there is a unique distribution of supporting reaction forces among the
three vertices of the contact region. For the cube on the right, however, the
reaction forces can be distributed in infinitely many ways among the corners
of the contact region.4 Featherstone's heuristic search algorithm can only be
used when the forces are uniquely determined, though he has developed some
4The only restrictions are that diagonally opposite corners must be given the same force, all
of the forces must be compressive or zero, and their sum must total mg.
93
Figure 5.9: Contact regions. For the triangular region, a unique supporting
distribution of reaction forces exists, while for the rectangular region infinitely
many distributions exist.
unpublished ideas [21] for the nonunique case.
Baraff[7] has recently suggested a technique that approximates the true prob-
lem by a linear rather than quadratic programming problem. While it is not intu-
itively clear that the chosen objective function closely approximates the intended
quadratic objective, the method appears to work well in practice. Construction
of situations for which the heuristic fails is not trivial and the resulting con-
trived configurations will probably be encountered only rarely. The method has
not been compared experimentally with the quadratic programming techniques,
which, when set up appropriately, may work just as well in practice.
Continuity arguments
At certain points in a simulation, such as at the very beginning and after velocity
discontinuities caused by impacts, no useful past history is available. For these
situations, determination of contact support sets relies on the methods outlined
above. However, throughout much of a simulation, when motion is continuous,
the system's ability to monitor contact positions, velocities, accelerations, and
94
reaction forces, can make the complex search techniques unnecessary.
First, since kinematic relationships do not usually change very often, an al-
gorithm will normally attempt to use, for each contact, the same support points
as for the previous time instant. Furthermore, by continually monitoring con-
tact reaction forces and their derivatives, information about which contacts will
break and when they will break can be inferred. For instance, if at time t + ?t,
while checking the consistency of the support set used at time t, it is noticed
that exactly one reaction force is negative (i.e. tensile, and so, illegal), then it
is clear that that reaction force is zero at some moment in the interval. An ex-
ceptional event can be generated in order to isolate the zero point, after which
the appropriate constraint changes, such a removal of a temporary hinge, can be
made. Gilmore examines similar continuity arguments in his work, though he
concentrates on two-dimensional problems[27,28,26].
5.4 Friction
Contact problems involving friction can be handled by methods similar to those
described in the previous section for dealing with frictionless contacts. L5tstedt
[46,47] includes consideration of friction in his formulation of quadratic program-
ming and linear complementarity problems for modeling contact, and discusses
some of the particular difficulties of handling friction.
In Newton, only a simple Coulomb friction model has been developed. Tem-
porary hinges modeling point-on-plane contact with friction have two associated
sets of constraints; one set is used when sticking occurs, the other during sliding.
95
Based on the material properties of the objects in contact, there is a static coef-
ficient of restitution, ?static, that limits frictional forces during sticking, and a ki-
netic coefficient of friction, (`ikinetic), that determines frictional forces during slid-
ing. The sticking constraint is the same as the spherical joint constraint described
in Chapter 4. The sliding constraint combines the frictionless contact constraint
of the previous section with a friction constraint, ffriction = ?kineticfnormaI, that
determines the magnitude of the friction force at the contact.5
The event handler can monitor conditions on the constraints to detect when
transitions from one constraint set to the other occur, and can then alter the
active constraint set accordingly. The transition from sticking to sliding occurs
when the frictional component of the contact reaction force becomes greater than
the product of the static friction coefficient and the component of the reaction
force along the contact surface normal, i.e. when ffriction > ?staticfnormai The
transition from sliding to sticking occurs when the velocity of the contact along
the plane returns to zero.
Figures 5.10 and 5.11 show a simple simulation in which a block slides on a
table. The mass of the block is 48 kilograms, and the static and kinetic coefficients
of friction are .9 and .6, respectively. Initially the block is fixed and gravity is the
only external force affecting it. After one second of simulation, an applied force
whose magnitude increases exponentially with time begins to push on the block in
the x direction. Throughout the simulation the normal direction reaction force is
mg = 48*9.81, or 470.88 newtons. ?Vhen the static friction force limit is reached,
5The direction of the force is the opposite to the direction of the velocity of the contact along
the plane.
96
the block begins to slide. This is seen at approximately the 2.7 second point
in the CTotal frictional force'' graph, where the friction force changes
discoutinuously from a value of .9* rng to a value of .6 * rng. After 3.25 seconds,
the applied force becomes constant and at time 3.75 it is removed. Following
removal of the force, friction causes the block's velocity to decrease linearly. The
block comes to a stop and reverts to the sticking state after about 5.6 seconds.
97
Siflu1ation
z
Figure 5.10: Simple simulation involving friction.
0.0 1??63 3.121 4.630 6.21
1.163			3.121
Ii			ction fo
00k - 0 pos?t?on			oto norna?--H
Block--HXvelocity			Totalfrictionalforce
11.0 --H			900.0
7.1 --H			???o?			?!!
0.0 .0 I
0.0			i?163			3?121			4.600			6.21			0.0			1.163			3?12
Olock - 0 ?cceler4tion			o+?pplieda t?rce on bl?cK
Figure 5.11: Graphs for the simple friction simulation
Chapter 6
Motion Control
Motion control is an important aspect of the design and analysis of complex
mechanisms. This chapter begins with a development of interface mechanisms
that allow motion control to be integrated with dynamic analysis. Approaches
to high-level control of complex objects are then examined in the context of some
example simnlations.
6.1 Introduction
One of the primary goals of the development of Newton was to build a system
that provided support for experimentation with design. The ease with which
users can modify the geometric specifications of simulated objects and the kine-
matic relationships between them is one aspect of such support. For observing
the behavior of simple systems, such as pendulums or bouncing balls, this might
be enough. However, complex objects, such as humans or robots do not just
98
99
react passively to forces in the surrounding environment. Their complex, inter-
esting motions are a result of an internal control mechanism that supplies control
forces and torques in order to achieve particular behavior. Standard robot arms,
for example, have controllers that compute joint torques that are then applied
by motors in the joints to achieve a specified end-effector trajectory or other
commanded motion.
Newton's control subsystem supports a variety of methods which experi-
menters can use to influence or direct the motion of simulated mechanisms.
6.2 Control quantities in Newton
Chapter 4 developed the formulation of the set of equations describing the mo-
tion of simulated objects. For a mechanism composed of n primitive objects and
m hinges each constraining hi degrees of freedom, this set consists of 6n + ?2m=1 hi
equations. By choosing as unknowns the 6n translational and angular acceler-
ations and the zm?=1 hi hinge reaction forces and torques, a determinate system
of simultaneous linear equations can be formed and solved. The computed ac-
celeration, force, and torque values are then used by the integration module to
produce object states for a new time instant.
In order to control the motion of objects, degrees of freedom must be reintro-
duced into the formulated set of equations. This is achieved in Newton through
the definition of control quantities. Control quantities are simply standard force
or torque quantities (see Chapter 3) whose use is specifically related to motion
control. ?Vhen a user defines a control quantity and associates it with a par-
loo
ticular object or hinge, terms for the quantity are added into the appropriate
dynamics equations.
As a simple example, consider an unconstrained block in space. The two
vectorial equations
mr = mg
Jw + ?`? x Jw =0
represent 6 scalar equations in 6 scalar acceleration unknowns. To influence
the motion, a user could define a vectorial control force quantity, j?on?roi, and
associate it with the block. The equations then become
mr = mg + fcontroi
Jw +w x Jw = 0.
By viewing the control force as an additional unknown, three degrees of freedom
are added back into the system of equations.
Similarly, experimenters often wish to model the inclusion of actuators in the
joints of complex objects. For two objects constrained by a pin hinge, a control
torque quantity can be created and associated with the hinge. The motion equa-
tions of the two hinged objects are then modified to include terms representing
the actuator control. Let Zhinge be the hinge z axis, the axis about which relative
rotation between the hinged objects can occur. Then if %ontroi is a scalar control
torque quantity, a positive torque term, TcontrolZhinge, is added into the angular
acceleration equation of one object. The negation of that term, --H Tcontro(Zhjnge
`01
is added into the other object's acceleration equation. Thus, leaving out the 5
pin hinge constraint equations, the object equations become
miri = mig + fhinge
Jiwi + `.`.`i x J1w1 = TihingeXhinge + T2hin9eYhinge + Tcontro?Zhinge
m2i:2 = m2g --H fhinge
J2w2 + W2 X J2w2 = TihingeXhinge --H T2hingeYhinge --H TconirojZhinge.
Once again, by viewing the control quantity as unknown, degrees of freedom
(here actually just one), are reintroduced into the equations describing object
behavior.
6.3 Basic control of motion
Once degrees of freedom have been introduced into the dynamics equations, two
basic approaches to motion control can be taken. The approaches are nearly
equivalent; either approach, however, might have advantages over the other for
use in particular applications.
During every iteration of the top-level simulation loop, a call to the inter-
face function run-control-phase is made, passing execution control to Newton's
control subsystem. The control subsystem then executes one or more control pro-
grams implementations of users' control strategies or algorithms.' These pro-
grams can access any relevant simulation state information by using the generic
`Control programs are best thought of as coroutines.
102
quantity access functions described in Chapter 3. Thus, control algorithms can
be based on values of object positions, velocities, accelerations or other simula-
tion parameters. For complex control schemes, programs can even request a set
of equation schemata describing the behavior of a mechanism.
Suppos?, then, that the addition of control quantities adds f degrees of
freedom to the basic dynamics equations, yielding an underdetermined set of
6n + ?m?=1 hi equations in 6n + ?1m=1 hi + f unknowns. The two approaches to
control both work by restoring this system to a completely determined set of
equations. In the first method control programs directly assign values to enough
of the unknowns to leave a determined system. The other method adds enough
control equations to the dynamics equations to again leave a determined set.
6.3.1 The direct control approach
In the first approach to control, the control programs are responsible for convert-
ing f of the unknowns into knowns by directly assigning values to them. During
each iteration of the top-level simulation loop the control programs can examine
the dynamics equations and the simulation state and then use the quantity access
functions to set the values of f of the unknowns. By doing so, the equations will
once again be reduced to a fully determined set of simultaneous linear equations.
Normally, for what we call forward dynamics control, the control programs
provide values for the f control quantities, leaving the object accelerations and
hinge reaction forces as unknowns. In the block example from above, the control
programs might assign, at each time instant, the value --Hrng to the control force
103
quantity. This would result in the object remaining fixed throughout the simu-
lation, since the applied control force would exactly balance the force of gravity.
Of course, the supplied values need not be constant over time. If the control
program supplied the value --Hmgt at each instant, where t represents the current
time, the block would move downward initially but eventually reverse direction
and accelerate upward indefinitely.
Control programs are not limited to supplying values for control force and
torque quantities; they can instead supply values for certain object accelera-
tions or even for hinge reaction forces. This is called inverse dynamic control.
Considering the block example once again, the control programs could set the
translational acceleration of the block to be zero at each time instant. The re-
sulting linear system of equations would then be solved to yield the control force
value, --Hmg, required to achieve this acceleration.
Note that a mixture of forward and inverse dynamic control is clearly sup-
ported. The control programs can choose to assign values to any f of the un-
knowns so long as a nonsingular set of 6n + ?m?=1 hi equations results.
6.3.2 The control equation approach
The second approach to control involves the use of control eqitations. Instead
of directly supplying values for particular unknowns, control programs define
equations to be added to the basic dynamics equations. By appropriately defining
f equations relating the unknowns of the dynamics equations, values for the
unknowns will become determined once again. For the block example, the control
104
equation method can be used to fix the block's position by adding either fcontrot =
--Hmg or = 0 to the basic equations.
In some situations, a combination of the direct control and control equations
approaches is useful. A control program can define a new abstract control quan-
tity and create an equation relating it to other quantities. The direct control
mechanism can then be used to assign values to the new quantity. For instance,
an experimenter might wish to control the acceleration of the center of mass of
a composite object. After defining a control quantity to represent the compos-
ite center of mass, the program can add a defining equation relating the new
quantity to the centers of mass of the component objects. Programs can then
use the direct control method to control the composite center of mass simply by
assigning values to the quantity.
6.3.3 Comparison of the two control mechanisms
The two control mechanisms are nearly equivalent. As in the block example,
the same influence on motion can often be achieved by either method. Both
approaches support pure forward and pure inverse dynamic control as well hybrid
schemes combining both forward and inverse dynamic control.
However, the control equation approach is somewhat stronger than the direct
control mechamsm since it directly affects the equation set integrated by the
dynamic analysis system. The control equations become an integral part of
the set of equations representing object behavior and are integrated as such. For
some purposes, such as simple animation, this is fine. However, for experimenting
105
with the design of real control systems, like those used in robot arms, the use
of control equations is too strong. Robot control systems have their 0Wll clocks
and a corresponding control cycle time. The controllers are only able to output
values influencing motion at time intervals that must be considered large with
respect to the size of the time steps taken by the integration module during
dynamic analysis. Thus realistic control systems can be only modeled using the
direct control method. Control programs can be given their own clocks that can
dictate, independently of the clock used in the integration module, when new
control values may be assigned to the controlled quantities.
6.3.4 Underdetermined equations
Both of the control mechanisms described above work by converting the under-
determined system of 6n + ?m?=? hi equations in 6n + ?:m.=i hi + f unknowns back
into a fully determined system. In many situations, requinng explicit elimination
of all f degrees of freedom puts an unreasonable burden on the experimenter. If,
using the direct method, fewer than f quantities are assigned values, or, using
the control equations method, fewer than f equations are added, the equations
remain underdetermined. In this case infinitely many solutions for the unknown
values exist, corresponding to many different simulated motions. For such situ-
ations, Newton allows control algorithms to guide the selection of a motion by
providing a quadratic cost function relating the unknowns. A solution to the
equations is then chosen that minimizes this cost function. Later sections of this
chapter provide examples of this technique.
106
6.4 Control and exceptional events
Control strategies, like other aspects of simulation, can be discontinuous. For
complex tasks, such as human walking, control algorithms generally consist of
a number of phases. During each phase different continuous motion control
strategies are active. To accommodate such complex control algorithms, Newton's
event handler can respond to exceptional events generated by the control module.
Control programs using the control equation method might want to alter the
set of equations that they had previously added to the basic dynamics equations.
They can do so by signaling an exceptional event to the event handler, which
then allows modification of the equation set. The event handler is able to treat
such control events just as it treats other kinds of events causing velocity or
(in this case) acceleration discontinuities. It coordinates handling of the control
events with the handling of other kinds of events. In addition, it can cause
restoration of the previous state and reintegration to a new time instant in order
to accurately isolate the time of event occurrence. Finally, it will make sure that
the integration module gets reinitialized to account for the motion discontinuity.
6.5 Control programming
The approaches to control described above simply provide two facilities for ex-
plicitly or implicitly specifying the values for the extra unknowns of the dynamics
equations. Clearly, a remaining difficulty involves the development of the control
programs that ultimately spedfy these values. Newton provides no special assis-
107
tance to control programmers beyond the facilities allowing simple uniform access
to quantities and other relevant simulation parameters. Control programs, then,
are wntten in Lisp and control algorithm phase information is usually maintained
using global variables. At the beginning of a simulation, the user places one or
more such control programs on a list maintained by the control interface function
run-control-phase. Each control program will be executed once during each
iteration of the simulation and may, as part of their execution, add or remove
other control programs from the list.
Sophisticated control tasks such as human walking will require complex con-
trol programs. The following section presents a number of examples of the kind
of control programming that has been done using Newton. The development of
languages and system components facilitating high-level programming of similar
control problems is an interesting area for future research.
6.6 Control examples
6.6.1 Simple low-level controllers
Complex control algorithms often depend on a combination of many simple lower-
level controllers. Users of Newton frequently define PD controllers2 and curve-
fitting controllers to control the "trajectory" of defined quantities. A small library
of such controllers has been developed.
In the biped walking example described in the next section, a number of the
simple controllers are used to control orientations, positions and joint angles. For
2?PD" stands for ?proportional, derivative?. For a discussion of PD, PID and similar con-
trollers see [5,16].
108
instance, quintic interpolation is used to plot the trajectory of the biped's heel,
and a PD controller is used to orient the foot before it struck the ground. Each
controller adds to the dynamics equations an equation that defines the second
derivative of the relevant quantity in terms of the first derivative and the quantity
itself.
The procedure in Figure 6.1 produces accelerations to move an object to
within 1% of a position p-desired within a given time delta-tirne.3The values
of p, v and a are quantity data structures representing the position, velocity, and
acceleration, respectively, of the controlled object. These data structures are
used by the add-control-equation function to create the appropriate equation.
Execution of the procedure position-with-PD causes a named equation to be
added to the set of dynamics equations. This equation will continue to affect the
motion of the object until it is explicitly removed by the control algorithm.
A list exhibiting the names of other controllers used by Newton control al-
gorithms is shown in Figure 6.2. Those with quintic in their name do quintic
interpolation to achieve the desired position and velocity in the desired time.
Quintic interpolation was chosen over cubic interpolation to eliminate jerk from
the beginning and end of the trajectory.
6.6.2 Control of biped walking
Newton has been used to experiment with the control of biped walking. The fol-
lowing describes an approach based largely on inverse dynamics. The approach
3As mentioned earlier, Newton users write control programs in Lisp. Pseudo?code is used here
for clarity.
109
procedure position-with-PD( constraint-name, object,
p-desired, delta-time )
var p, v, a: quantity
r: real
Begin
get-quantityC object, position )
get-quantityC object, velocity)
get-quantity( object, acceleration)
T = - delta-time / log( .01 )
add-control-equationC constraint-name,
a+2rv+W(p--Hp?desired)=O
end
Figure 6.1: PD controller used for positioning
position-with-PD( constraint-name, object, xd, At )
position--Hpoint-with-PD( constraint-name, object, point-on-object, xd, At )
orient-with-PD( constraint-name, object, ?d At )
set-angle-with-PD( constraint-name, joint, 0d At )
position-with--Hquintic( constraint-name, object, xd, Vd, At )
position-point-with-quintic( constraint-name, object, point-on-object, xd, vd, At )
orient-with-quintic( constraint-name, object, ?d? ?d? At )
set-angle-with-quintic( constraint-name, joint, 0d, 0d? At )
Figure 6.2: Low-level controllers
110
successfully demonstrates the Newton system's ability to support experimenta-
tion with motion control. However, because the method relies on inverse dynam-
ics and the assumption of a precisely known dynamic model of the controlled
object, it is not the kind of approach that would be used to design real control
algorithms for walking bipeds. Forward dynamics techniques that don't assume
perfect knowledge of object dynamics are currently being investigated.
The simple biped model
The simulated biped is composed of a torso, two legs with knee joints and two feet
with toe joints. This model was adapted from a description in [48] and is shown
in Figure 6.3. The hips and ankles are three-degree-of-freedom spherical joints,
while the knees and toes are one-degree-of-freedom revolute joints, making a total
of sixteen degrees of freedom. The biped is about six feet tall with moments
approximating those of a human being.
Eventual improvements to the model will incorporate joint limits and elastic
tendons. McMahon suggests that, during walking, energy is stored in stretched
tendons and is released when the stretched leg swings forward [48]. This idea
might be used to simplify the walking algorithm described in the next section.
To simplify this problem, impact and contact between the biped's feet and
the ground, were not handled by the normal impact and contact mechanisms.
Instead, the impacts and contacts are handled as part of the control problem by
adding external forces and torques to the feet. The torques and forces hold the
feet level with the ground until they are released by an explicit command from
111
fine 9.7
??ion
Figure 6.3: Simulated biped model
the control algorithm. The result is that the biped walks as if it was wearing
magnetic shoes on a steel floor.
The walking algorithm
An abbreviated version of the walking algorithm is shown in Figures 6.7, 6.8, and
6.9. The algorithm cycles through a set of six states: swing the right leg, land the
right foot, lift the left foot, swing the left leg, land the left foot, lift the right foot
and then repeat the cycle. In the swing phase, a quintic trajectory is plotted for
the swing foot with move-heel-to-target, while the stance leg is stiffened with
set-angle-with-PD and the foot is oriented for landing with orient-with-PD
(shown under START in Figure 6.9). In the landing phase, the leading leg is
stiffened as the foot nears the ground. Following this, the takeoffphase flexes the
trailing leg, causing the trailing foot to lift from the ground. Once the trailing
112
Constraint Name			dof			Constrained Item
TORSO-CONSTRAINT			3
L-KNEE-ANGLE			1
R-HEEL-TRAJ			3
R-FOOT-ORIENTATION			3
R-TOE-ANGLE			1
torso orientation in 3 dim
angle of revolute knee joint
heel acceleration in 3 dim
foot orientation in 3 dim
angle of revolute toe joint
Figure 6.4: Swing phase constraints
toe is bent to toe-joint-limit the flexing constraint is removed and the swing
phase begins for the trailing leg.
The largest number of constraints are applied during the swing phase, as
shown in Figure 6.4. Since the biped has sixteen degrees of freedom (dof) it
remains underconstrained at all times. A quadratic cost function is therefore
defined (in initialize of Figure 6.8) in order to fully determine the motion of
the biped. The cost function is a weighted sum of the translational and angular
accelerations, and of the difference between the torso translational acceleration
and some acceleration defined by a function F which tries to keep the torso half
way between the two feet.
We found that a cost function which minimizes trauslational and rotational
acceleration usually produces smooth motion. In the case of the simulated biped,
the cost function causes the constrained heel acceleration to be achieved by a
linear combination of small accelerations of many components of the body, rather
than a few large accelerations of those components which are near the heel. The
combination of many small accelerations yields more stable motion than large,
local accelerations.
113
?RLKI? t fo?t 1..d. fr
Figure 6.5: ?Valking cycle
Figure 6?5 shows ten frames in which the biped completes a full cycle of
the six phases described above. The full simulation consisted of twenty seconds
of straight-line walking on a flat surface and generated the statistics shown in
Figure 6.6.
6.6.3 Control of redundant manipulators
Redundant manipulators are robot arms that have more degrees of freedom than
necessary. Control algorithms can use these degrees of freedom to avoid obsta-
cles and singularities. Ltsing Newton's constraint equation capabilities we have
constructed a control algorithm which moves the end effector of a robot arm
along a specified trajectory. The algorithm uses the ideas presented in [41] in the
framework of the Newton system.
The redundant arm on which the algoritlim was tested is shown in Figure 6.10.
It has five degrees of freedom: the base rotates about a vertical axis and the four
114
,ft ?.. ooo?tfoo			tg?t ?..I ooo?t1.o
:;:t:I\I\.i\I\			____________
LI LO ?s ?..			?0 LI LO .1 0.0
.rt ?..I clootty (oh..d)			tI?t ?ooI .loo?t? c.h.?a)
(v)			(v)
Figure 6.6: Newton statistical output for walking simulation
115
time-in-air			=			0.5 s
stride			--H			0.5 m
direction			=			(1 0 0)
inside-ste?fraction			=			20 %
heei-Y-stnkespeed			=			--H0.05 m/s
heei-X-strik?speed			=			0.02 m/s
foot-strik?orientation			=			loo about (0 0 1)
tors?orientation			=			l0? about (0 0 1)
const
var phase: (start r-swing r-iand i-lift I-takeoff i-swing I-land r-lift r-takeoff)
procedure moveheel-t,>,target( constraint-name, foot, other-foot, hip, other-hip)
var target-x, target-v, hip?tohip: vector
begin
hi?t?hip = get-position( TORSO, hip ) - get-position( TORSO, other-hip)
target-x = get-position( other-foot, HEEL ) + stride x direction
+ inside-ste?fraction x hi?t?hip
target-v = heei-Y-strike-speed x (0 1 0) + heel-X-strike-speed x direction
position-point-with-quintic( constraint-name, foot, HEEL, target-x, target-v, time-in-air)
end
Figure 6.7: Definitions for the walking algorithm
procedure initialize
let F = K?(--H21(ri?100t + rr?joot) --H rtorso) +Kv(--H21(ri?j00t + ?r??oot) --H rtorso)
begin
quadra?ic-cost =
phase = START
end
+ ?r2 + 20(r"t orso --HF)2
Figure 6.8: Walking algorithm initialization procedure
116
procedure controller( time )
begin
case phase of
START:
phase = R-SWING
orient-with-PD( TORSO-CONSTRAINT, TORSO, tors?orientation, 2.0 s
move-hed-t?t&get( R-HEEL-TRAJ, R-HEEL, L-HEEL, R-HIP, L-HIP)
set-angle-with-PD( L-KNE?ANGLE, L-KNEE, 1750, 0.1 s)
orient-with-PD( R-FOOT-ORIENTATION, R-FOOT, foot-strike-orientation, time-in-air)
set-angle-with-PD( R-TO?ANGLE, R-TOE-J OINT, 00, time-in-air)
R-SWING:
if distance-to-target( R-FOOT ) < 0.01 m then
phase = R-LANDING
remove-constraint( R-HEEL-TRAJ )
set-angle-with-PD( R-KNEE-ANGLE, R-KNEE, 1750, 0.05 s)
R-LANDING:
if heeI-has?touched( R-FOOT) then
phase = L-TAKEOFF
remove-constraints( R-FOOT-ORIENTATION, R-TO?ANGLE, L-KNEE-ANGLE)
set-angle-with-PD( L-KNEE-ANGLE, L-KNEE, 1600, 0.1 s)
L-TAKEOFF:
if joint-angle( L-TOE-JOINT ) > 100 then
phase = L-SWING
remove-constraint( L-KNEE-ANGLE )
move-heel-t?target( L-HEEL-TRAJ, L-HEEL, R-HEEL, L-HIP, R-HIP)
orient-with-PD( L-FOOT-ORIENTATION, L-FOOT, foot-strike-orientation, time-in-air)
set-angle-with-PD( L-TOE-ANGLE, bL-TOE-JOINT, 1800, time-in-air)
Cases L-SWING, L-LANDING, and R-TAKEOFF
are anaiogo?s to the preceding three cases.
end
end
Figure 6.9: Abbreviated waiking aIgorithm
117
distal links rotate about parallel axes that lie in the horizontal plane. The func-
tion of the control algorithm is to move the end effector on a specified trajectory
while using the two extra degrees of freedom to keep the arm away from singu-
larities. We do not consider here the problem of obstacle avoidance.
In order to move the end effector, a time-varying acceleration constraint is
added to the system of equations. To stay out of singularities a quadratic cost
function is defined such that the solution of minimum cost is the solution that
steers away from singularities. We tried two naive choices for this cost function:
the sum of the squares of the joint torques and the sum of the squares of the
joint accelerations. While neither function caused the arm to stay away from
singularities, both resulted in far smoother movement than was the case when
no cost function was provided.
A better idea is found in [41], where it is stated that, for a given position of
the end effector, the manipulability measure of the robot arm, det(JJT) where
J is the Jacobian of the arm, is maximal when the joint angles of the arm are
equal.
Our algorithm, designed with this in mind, consists of two phases. The first
component takes a randomly oriented arm and moves it into a configuration of
equal joint angles without moving the end effector; the second phase moves the
end effector while maintaining equal joint angles.
To achieve equal joint angles individual PD controllers are created for each
of the three distal joints. Each controller yields a desired joint acceleration
that would bring the arm toward a configuration with equal joint angles. The
118
combination of these joint accelerations, however, might also cause movement of
the end effector.
To handle this, a quadratic cost function is chosen which is proportional to the
square of the difference between the allowed actual joint accelerations (i.e. those
accelerations that keep the end-effector in place) and the accelerations desired by
the PD controllers. In this manner, the desired accelerations are projected into
the space of accelerations that keep the end-effector fixed. Thus, the minimum-
cost solution to the motion equations will steer the arm into an configuration
with equal joint angles while maintaining the position of the end effector.
:Y???? t?:%
DcA;yyI			ffi4?tym
Figure 6.10: NIovement of a redundant robot arm
The equations bdow comprise the PD controller for the jth joint. These are
not added to the motion equations as are constraint equations; they are just
used by the control algorithm to calculate the desired joint acceleration. The
joint angle error and its derivatives are ??, ?j, and &i, while the individual joint
angles and their derivatives are ()j, ()j, and ()j. The average joint angle and velocity
119
= 0 --H
+ kv ?i + kp ?j = 0
are 0 and 0, respectively.
Shown below is the quadratic cost function which minimizes the error in the
joint angle.
5..
cost = Z(0i --H
i=3
Once the joint angles and joint velocities become equal, ?? will be zero. Under
this condition the cost function is minimized when 03 = 04 = 0?. Subsequent
movements of the end effector maintain equal joint angles (this is proven in [41]).
Simulations using this quadratic cost function resulted in the arm moving into
a configuration with equal joint angles. Subsequent arm movement maintained
equal joint angles and singularities were only encountered at the fully extended
configurations, where they are unavoidable.
Chapter 7
Conclusions and Continuing
Research
This chapter reviews the progress made by this work toward developing funda-
mentally more powerful simulation and design systems, and outlines the areas
that require further research.
7.1 Contributions
This thesis has proposed a number of basic architectural ideas enabling the de-
velopment of powerful design and simulation systems. In particular, the Newton
architecture and details of an initial implementation were presented. The Newton
architecture is novel in a number of ways. Unlike traditional simulation systems,
which concentrate mainly on integrating an unchanging set of equations, and
most current CAD systems, which concentrate on geometry specification but have
120
121
little in the way of simulation and analysis tools, Newton was designed provide
a level of automatic support that encourages the kind of experimentation nec-
essary for successful design. By using a model-based object representation and
fully integrating geometric modeling techniques it was possible to incorporate
into Newton a general mechanism to deal with events that cause discontinuities
in object behavior. Thus, Newton can automatically and incrementally modify
its internal description of mechanism behavior as relationships between objects
change due to events such as impacts, contact breakages, or changes in control
algorithm states. Such a facility greatly increases the power and flexibility of a
simulation system.
The primary contributions of the work in this thesis are
o+ the integration of geometric modeling techniques with traditional mechan-
ical simulation.
Previous systems have, for the most part, avoided geometric issues. They
place the burden of geometric reasoning on the experimenter by requiring
explicit user input of geometric quantity values such as volumes, inertias,
and potential impact points. The integration of geometric modelers with
other simulation components relieves users of such burdens and allows com-
prehensive automatic analysis mechanisms to be built.
o+ the development of significant automatic aids to mechanical analysis.
The analysis component of the Newton architecture relies on a general pur-
pose mechanism for handling except?ona1 events. Events such as impacts,
122
contact formation or breakage, and control algorithms state changes, were
handled in previous systems in ad hoc ways or not at all. Newton's event
handler allows a variety of such events to be handled in a natural and
uniform manner. As part of this, the dynamic analysis component of the
system works together with the event handler to provide incremental and
automatic reformulation of the sets of equations describing object behavior
in accordance with exceptional events.
o+ a demonstration of the flexibility provided by a system that uses a rich
structured representation to model the simulated objects. By explicitly or-
ganizing information about objects into a number of domain-specific mod-
els and having separate system components responsible for manipulating
information about particular domains, the system becomes flexible and ex-
tensible. Such a system, by allowing different users to emphasize different
aspects of object behavior, encourages experimentation that is necessary
in design and other tasks.
The goal of our work is quite broad - to develop a wide range of modeling
and simulation techniques that can be used together to fundamentally change
the coverage of design, analysis and simulation systems. The Newton architec-
ture provides the basis for such development. The initial implementation of the
architecture demonstrates the viability of our approach, but represents only a
first step; further research in a number of areas will be required to fully realize
the goal.
123
7.2 Continuing research
One of the examples mentioned in the introductory chapter involved testing de-
signs and control algorithms for anthropomorphic robot grippers such as the Sal-
isbury hand. While presenting no fundamental difficulty for the Newton system
architecture, realistic, accurate electronic prototyping of complex mechanisms
such as the hand will necessitate further development of a number of Newton's
components.
Geometry and object representations
The geometric modeling systems used in Newton have been quite simple, in-
cluding coverage of only a few primitive object geometries and limited boolean
operation and interference detection algorithms. Because of the interface between
Newton's other analysis components and the geometric modeler, substitution of
more sophisticated modelers is, for the most part, a simple matter of rewriting
the interface to work with a new modeler. However, available general purpose
geometric modeling systems often lack features required by Newton. Those that
do boolean operations efficiently are often not very robust or lack operations for
calculating mass properties. The robust systems are slow, and so on. A strong
geometric specification and editing facility is also necessary, and a well-developed
mechanism for referencing and reasoning about geometric features would be de-
sirable. Ultimately, we would like to see general purpose object editors that
support specification of the entire range of object properties: geometry, dynam-
ics, material properties, electrical properties, etc. By providing rich uniform
124
definition and representation facilities, such systems would provide a strong core
around which to base Newton's modeling components.
Mechanical analysis
Dynamics coverage should be extended to deal with nonrigid bodies, including
use of partial differential equation solvers and finite element techniques. Prelimi-
nary work on incorporating flexible bodies into Newton has been done at Purdue
University in an experiment that involved bouncing a ball off of a flexible beam.
The assumption of the ideality of simulated objects is another restrictive fea-
ture. Solving real-world robotics problems, for instance, requires that uncertainty
be accounted for. Object dimensions, locations, velocities, and material proper-
ties, are never known precisely. Furthermore, joints are not perfect; they usually
have some "play" in them. Dealing with such issues is an important, relatively
undeveloped, but active research area (see e.g. Donald[19]).
Furthermore, as alluded to in the thesis introduction, we would like future
systems to be able to handle an interesting variety of domains and thus be able
to handle in detail the modeling and analysis of mechanisms as complex as au-
tomobile engines or even spacecraft. Such modeling would involve combining
mechanical analysis with electrical, fluid, magnetic, and other analyses. N\'hile
the work of this thesis does not directly address the way in which such a goal
might be realized, Newton's use of a rich model-driven object representation to
combine previously independent analysis techniques suggests that such goals are
attainable.
125
Contact, impact, and friction
The methods described in Chapter 5 for handling impact, contact, and friction do
not meet the accuracy requirements of many engineering applications. Efficient
algorithms for the more sophisticated models studied by researchers in mechanics
and mechanical engineering are needed. Impact models need to better account
for the local deformations that occur even in nearly rigid objects. This thesis
described two particular schemes for determining the influence of impact-caused
impulses on other parts of a mechanism. The first assumed simultaneous oc-
currence and resolution of impulses on all contacts, while the second postulated
and solved a instantaneous sequence of impacts. Both methods have advantages
and model certain situations well, but more sophistication is clearly needed. In
reality, forces due to impacts do not propagate instantaneously, making the use
of impulses, or at least the simple techniques for impulse resolution described
here, unsuitable for some situations.
The contact algorithms that were presented also make unrealistic assump-
tions. Modeling a surface-surface contact as a finite set of points (the corners
of the convex hull of the intersection polygon in the polyhedral case) is suitable
maiuly for simulation of ideal frictionless rigid bodies. ?Vith the inclusion of fric-
tion, nonuniform surfaces, and curved objects, things are much more difficult.
Researchers in mechanics have done a great deal of work analyzing impact and
contact problems. Development of efficient computational models of such work
is clearly an area of important future research.
126
Definition and manipulation of equations and quantities
In order to encourage experimentation, users of a system like Newton must have
available the power to tailor the system to their own particular requirements. It
must be possible to easily define new quantities, constraints, hinges, and other
simulation entities. Equations for defining constraints or for use in control al-
gorithms should be easily manipulated. Toward this, powerful symbolic algebra
techniques and systems, on the order of MACSYMA, will probably be a necessary
component of general-purpose systems like Newton. As simulation complexity in-
creases, users of a simulator, as well as the system's analysis components, will
need to use these symbolic algebra techniques to reason about the simulated
processes.
Initial positioning of objects
Another problem not addressed by this thesis, but crucial to the ease of use of
large design systems, is that of positioning objects during design or before simu-
lation. It should be possible to develop algorithms and user interfaces that under-
stand natural descriptions of object configurations containing minimal amounts
of numerical data. Current systems, even if they don't require exhaustive numeri-
cal specification of object states, use a variety of numerical techniques (relaxation
or Newton-type methods), or, in limited cases, inverse kinematics techniques, to
provide the user with some help by requiring only a good approximation of the
complete intended state. Investigation into the use of information contained in
the description of an object, its topological structure and the kinematic relation-
127
ships between its parts, will lead to more powerful and easily used positioning
methods.
System robustness and efficiency
To be generally useful, simulation systems must be robust and efficient. Problems
of numerical instability make dealing with certain issues, contact in particular,
quite difficult. A partial solution is to use increasingly more sophisticated in-
tegration techniques. However, also warranted is an investigation into the use
of symbolic reasoning mechanisms for making decisions based upon numerical
tolerances. For instance, decisions about the nature of impact and contact rela-
tionships are based upon whether or not particular velocities and accelerations
are zero. Typically, simulation systems have a variety of such decision problems,
and accordingly, build up an a set of numerical tolerances for use in making the
decisions. The empirical manner in which the values are normally chosen leads
to a difficult-to-understand and not truly robust system. ?Vork examing the use
of symbolic reasoning techniques for solving these problems in the geometric do-
main is outlined by Hoffmann et al. [40]. Development of similar techniques for
mechanical simulation would be quite valuable.
The prototype implementation of the Newton system is slow. For interesting
simulations, the system generally requires between one hundred and one thousand
seconds to simulate one second of motion. The biped walking and redundant arm
examples of Chapter 6 each required a few hours of execution to produce a few
seconds of simulated motion. A good deal of Newton's slowness is due to a lack
128
of emphasis on speed in the initial implementation. Numerical computations
take up a large portion of simulation time. The matrix equations we solve are
very large but also very sparse. State-of-the-art solution techniques have not yet
been used in Newton and should yield significant speed-up for these numencal
computations.
In certain simulations, especially when contacts are made and broken very
quickly, the amount of time spent in equation reformulation and manipulation can
become significant. In many situations, however, constraint changes are relatively
rare. The support provided by Newton for efficient reformulation of equations
is then not really necessary since the time spent in numerical integration work
between motion discontinuities easily dominates. Under these circumstances,
dynamics techniques that use extra time to formulate a more efficiently solvable
set of equations would be beneficial, even if the formulation was not well suited
to incremental equation modification. Ideally, the simulation system itself will
be able to make such decisions and choose formulation methods based on the
structures of the simulated objects and perhaps some user-specified simulation
parameters.
Finally, it is clear that the next generations of simulation systems will benefit
from the development of parallel algorithms. Both in the geometry and in the
dynamics domains, there are many algorithms that are well suited to parallel
implementations. Ongoing work in this area will allow design and analysis sys-
tems combining algorithms from the various domains to work at near real-time
speeds.
Bibliography
[1] 5. Abhyankar, C. Bajaj, and C. Hoffmann. Computational algebraic ge-
ometry and geometry modeling. ACM SIGGRAPH-88, Notes for Tutorial
26, August 1988. Partial contents: Rudiments of algebraic geometry, Pa-
rameterization and implicitization of rational algebraic curves and surfaces,
Applying algebraic geometry to surface intersection evaluation.
[2]
A. Ambler. Robotics and solid modeling: A discussion of requirements
robotic applications put on solid modeling systems. In Robotics Research,
The Second International Symposium, pages 361--H367. MIT Press, 1985.
[3] W. W. Armstrong and M. W. Green. The dynamics of articulated rigid
bodies for purposes of animation. The Visual Computer, 1:231--H240, 1985.
[4] W. W. Armstrong and M. W. Green. Near-real-time control of human figure
models. IEEE Computer Graphics and Applications, 7(6):52--H61, June 1987.
[5] II. Asada and J.-J. E. Slotine. Robot Analysis and Control. John ?Vlley and
Sons, 1986.
[6] N. I. Badler and K. II. Manoochehri. Articulated figure positioning by
multiple constraints. IEEE Computer Graphics and Applications, 7(6):28--H
38, June 1987.
[7] D. E. Baraff. Analytical methods for dynamic simulation of non-penetrating
ngid bodies. Submitted to ACM SIGGRAPH `89.
[8] L. A. Barford. A Graphical, Language-based Editor for Generic Solid Models
Represented by Constraints. Ph.D. dissertation, Cornell University, March
1987.
[9] R. Barzel and A. H. Barr. Dynamic constraints. ACM SIGGRAPH-87,
Notes for Tutorial 17: Topics in Physically Based Modeling, July 1987.
129
130
[10] R. Barzel and A. 11. Barr. A modeling system based on dynamic constraints.
In Computer Graphics (SIGGRAPH 88), pages 179--H188. ACM, August 1988.
[11] J. Baumgarte. Stabilization of constraints and integrals of motion in dy-
namical systems. Computer Methods in Applied Mechanics and Engineering,
(1):1--H16, 1972.
[12] A. Borning, R. Duisberg, B. Freeman-Benson, A. Kramer, and M. `Vodf.
Constraint hierarchies. In Proceedings of the 1987 ACM Conference on
Object-oriented Programming Systems, Languages and Applications (0OP-
SLA 87), pages 48--H60. ACM, October 1987
[13] L. 5. Brotman and A. N. Netravali. Motion interpolation by optimal control.
In Computer Graphics (SIGGRAPH 88), pages 309--H315. ACM, August 1988.
[14] M. Chace. Methods and experience in computer aided design of large-
displacement mechanical systems. In E. Haug, editor, Computer Aided Anal-
ysis and Optimization of Mechanical System Dynamics[31], pages 233--H259.
Springer-Verlag, 1984. NATO ASI Series. Series F: Volume 9.
[15] M. Chace. Modeling of dynamic mechanical systems. Presented at the
CAD/CAM Robotics and Automation Institute and International Confer-
ence, Tuscon, Arizona, February 1985.
[16] J J. Craig. Introduction to Robotics: Mechanics and Control. Addison Wes-
ley, 1986.
[17] J. J. Craig. Issues in the design of off-line programming systems. In Robotics
Research, The Fourth International Symposium, pages 379--H389, 1988.
[18] R. Deyo, J. A. Biggs, and P. Doenges. Getting graphics in gear: Graphics
and dynamics in driving simulation. In Computer Graphics (SIGGRAPH
88), pages 317--H326. ACM, August 1988.
[19]
[20]
B. R. Donald. Error Detection and Recovery for Robot Motion Planning with
Uncertainty. Ph.D. dissertation, Niassachusetts Institute of Technology, July
1987.
P. E. Ehie and E. J. 11aug. A logical function method for dynamic and
design sensitivity analysis of mechanical systems with intermittent motion.
ASME Journal of Mechanical Design, 104:90--H100, January 1982.
[21] R. Featherstone. personal communication.
131
[22] R. Featherstone. The dynamics of rigid body systems with multiple concur-
rent contacts. In O. D. Faugeras and G. Giralt, editors, Robotics Research:
The Third International Symposium, pages 191--H196 The MIT Press, 1986.
[23] R. Featherstone. Robot Dynamics Algorithms. Kiuwer, 1987.
[24] W. C. Gear. Numerical linitial Value Problems in 9rdinary Differential Equa-
tions. Prentice Hall, Englewood Cliffs, NJ, 1971.
[25] W. C. Gear. Differential-algebraic equations. In E. 11aug, editor, Computer
Aided Analysis and Optimization of Mechanical System Dynamics [31], pages
323--H334. Springer-Verlag, 1984. NATO ASI Series. Series F: Volume 9.
[26]
B. Gilmore and D. Streit. A rule-based algorithm to predict the dynamic
behavior of mechanical part orienting operations. In Proceedings of the 1988
IEEE International Conference on Robotics and Automation, pages 1254--H
1259, April 1988.
[27] B. J. Gilmore. The Simulation of Mechanical Systems with a Changing Topol-
ogy. Ph.D. dissertation, Purdue University, August 1986. Mechanical Engi-
neering.
[28]
B. J. Gilmore and R. J. Cipra. The automatic determination of constraint
changes in dynamic multibody systems. Technical Report 86-DET-158,
ASME Paper, 1986.
[29] NI. Girard. Interactive design of 3-d computer-animated legged animal mo-
tion. IEEE Computer Graphics and Applications, 7(6):39--H51, June 1987
[30] J. K. Hahn. Realistic animation of rigid bodies. In Computer Graphics
(SIGGRAPH 88), pages 299--H308. ACM, August 1988.
[31] E. J. 11aug, editor. Computer Aided Analysis an(i Optimization of Mechani-
cal System Dynamics. NATO ASI Series. Series F: Computer and Systems
Sciences, Volume 9. Springer-Verlag, 1984.
[32]
E. J. 11aug. Elements and methods of computational dynamics. In E. 11aug,
editor, Computer Aided Analysis and Optimization of Mechanical System Dy-
namics[31], pages 3--H38. Springer-Veilag, 1984. N?TO ASI Series. Series F:
Volume 9.
[33] E. J. 11aug. Computer Aided Kinematics and Dynamics of Mechanical ys-
tems, Volume I: Basic Methods. Allyn and Bacon, 1989.
132
[34] E. J. Haug and G. M. Lance. Developments in dynamic system simulation
and design optimization in the center for computer aided design: 1980-1986.
technical report 87-2, University of Iowa, February 1987. The report is 893
pages long, containing approximately 50 of the group's papers.
[35]
E. J. 11aug, 5. C. Wu, and 5. M. Yang. Dynamics of mechanical systems
with coulomb friction, stiction, impact and constraint addition-deletion-i
theory. Mechanism and Machine Theory, 21(5):401--H406, 1986.
[36] E. J. 11aug, 5. C. Wu, and 5. M. Yang. Dynamics of mechanical systems
with coulomb friction, stiction, impact and constraint addition-deletion-ii
planar systems. Mechanism and Machine Theory, 21(5):407?16, 1986.
[37] E. J. 11aug, 5. C. Wu, and 5. M. Yang. Dynamics of mechanical systems
with coulomb friction, stiction, impact and constraint addition-deletion-iii
spatial systems. Mechanism and Machine Theory, 21(5):417--H425, 1986.
[38] C. 11offmann and J. 11opcroft. Automatic surface generation in computer
aided design. The Visual Computer 1, pages 92--H100, 1985.
[39] C. M. 11offmann and J. E. 11opcroft. Simulation of physical systems from
geometric models. IEEE Journal of Robotics and Automation, RA-3(3):194--H
206, June 1987.
[40]
[41]
[42]
[43]
C. M. 11offmann, J. E. 11opcroft, and M. 5. Karasick. Robust set operations
on polyhedral solids. technical report 87-875, Cornell University, October
1987.
P. 11su, J. 11auser, and 5. Sastry. Dynamic control of redundant manipula-
tors. In Proceedings of the 1988 IEEE International Conference on Robotics
and Automation, pages 183--H187, April 1988.
P. M. Isaacs. Controlling computer generated motion with dynamics kine-
matics, and behavior functions. Masters dissertation, Cornell University,
January 1988.
P. M. Isaacs and NI. F. Cohen. Controlling dynamic simulation with kine-
matic constraints, behavior constraints and inverse dynamics. In Computer
Graphics (SIGGRAPH 87), pages 215--H224. ACM, July 1987.
[44] M. 5. Karasick. On the Representation and Manipulation of Rigid Solids.
Ph.D. dissertation, McGill University, 1989.
133
[45] J. U. Korein. A Geometric Investigation of Reach. Ph.D. dissertation, Uni-
versity of Pennsylvania, 1984.
[46] P. L5tstedt. Mechanical systems of rigid bodies subject to unilateral con-
straints. SIAM Journal of Applied Mathematics, 42(2):281--H296, 1982.
[47] P. L5tstedt. Numerical simulation of time-dependent contact and ftiction
problems in rigid body mechanics. SIAM Journal of Scientific and Statistical
Computing, 5(2):370--H393, 1984.
[48] T. A. McMahon. Mechanics of locomotion. The International Journal of
Robotics Research, 3(2):428, 1984.
[49]
M. Moore and J. Wilhelms. Collision detection and response for computer
animation. In Computer Graphics (SIGGRAPH 88), pages 289--H298. ACM,
August 1988.
[50] M. E. Mortenson. Geometric Modeling. John ?Viley and Sons, New York,
1985.
[51]
[52]
[53]
J. I. Neimark and N. A. Fufaev. Dynamic of Nonholonomic Systems. Trans-
lations of Mathematical Monographs, Volume 33. American Mathematical
Society, Providence, RI, 1972.
N. Orlandea, M. A. Chace, and D. A. Calahan. A sparsity-oriented approach
to the dynamic analysis and design of mechanical systems - parts i and
u. Journal of Engineering for Industry, Transactions of the ASME, 99:733--H
779,780--H784, 1977
G. P. Ostermeyer. Numerical integration of systems with unilateral con-
straints. In E. Haug, editor, Computer Aided Analysis and Optimization
of Mechanical System Dynamics[31], pages 415--H418. Springer-Verlag, 1984.
NATO ASI Series. Series F: Volume 9.
[54] R. 5. Palmer. Computational Complexity of Motion and Stability of Polygons.
Ph.D. dissertation, Cornell University, January 1987
[55] A. Paoluzzi. Integration constraints in parametric design of physical objects.
technical report 87--H804, Cornell University, January 1987.
[56] B. Paul. Analytical dynamics of mechanisms - a computer-oriented
overview. Mechanism and Machine Theory, 10:481--H507, 1975.
134
[57] B. Paul. Kinematics and Dynamics of Planar Machinery. Prentice-Hall,
Englewood Cliffs, NJ, 1979.
[58] J. C. Platt and A. 11. Barr. Constraint methods for flexible models. In
Computer Graphics (SIGGRAPH 88), pages 179--H188. ACM, August 1988.
[59] C. H. Suh and C. W. Radcliffe. Kinematics and Mechanisms Design. John
Wiley and Sons, 1978.
[60]
[61]
D. Terzopoulos and K. Fleischer. Modeling inelastic deformation: Viscoelas-
ticity, plasticity, fracture. In Computer Graphics (SIGGRAPH 88), pages
269--H278. ACM, August 1988.
D. Terzopoulos, J. Platt, A. Barr, and K. Fleischer. Elastically deformable
models. In Computer Graphics (SIGGRAPH 87), pages 205--H214. ACM, July
1987.
[62] D. Terzopoulos and A. Witkin. Physicall-based models with rigid and de-
formable components. In Proc. Graphics Interface `88, pages 146--H154, 1988.
[63] D. Terzopoulos, A. Witkin, and M. Kass. Energy constraints on deformable
models: Recovering shape and non-rigid motion. In Proc. AAAF87, pages
755--H760, 1987
[64]
[65]
[66]
M. Vanderploeg and G. M. Lance. Controlled dynamic systems modeling. In
E. Haug, editor, Computer Aided Analysis and Optimization of Mechanical
System Dynamics[31], pages 401--H414. Springer-Verlag, 1984. NATO ASI
Series. Series F: Volume 9.
5. E. Walter. Collision Detection of Polygonal Objects at Interactive Rates.
Ph.D. dissertation, Cornell University, August 1985. Mechanical Engineer-
ing.
R. A. Wehage and E. J. 11aug. Dynamic analysis of mechanical systems
with intermittent motion. ASME Journal of Mechanical Design, 104:778--H
784, October 1982.
[67] R. A. Wehage and E. J. 11aug. Generalized coordinate partitioning for
dimension reduction in analysis of constrained dynamic systems. ASME
Journal of Mechanical Design, 104:247--H255, January 1982.
[68] J. Wilhelms. Using dynamic analysis for realistic animation of articulated
figures. IEEE Computer Graphics and Applications, 7(6):12--H27 1987
135
[69] J. ?Vilhelms and B. A. Barsky. Using dynamic analysis to animate articu-
lated bodies such as humans and robots. In Proceedings: Graphics Interface
`85, pages 97--H104, 1985.
[70] J Wilhelms, D. Forsey, and P. Hanrahan. Manikin: Dynamic analysis for
articulated body manipulation. Technical Report UCSC-CRL-87-2, Univer-
sity of California at Santa Cruz, April 1987.
[71]
A. Witkin, K. Fleischer, and A. Barr. Energy constraints on parameterized
models. In Computer Graphics (SIGGRAPH 87), pages 225--H232, July, 1987
ACM.
[72] A. Witkin and M. Kass. Spacetime constraints. In Computer Graphics
(SIGGRAPH 88), pages 159--H168. ACM, August 1988.
[73] J. ?Vittenburg. Dynamics of Systems of Rigid Bodies. B. G. Teubner
Stuttgart, 1977
[74]
J. Wittenburg and U. Wolz. The program mesa verde for robot dynamics
simulations. In 0. D. Faugeras and G. Giralt, editors, Robotics Research,
The Third International Symposium, pages 197--H204. MIT Press, 1986.
[75] D. Zeltzer. Motion control techniques for figure animation. IEEE Computer
Graphics and Applications, 2(9):53--H59, 1982.
[76]
D. Zeltzer. Knowledge-based animation.
tion and Perception, pages 187--H192, Toronto,
GRAP11/SIGART Interdisciplinary Workshop.
In MOTION: Representa-
Canada, April 1983. SIG-
[77] D. Zeltzer. Towards an integated view of 3-d animation. The Visual Com-
puter, 1:245--H259, June 1985.