Pan Gun Park Emmanuel Witrant
email: emmanuel.witrant@ee.kth.se.

Mikael Johansson

Carlo Fischione

Karl Henrik Johansson

Predictive control over wireless multi-hop networks

HAL is

Introduction

There are major advantages in terms of increased productivity and reduced installation costs in the use of wireless communication technology in industrial control systems. Several wireless technologies are currently reaching the plant and workshop floor addressing a range of applications. A multi-hop wireless network is one such technology that offers great flexibility, but also challenges when it comes to predictability. These networks are used in emerging applications in industrial automation and dedicated algorithms, such as randomized multi-hop routing protocols [START_REF] Bonivento | Randomized protocol stack for ubiquitous networks in indoor environment[END_REF], were recently shown to be useful in these applications. From the perspective of designing control systems utilizing wireless multi-hop networks, one of the most significant characteristics is the need of hiding the system complexity through suitable abstraction of the networks and their nodes. One simple but important way of doing this is by interpreting the network as a communication channel with time-varying delay.

Experimental results on telerobotics over TCP/IP networks illustrate the importance of varying transmission delays in the control loop [START_REF] Lelevé | Telerobotics over IP networks: towards a low-level real-time architecture[END_REF], where the authors combine a generic estimation/prediction control scheme with a buffering of the received packets, which compensates for the delays jitter. While the buffer may decrease the transient performances of the closed-loop system [START_REF] Witrant | On the use of state predictors in networked control systems[END_REF], this work motivates the use of estimation/prediction schemes to achieve robust and efficient control over networks. Another interesting experimental result is provided in [START_REF] Ploplys | Closedloop control over wireles networks[END_REF], where a peer-to-peer wireless network (802.11b) and a UDP protocol are used to transmit the sensors and actuators signals.

The main contribution of the paper is a new scheme for control over wireless multi-hop networks. The characteristics of the routing protocol together with lowerlayer network mechanisms give rise to a delay process with high variance and stepwise changing mean, which is particularly difficult to model. Here we use a loworder model proposed in [START_REF] Jacobsson | Some modeling and estimation issues in traffic control of heterogeneous networks[END_REF] and propose a new predictive control scheme with a delay estimator. The estimator is based on a Kalman filter with a change detection algorithm [START_REF] Jacobsson | Some modeling and estimation issues in traffic control of heterogeneous networks[END_REF]. It is able to track the delay mean changes but efficiently attenuate the high frequency jitter. Compared to similar works such as [START_REF] Seuret | Implementation of an internet-controlled system under variable delays[END_REF], where an output feedback scheme uses the maximal value of the delay and GPS synchronization for remote control over Internet, we explicitly compensate for the time-varying delay effect. Network data from an experimental setup are used to illustrate the efficiency of the approach.

The outline of the paper is as follows. Section II defines the considered problem of control over a wireless multi-hop network. The wireless network and its randomized routing protocol is described in Section III. The time-delay estimator is presented in Section IV. Section V presents the controller and its delay compensator. Experimental results are given in Section VI followed by the conclusions in Section VII.

Problem description

The problem considered is depicted in Fig. 1, where a linear system is remotely controlled over a wireless multi-hop network. We suppose that round-trip time or end-to-end delay measurements (i.e., between the sensors and controller, which also implies a proper clock synchronization) are available. A predictive control law uses the filtered measurements to compensate for variations induced by the multi-hop network. More precisely, we consider the remote control of linear systems that write as

ẋ(t) = Ax(t) + Bv(t), x(0) = x 0 (1) y(t) = Cx(t) (2)
where x ∈ R n is the internal state, v ∈ R is the system input, y ∈ R m is the system output, and A, B, C are matrices of appropriate dimensions. The pairs (A, B) and (A,C) are assumed to be controllable and observable, respectively, but no assumption is made on the stability of A. The delays introduced by the network dynamics are defined as

• τ 1 (t): delay between the sensors and the controller,

• τ 2 (t): delay between the controller and the actuator.

The actuator is event-driven and uses the last received signal. Packet losses influence the actuation as

w(t k ) = u(t k -τ 2 (t k )) if packet k is received w(t k ) = w(t k-1 ) else ( 3
) where w(•) is the signal received by the actuator from the network, u(•) is the controller output and v(t) = w(t k ) ∀t ∈ [t k ,t k+1 ). The chronology of the signals is indicated in the transmitted data in order to reconstruct their history on the receiver site. The receivers consider only the most recent information, meaning that a packet is also considered as lost if a more recent one is received first.

Wireless multi-hop networks

In many scenarios of relevant interest, wireless multi-hop networks do not have fixed communication paths, but the end-to-end path followed by packets happens according to a dynamic selection of hops [START_REF] Greuen | Adaptive sleep discipline for energy conservation and robustness in dense sensor networks[END_REF][START_REF] Xu | Geographyinformed energy conservation for ad hoc routing[END_REF]. A major technique to ensure power savings and longer network lifetime is to turn off a node whenever its presence is not strictly required for the correct operation of the network [START_REF] Xu | Geographyinformed energy conservation for ad hoc routing[END_REF][START_REF] Chen | Span: An energy-efficient coordination algorithm for topology maintenance in ad hoc wireless networks[END_REF][START_REF] Bonivento | Randomized protocol stack for ubiquitous networks in indoor environment[END_REF], which means that the network topology is changing randomly.

Consider the multi-hop network reported in Fig. 1, where the transmitter to receiver path is composed of a time-varying number of hops h(t). For each node, denoted by i = 1 . . . h(t), we can consider the following sources of delay:

• the time α i to wait before sending a data packet, which is typically a random variable;

• the time F to forward a data packet once the connection with the next node is set up, which is constant and contains the propagation and transmission delays;

• possibly the time β i induced by an Automatic Repeat reQuest (ARQ) mechanism, which retransmits a packet for a predefined number of time if no acknowledgement is received.

The end-to-end delay can then be expressed as follows:

τ(t) = h(t)F + h(t) ∑ i=1 (α i + β i ) (4)
If one looks at the outcomes of (4), there are three possible behaviors. The low frequency delays are due to changes in the number of hops h(t). The middle frequency ones depend on the selection of the next awaken hop for a fixed number of clusters (related to µ c,i in the example considered), and the high frequency behavior is related to the transmissions between nodes.

Time-delay estimator

A critical component of the proposed control strategy is the delay estimator, whose output is used to set the time-varying horizon of the predictive controller. Examining the real delay traces from our network, we see that the delay exhibits a large variability and persistent changes in the mean. Although these changes can be traced back to a particular feature of our protocol (discussed in the previous section) similar effects have been reported for a wide range of network quantities and network technologies [START_REF] Kim | Mobile network estimation[END_REF]. Since the controlleractuator latency is not immediately available to the controller, but has to be deduced from the (delayed information about the) controller-actuator round-trip time, we will look for a delay estimator that allows to smooth out the fast variations while being able to react quickly to the changes in mean.

To this end, we propose to use a combination of Kalman filter and CUMSUM change-detection (cf., [START_REF] Jacobsson | Some modeling and estimation issues in traffic control of heterogeneous networks[END_REF][START_REF] Gustafsson | Adaptive Filtering and Change Detection[END_REF]). The underlying signal model is

τ k+1 = τ k + v k r k = cτ k + w k
where τ k is the one-way latency, r k is the measured time, c is the number of delays included in r k (i.e., 1 if the end-to-end delay of the channel is measured and 2 if it is the round trip time), w k is a white-noise sequence with covariance R k and v k = 0 with probability 1q, ν with probability q, where Cov(ν k ) = q -1 Q k models the abrupt changes in the mean. The associated Kalman filter updates are given by

τk+1 = τk + K k (r k -c τk ) K k = cP k-1 R k + c 2 P k-1 P k = P k-1 - c 2 P 2 k-1 R k + c 2 P k-1 + Q k
The tracking ability of the filter is proportional to the variance of the process noise Q k , but increasing Q k also decreases the smoothing effect that we are looking for.

To improve performance, we run a CUMSUM detector in parallel with the Kalman filter. The detector looks at the prediction errors ε k = r kc τk , sums up the negative and positive prediction errors

g + k = max(g k-1 + ε k -κ, 0) g - k = max(g k-1 -ε k -κ, 0)
and monitors if g + k or g - k exceeds a threshold g thr . Here, the parameter κ is a design parameter that specifies the negative drift (see [START_REF] Gustafsson | Adaptive Filtering and Change Detection[END_REF] for design guidelines). At these threshold crossings, Q k is set to a large value, g + and g -are reset to zero, τk is set to c -1 r k and the Kalman filter iterations resumed. The combined Kalman-filter and CUMSUM test allows to achieve both fast tracking of the changes in the mean and good suppression of the high-frequency variations.

Predictive control approach

The control approach proposed here is based on the state predictor with a time-varying horizon first proposed in [START_REF] Olbrot | On controllability of linear systems with time delays in the control[END_REF]. The main advantage of such method compared to other delay compensation strategies [START_REF] Richard | Time delay systems: An overview of some recent advances and open problems[END_REF] is to allow for a finite spectrum assignment on systems with time-varying delays. An overview on the use of state predictors in networked control systems is provided in [START_REF] Witrant | On the use of state predictors in networked control systems[END_REF], where it is compared with more classical control approaches.

Ideal delay compensation

We use the results established in [START_REF] Witrant | Remote stabilization via time-varying communication network delays: Application to TCP networks[END_REF] and subsequent works, where the delay induced by the network dynamics is explicitly taken into account in the control setup. For simplicity and to focus on the network effect on the closed-loop system, we suppose that the full delayed state 1 x(t -τ 1 ) is available to set the control law. When there is no packet loss and the delay is fully known (and can be predicted), the main theorem in [START_REF] Witrant | Remote stabilization via time-varying communication network delays: Application to TCP networks[END_REF] can be extended to the case where both communication channels experience a delay as follows.

Theorem 5.1 Consider the system ẋ(t) = Ax(t) + Bu(t -τ 2 ), x(0) = x 0 with (A, B) a controllable pair. Assume that the network dynamics is such that the following holds for τ 1 and τ 2

A1) τ 1,2 ∈ C (R + , [0, τ max ]) ∀t, A2) τ1,2 < 1 ∀t ≥ t 0 .
Then the feedback control law

u(t) = -Ke Aδ e Aτ 1 x(t -τ 1 ) +e At t+δ t-τ 1 e -Aθ Bu(θ -τ 2 (θ ))dθ , δ (t) = - λ 1 -dτ 2 (ζ )/dζ δ + dτ 2 (ζ )/dζ + λ τ 2 (ζ ) 1 -dτ 2 (ζ )/dζ ,
with ζ = t +δ , λ a positive constant, δ (0) = δ 0 , ensures that the system trajectories x(t) converge exponentially to zero.

If A1-A2 hold (which correspond to a lossless network description), then the closed-loop behavior can be inferred from the time-shifted dynamics dx/dζ = (A -BK)x(ζ ) and its spectrum coincides with the spectrum of the matrix A -(1 + δ )/( τ(t + δ )) e -Aδ BK, where ( •) denotes the differentiation with respect to time t.

Considering the worst-case delay properties, the feedback gain has to be chosen such that the matrix

A - 1 1 -τmax e -Aτ max BK (5)
is Hurwitz.

Use of the delay estimator

We now consider that the control law is set using the delay estimator described in section 4. Both delays τ 1 and τ 2 are then replaced by τ in the control setup description introduced in Theorem 5.1. This is motivated by the fact that the predictor architecture in used to compensate for the time-variations of the latency (i.e., the slow and medium dynamics of the delay), which are estimated with the CUMSUM Kalman filter. It is reasonable to suppose that the average values of the delays in both communication channels are close and provided by the estimator. The control law is then given as

u(t) = -Ke Aδ e A τ x(t -τ 1 ) +e At t+δ t- τ e -Aθ Bu(θ -τ(θ ))dθ . (6)
The predictor horizon is computed using the actual values of the estimated delay and delay variation directly instead of the predicted ones (at time ζ ). We then approximate δ with δ , which dynamics is given by

δ (t) = τ(t) + λ ( τ(t) -δ (t)) 1 -τ(t) . ( 7
)

Control algorithm

Special care has to be taken in the computation of the control law ( 6)- [START_REF] Greuen | Adaptive sleep discipline for energy conservation and robustness in dense sensor networks[END_REF], especially concerning the integral term used in the predictor's part. Indeed, the discretization of the integral may introduce some numerical instabilities [START_REF] Mondié | Finite spectrum assignment of unstable time-delay systems with a safe implementation[END_REF]. The effect of these instabilities is further increased if the communication channel experiences packet losses and bandwidth limitation. A solution to this problem is to include the discretized integral in a dynamic control law, as described in this section.

Considering the integral term I (t)

. = e At t+ δ t-τ e -Aθ Bu(θ -τ(θ ))dθ computed at the time instant denoted by k, we have that

I k = e A τk τk + δk 0 f (µ)dµ with f (µ) = e -Aµ Bu(min(µ + t k -τk -τ(min(µ + t k - τk ,t k )),t k ))
to ensure causality. I k is discretized using the results proposed in [START_REF] Mondié | Finite spectrum assignment of unstable time-delay systems with a safe implementation[END_REF], with a trapezoidal rule and included into the dynamic controller

               ż(t) = -az(t) + bKe A( δk + τk ) x(t -τ 1 ) + t s 2 n k ∑ i=1 f ((i -1)t s ) + f (it s ) , u k = z(t k ) (8)
where n k . = ( δk + τk )/t s , a and b are some positive design parameters, t s is the sampling time set by the controller clock and n k is the number of steps used to discretize the integral, which is time-varying. Note that we chose the controller sampling time as the integrand step, since it sets the buffering of the past estimated delays and controller outputs. An upper bound τthr is set on the maximum estimated delay variation and the dynamics of the predictor horizon is set with

δ (t) = τk + λ ( τk -δ (t)) 1 -τk , (9)
where τk = ( τkτk-1 )/t s if τ < τthr , and τk = τthr else. The issues considered in this section are summarized with the following result.

Result 5.1 Using the delay estimation provided by the CUMSUM Kalman filter described in section 4, the proposed predictive control approach is set with the dynamic controller [START_REF] Xu | Geographyinformed energy conservation for ad hoc routing[END_REF], where the predictor horizon is set by the dynamics [START_REF] Chen | Span: An energy-efficient coordination algorithm for topology maintenance in ad hoc wireless networks[END_REF] and the controller gain is such that the closed-loop matrix ( 5) is Hurwitz.

Experimental and simulation results

In this section, the network setup is presented along with some experimental results, which are used to estimate the efficiency of the proposed control setup using appropriate simulation tools. More details on the experimental setup and protocol specifications can be found in [START_REF] Park | Randomized protocol design in wireless automation[END_REF]. We use the Randomized Protocol (RP) proposed in [START_REF] Bonivento | Randomized protocol stack for ubiquitous networks in indoor environment[END_REF] and briefly discussed in Section 3.

Network and experimental measurements

Our experiment is based on the RP using T-motes Sky nodes equipped with the radio controller Chipcon CC2420 (2.4 GHz) and satisfying the communication standard IEEE 802.15.4 [START_REF]Tmote Sky Data Sheet[END_REF]. The source is located at 30 m from the destination node connected to the PC and generates 20 pckts/s. Eighteen nodes are placed at regular intervals (1.58 m) along a straight line without obstacles. Surrounding objects are static, with minimal timevarying changes in the wireless channel due to multipath fading effects. Fig. 2 shows the average end-to-end delays and Packet Reception Rate (PRR). The number of clusters is set depending on the channel condition and to reduce the transmission power. Both the average endto-end delays and PRR depend on the number of hops. To save energy, the RP increases the number of clusters in case of high PRR, and vice versa.

Physical system and control setup

Consider the T-shape ECP inverted pendulum presented in [START_REF] Witrant | Remote stabilization via time-varying communication network delays: Application to TCP networks[END_REF], which is an unstable non-minimum phase system with system matrices

A =     0 1 0 0 0.80 0 12.56 0 0 0 0 1 -2.42 0 -8.33 0     , B =     0 4.57 0 0.38     .
The eigenvalues of A are λ 1,2 = 0.74 ± 2.08i, λ 3,4 = -0.74 ± 2.08i. The choice of an unstable open-loop plant emphasizes the efficiency of the controller to compensate for the transmission perturbations and the performance limitations.

Experimental delay measurements τ 1 and τ 2 are obtained by splitting the initial experimental data into two sets, one every other sample being attributed to each set (the channels are not symmetric but belong to the same network). Only the last received packets are considered and the packet losses are handled as described in (3). The bandwidth limitation is set by zero-order holds with a sampling time of 0.05s. The controller includes the time-delay estimator (based on end-to-end delay measurements), a reference trajectory block that sets the control objective, and the predictive control setup. The CUMSUM Kalman filter is set with R = 10 000, Q = 1, ν = 0.01 and g thr = 35. The reference trajectory, for the pendulum angle x(3), is a filtered square signal. The predictive control is set with λ = 10, a = b = 100 (the introduced controller dynamics is equivalent to a lowpass filter) and K = [0.48 0.36 2.20 0.57] sets the poles of ( 5) to (-3 + 0.1i; -3 -0.1i; -4; -4).

Simulation results

Fig. 3 present the actual and estimated time-delays (top), the reference and the system output (middle), and the control output (bottom). The transmissions are very noisy, in terms of jitter and packet losses (up to 60% of the signal). This figure illustrates the control scheme efficiency and the proper handling of the high frequency noise induced by the numerical integration problem. The effect of the number of clusters on the system response is particularly significant when a change in the reference sign occurs when there are three clusters. Indeed Fig. 3 depicts a situation where the bad quality of the communication channel introduces some oscillations, which are reduced as soon as the number of clusters is set back to one. Comparing the different simulation results and extra tests carried with ARQ, it appears that the system oscillations are induced by the low PRR rather than high end-to-end delay (i.e., comparing the different 3-clusters cases). This is due to the fact that the control setup explicitly compensate the delay effects while the robustness with respect to packet losses is not specifically addressed in the feedback design.

These simulations also show that an increase in the number of clusters has few effects on the system response when it is established in a steady state. A control oriented multi-hop wireless network protocol could efficiently use this information to minimize the emission energy (increase the number of clusters) according to some measurements of the physical state while ensuring satisfying performances.

Conclusions

We addressed the problem of remote control over a multi-hop wireless network. The constraints of such a network on the communication channel, such as timevarying delays, jitter, packet losses and bandwidth limitations are considered. The time-delay induced by the network is estimated thanks to a CUMSUM Kalman filter, that provides for a fast tracking of the mean and good suppression of the high frequency variations of the delay (jitter). This particular Kalman filter is used to design a dynamic predictive control law that compensates for the estimated part of the delay explicitly. The causality of this controller is set according to the estimator and specific considerations on safe implementation are detailed. A necessary condition on the assigned closed-loop spectrum is provided but the particular poles location is not specified, which provides for a degree of freedom for optimal control approaches and future works. Based on experimental delay and packet loss data, some simulations illustrate the proposed control scheme.

Figure 1 .

 1 Figure 1. Overview of the system and control setup.

Figure 2 .

 2 Figure 2. Average end-to-end network-induced delays (top) and packet reception rate

Figure 3 .

 3 Figure 3. Delays, system output and control signal without ARQ.

The time dependency of the variables δ (t), τ 1 (t) and τ

(t) are omitted in the notations from this point.

This work has been supported in part by the European Commission through the projects SOCRADES and HYCON, the Swedish Agency for Innovation Systems via the Nordite project VISA, the Swedish Research Council and the Swedish Foundation for Strategic Research.


Ronen Barzel

Alan H Barr

A Modeling System Based On Dynamic Constraints

Keywords: Modeling, Dynamics, Constraints, Simulation CR categories: 1.3.5-Computational Geometry and Object Modeling; 1.3.7-Three-Dimensional Graphics and Realism 1

We present "dynamic constraints," a physically-based technique for constraint-based control of computer graphics models. Using dynamic constraints, we build objects by specifying geometric constraints; the models assemble themselves as the elements move to satisfy the con straints. The individual elements are rigid bodies which act in ac cordance with the rules of physics, and can thus exhibit physically realistic behavior. To implement the constraints, a set of "constraint forces" is found, which causes the bodies to act in accordance with the constraints; finding these "constraint forces" is an inverse dynamics problem.

Introduction

Some of the most natural and graceful motion in computer anima tion has been achieved recently by simulating the physical behavior of objects. But physical simulation has not yet become the standard technique for modeling and animation, because of several limitations:

• Simulations are hard to implement: Typically, a special-purpose program is written to simulate the behavior of a given computer graphics model; the overhead for making new models is large.

• Simulations are hard to control: If "innate" behavior is pro grammed into models, it becomes hard to make the models do exactly what we want; the behavior is often determined indirectly by non-intuitive or non-orthogonal parameters.

• Simulations are slow: Physical simulation can be computationally intensive.

The goal of this work is to develop a modeling system in which it is easy to build and animate physically-based computer graphics models. To this end, our modeling approach is based on four features:

• Genemlity: A model is built from a collection of primitive physically-based elements.

• Geometric Constmints: A model is constructed by applying con straints to the objects, starting from an initial confi guration of the primitive elements. A model is also positioned and animated through constraints.

• Newtonian Mechanics: Each primitive element is a rigid body whose motion is due to the effects of inertia and forces and torques acting on the body. Many of the forces and torques are externally applied; other forces and torques, however, are derived from the geometric constraints.

• Equivalence of Modeling and Animation: The temporal behavior of physically based objects is bound up in the model itself.

To implement the constraints, we solve an inverse dynamics prob lem: given constraints on the behavior of the model, the problem is to determine the forces which result in an example of the constrained behavior we desire. Thus, we convert each constraint into a "constraint force"; as the model animates, the constr' aint forces are continuously computed, to continuously maintain the constraints.

Sec. 2 of this paper presents the modeling system, and provides implementation notes. Sec. 3 discusses the inverse dynamics problem.

Sec. 4 presents the technique for setting up and solving a "constraint force" equation. The simulation of Newtonian mechanics, derivations for various examples of constraints, and miscellaneous mathematical details are found in the Appendices.

Related Work [Witkin,Fleischer,and Barr 87] uses "energy'' constraints to assem ble 3D models, for changing the shape of parametrically-defined primitive objects.

This work is not concerned with dynamic me chanical simulation of models. [Platt and Barr 88] uses augmented lagrangian constraints in the physical simulation of flexible ob jects. [Isaacs and Cohen 87] does physical simulation of rigid bod ies, for the special case of linked systems without closed kinematic loops.

They share our emphasis on ease of modeling, and also use an inverse-dynamics formulation to control the models' behavior. [Wilhelms and Barsky 85] utilizes physically based modeling, but has a reduced emphasis on control.

2

The Modeling System

Modeling with the "Dynamic Constraints" system consists of instan tiating primitive bodies, connecting and controlling them with con straints, and influencing their behavior by explicitly applying external forces. The modeling system thus has three libraries:

• Primitive bodies: A collection of rigid• bodies, such as spheres, rods, torii, and more complex shapes, that are the component elements of models. The modeler specifies the body density, as well as specific parameters such as the. length and radius of a rod. Each body type defines the quantities needed for physical simulation, such as the rotational inertia tensor for that body type.

• External applied forces: Forces that the modeler can introduce into the model, including gravity, springs, and damping forces.

Each force has parameters specific to the force, such as damping coefficients or spring constants. • Constraints: Various types of geometric constraints, such as �'point-tO-nail" or "orientation" ate described below. To build a model, we make instances of objects, calculate the forces, and run the simulation. We also specify "timelines'' of events to take place, such as creating or removing instances, turning constraints or explicit forces on or off, or otherwise adjusting parameters.

2.1

Newtonian Mechanics Fig. 1 illustrates one of the simplest examples of a body obeying New ton's laws. This behavior is easy to simulate; Appendix A describes the general newtonian simulation procedure. All primi tive bodies in our system exhibit physically realistic be havior, in the sense that they respond correctly to forces and torques.

Constraints

We show some examples of constraints supported by our modeling sys tem.

• "Point-to-Nail" constraint (see Fig. 2): This fixes a point on a body to a user-specified location in space. The body may sw ivel and swing about the constrained point, but the constrained point may not move.

• "Point-to-Point" constraint (see Fig. 3): This forms a joint be tween two bodies. The bodies may move about freely, as long as the two constrained points stay in contact.

• "Point-to-Path" constraint (see Fig. 4a): We can require a point on an object to follow an arbitrary user-specified path; this al lows us to animate models by using standard kinematic keyframe techniques.

• "Orientation" constraint (see Fig. 4b): A constraint to align ob jects by rotating them.

• Other constraints (not illustrated): Other constraints include "point-on-line," which restricts a point to lie on a given line, and "sph ere-to-sphere," which requires two spheres to touch, but lets them slide along each other.

We can easily add new types of geometric constraints to our con straint library, by defining the constraint "deviation" function and de riving various required quantities, as described in Appendix C. The only restriction is that the "deviation" function be twi ce-differentiable (as is discussed in the appendix).

2.3

Constraint Forces

When we have built a model using dynamic constraints, the model is held together by constraint forces, as illustrated in Fig. 5. Thus the constraint forces are analogs of the internal forces which hold the parts of compound objects together. Constraint forces also assemble the models, pulling the com ponents into the proper configurations.

Thus constraint forces represent forces which could be used to assemble real-world ob jects.

For example, figures in the appendices show frames from animations demonstrating the self-assembly of space structures [B arr , Von Herzen, Barzel, and Snyder 87]. b) "Orientation" con straint. We rotate the rod to make its axis parallel to the slot axis.

Implementatiort

We describe here the high-level program structure; details of si mulat ing Newtonian mechanics are given in Appendix A, the procedure to calculate the constraint forces is given in Sec. 4.4.

Our modeling system is implemented in Common Lisp on Symbolics Lisp Machines, using Symbolics' object-oriented "Flavors" mechanism. The fundamental object classes we have defined are:

• rigid-body: a primitive body in the model. This class defines the functions and state variables needed for the dynamics calculation (see Appendix A), including a list of forces and torques acting on the body. There are subclasses for each type of body in our library; each subclass provides type-specific information, such as the rotational inertia tensor.

• control-point: a point on a body, or in space. A point on a body contains a reference to the body, as well as the position of the point in body-coordinates. A point in space defines its position, which can be constant, or a function of time. Fo rces and con straints are typically created by specifying the control-points at which they act.

• force: a force being applied to a body. Each force contains a reference to a control point at which it is applied. There are subclasses for each type of force in our library; each subclass provides a function that computes the force.

• constraint: any type of dynamic constraint. There are subclasses for each type of constraint in our library. Each subclass provides shows forces on the upper rod, (b) shows forces on the lower rod: F9 is gravity pulling down on rods. F1 is the "point-to-nail" constraint force on the upper rod, holding it at the naiL F2 is the "point-to-point" constraint force on the lower rod, holding it to the upper rod; -F2 is the reaction force on the upper rod.

the quantities needed to determine the constraint force (described in Appendix C). Each constraint also keeps references to the bodies being constrained, and associates the appropriate forces with the bodies. All objects handle a "draw" message, which displays the object in its current state. For debugging a model, we send the "draw" mes sage to all objects, including forces, control points, and constraints; for producing an animation, we send the "draw" message only to bodies.

Some examples of subclasses are:

• rod: a subclass of rigid-body. This class provides the values spe cifi c to rods, e.g. the rotational inertia tensor (see Appendix A).

The class also associates two control points with each rod, named "endl" and "end2", at the ends of the rod, and provides functions to access them. -

• nail: a control point fixed at a location in space.

• point-to-nail: a subclass of constraint. Provides the functions which calculate the terms needed for a "point-to-nail" constraint (see Sec. 4, Example 1). The addition of new types of bodies, forces, or constraints to the system merely requires the creation of an appropriate new subclass.

Currently, the user-interface is via the lisp environment; for example, the pendulum of Fig. 3 could be built via the series of commands:

; create bodies and control points (make rod "upper-rod") (make rod "lower-rod") (make nail "nail" 0 0 100) ; specify constraints (connect (endl "upper-rod") "nail") (connect (end2 "upper-rod") (end1 "lower-rod"))

; add e,;ternal forces (gravity-on) ; apply gravity to each body To animate a model once the instances are made, we simply iterate these steps:

• Simulate until end of frame (Appendix A).

• Send "draw" message to objects. The implementation makes heavy use of a home-grown pack age of numerical routines, which include linear-system solvers, dif ferential equation solvers, and the like; some useful references are [Press et.

3

Inverse Dynamics

If we are given the forces which act on a collection of objects, we can easily solve the fonoard dynamics problem-that of determining the objects' behavior-as described in Appendix A. However, to meet constraints, we must solve the reverse problem: Given a partial de scription of the desired behavior, we must determine forces which will yield an appropriate behavior. This inverse dynamics problem, sum marized in Fig. 6, consists of two parts: (a) fi nding forces to meet a constraint, and (b) fi nding forces to maintain a constraint. Meeting A Constraint Fig. 7 shows a constraint force being used satisfy to part (a) of the in verse dynamics problem, that of moving the objects to meet an initially unmet constraint. Notice that this part of the problem is actually very loosely specified: How quickly should the constraint be met? Along what path should the object move? For our solution, as we shall see in S ec. 4.2, the "deviation" of the constrained point decays exponentially, with a user-specified time constant.

Maintaining A Dynamic Constraint Fig. 8 shows a dynamic constraint force adapting to satisfy part (b) of the inverse dynamics problem, that of keeping a constraint met de spite motion and other forces. There is typically a unique solution, in which the dynamic constraint forces provide the internal forces that hold together an object. In this section, we present the technique for computing the con straint forces. The presentation is in several sections:

1. Definition of several mathematical quantities associated with each constraint.

2. Construction of a "constraint-force" equation for a dynamic con straint; if the constraint force is chosen such that this equation holds, the constrained objects will behave in accordance with the dynamic constraint.

3. Grouping the constraint-force equations for multiple constraints into a single multidimensional constraint-force equation for all the dynamic constraints.

4. Setting up and solving the dynamic constraint-force equation. A' A d X 3 matrix -coiTesponding to the net torque on body i;

Definitions

we have one such for each body in the constraint.

{3

The part of Dl2l independent ofF and T; ad-dimensional vector.

I

The number ot <legrees ot lree<lom m the constramt force.

G•

A 3 X J matrix. The constraint force on body i is given by G; Fe

We have one such matrix for each body in the constTaint.

H• A 3 X f matrix. The constraint torque on body i is given by Hi Fe

We have one such matrix for each body in the constraint. We also define f c

The unknown "constraint force," an /-dimensional vector.

ij:

The net external applied force on the i-th body

TE

The net external applied torque on the i-th body

(1)

Note that strictly speaking, f . is not necessarily a f�: n: .>e, but rather is a quantity that determines the constraint force(= G' Fe) and constraint torque(= Hi f e) on the constrained bodies; colloquially, however, we refer to Fe as the "constraint force." The vectors f � and Tj, are due to the external forces, such as gravity, that act on the i-th body body.

The net force (torque) on a body is the sum of the net external force (torque) and the constraint forces (torques):

F' = ( L G}Fe;) + fl: 'fi = (c•• • f";H�FeJ + 'f1
(2)

CODI'I;rain'h j

In the above equations, we label terms for the ith body with super script i's, and the jth constraint with subscript j's.

See Example 1 for a description of terms needed for a "point to nail" constraint.

4.2

Constraint-force Equation For each constraint, we describe the desired behavior of the constraint "deviation" by linearly combining D, jjC1), and i)C2l:

flC2l(Y,F,T,t) + �jj(l)(Y,t) + 1 2 D(Y,t) = 0, t �to ( 3 )
-r r

This equation i� equivalent to the differential equation in Fig. 12; its solution brings D down to 0, then holds it at 0.1 We will substitute for D ( 2 ) , D ( l ) and D in Eqn. 3 to produce a linear system of equations in which we solve for f •. If we continually adjust the force so that Eqn. 3

is met, we will be solving the inverse dynamics problem.

Thus, we expand D(2) in Eqn. 3, using the definitions in Fig. d.

e.

f. g.

a.

t =to t=to+0.3r t=to +0.6r t=to+0.9r f=fo+2.0r t ;<:to+ 6.07" Example 1: "Point-to-Nail" Constraint. We choose jj to be the vector from the constrained point, at Xp, to the "nail", at Xa (see Fig. 9).

We thus ha.

G = i FE = Fg Substituting into Eqn . 5 gives: 1 - 1 rt 2 1 - - -Fe+ -r9+-ii+ 2 (X-Xo) = 0 m m -r T
We easily solve for the constraint force: (5 )

Notice that, to meet this one constraint, we must take into account the effect of all the constraint forces.

See Example 2 for a sample derivation of the constraint forces.

Multiple Constraints

Each constraint resu lts in a version of Eqn. 5; with several constraints, we have a set of simultaneous equations which must be solved. We duplicate Eqn. 5, for each constraint in the model:

L ( L (ri G)+ AiHj))Fe 1 ) =0,{ for all constraints k (6)
where we label terms for the kth constraint with subscript k's. Writ ing this system of equations more compactly, as a multidimensional vector equation, we have the constraint force equation for the model:

+!i + _1_ jj(l) + 1 jj k 'TJc k � k (7)
Fig. 13 illustrates collecting individual constraint equations into the multidimensional vector equation. Notice that each element of M is a matrix, and each element of Fe and of B is a vector.

4.4

Solving the Constraint-Force Equation Fig. 15 outlines the procedure to set up Eqn. 7, as well as solve it and compute the net force and torque on each body.

In step 3 of Fig. 15, we call a standard linear-system solver to solve Eqn. 7. There are many well-known methods for solving linear sys tems (see [Press et. al. 88,Ralston and Rabinowitz 78]). We note some characteristics of M that should be taken into account when choosing a solution method:

• M is typically sparse. The [k, j]th entry in M is non-0 only if some body is influenced by both constraint k and constraint j. Typically, most of the elements are zero.

• M is not necessarily square. A constraint may have d oJ. f; for example, the "orientation" constraint (Appendix C.3) has d = 1 and f = 3, yielding a matrix which is "wider" than it is "tall." • M may be singular, implying that Eqn. 7 is overconstrained or • underconstrained? We most often use singular-value decomposition (SVD) to solve Eqn. 7, because it robustly handles singularity and near-singularity, as well as non-square systems. However, SVD does not take advantage of sparseness, and is a relatively slow technique.

Underconstrained Equations

Constraint-force equation Eqn. 7 will sometimes be underconstrained, thus having many solutions. This can occur, for example, when there are several constraints acting on a single body; it may be possible to vary some of the individual constraint forces without affecting the net torque or force on the object. An example is shown in Fig. 14a, in which the pair of forces labeled "V" yield the same net force ( = 2Fv) and torque(= 0) as the pair labeled "W''.

There is no difficulty caused by having many solutions to Eqn. 7; we could use any solution, since they will all yield satisfactory behavior. We might wish to use the solution which is smallest in magnitude, to avoid numerical difficulties; SVD yields this solution.

Overconstrained Equations eats In Fig. 14b, the user has specifi ed constraints which can not be met; there is no "correct" constraint force to be applied. In Fig. 14c, the specifi ed constraints can be met, but not by moving the constrained point in a straight line; however, Eqn. 3 requires that the point move in straight line if the constrained point is initially at rest.3

For overconstrained systems, using the least-squares solution for the constraint forces typically yields "reasonable" behavior -the ob ject typically assumes some intermediate confi guration, for the case of Fig. 14b, or moves along the feasible path, for the case of Fig. 14c. SVD computes the least-squares solution for overconstrained systems.

Summary

We have developed a modeling system featuring constraint-based con trol of rigid bodies. The bodies' behavior is determined by simulation of Newtonian mechanics. We compute dynamic "constraint forces" to apply to the bodies such that they behave in accordance with user specifi ed geometric constraints; the computation of these forces is an inverse dynamics problem.

The modeling system supports various types of geometric con straints, such as "point-to-nail" and "point-to-point." The modeler builds objects by using constraints to connect primitive components; the constraint forces cause the components to assemble themselves into the model, and ensure that the model stays assembled as it animates.

2Unfortuna.tely, we have some overloading of the word "constrain": 11overcon strained" and "underconstrained" refer to the linear system of equations Eqn. 7, rather than to the constraints themselves.

3 A solution to the probl<m1 of unreal izable paths is to use scalar constraint mea sures (d = 1). For example:, the "point-to-nail" constraint could be redefined so that D is the distance from the point to the nail, rather than the vector separating the point and the nail. We have developed a technique to compute the constraint forces by setting up and solving a "constraint-force equation." The constraint force equation is a multidimensional linear equation of the form MFe + B = 0, where Fe is the collection of unknown constraint forces.

Each constraint is described by a "deviation" measureD, such that D = 0 when the constraint is met. D must be a twice differentiable function of the positions and orientations of the constrained bodies.

Appendix C derives D for several types of constraints.

6

Future Work Further work we are interested in pursuing includes:

• Expanding the constraint library. Deriving new constraint "devi ation" functions, as described in Appendix C.

• Interactive graphical modeler.

• Object Intersection. Development of non-interpenetration con straints

• Flexible bodies.

Incorporation of fl exible-body simulation with dynamic constraint control [Platt and Barr 88, Te rzopoulos, Platt, Fleischer, and Barr 87].

• Special-case models Direct implementation of the equations of motion for common objects, such as the linked systems of • Constraints on velocity or acceleration.

We are also looking forward to using the dynamic constraints mod eling system as a tool in other research areas, such as molecular biology [Lengyel 87] and robotics. !!._Q = !.;; Q dt 2 We then define R to be an a uxiliary variable, which is computed from Q as discussed in [Shoemake 8 5).

• Rotational Inertia Ten&or I: I determines the xota.tiona.l behavioi of a body.• For a rigid body, Ibod> is constant. Note also that in Fig. 16 <t. A discussion of the cha.ra.cterstics of I is beyond the scope of this pe.perj see [Fox 67,Goldstein 83]. [Lien and Kajiya 84] gives an algorithm to compute I for arbitrary nonconvex polyhedra.

'ii' E b X a b'T ii' Eii' X b b'Tb = 0 B.2
Behavior of a Point Consider a point "P" which is fixed relative to a rigid body (Fig. 17). We define bbody to be the vector from the center-of-mass of the body to P, ex pressed in the body's home coordinates; since the point is fixed, bbodll is cons tant. We would like to derive expressions for the position, velocity, and acceleration of P.

We will need to know the derivative of I-1 • Remember that since the body is rigid, x;.;..

We have again substituted w'R for AR.

Finally, we can express the positio'if, velocity, and acceleration of point P in terms of the state of the body and the net force and torque on the body: For each type of constraint, we must derive expressions for the various quan tities defined in Fig. 11. The steps we follow are: 1. Choose a simple "deviation" measure D. D is a function of the po sitions (X) and orientations (R) of the constrained bodies, and may optionally depend on t. 2. Differentiate D, to derive j)('l(y, t); Substitute ;; and w' R for the f. X and f,R terms which will arise (see Fig. 11).

X p X+b

3. Differentiate ag ain, to derive i5<2l(Y, I, T, t). Replace f,p and f,l terms with I and t', thus giving rise to the linear dependence of j)C•l on the fo rces and torques. Define the d X 3 matrices r, A, and the d-vector if 4. Choose where to apply the constraint forces needed to meet the con straint. Most often, we apply a vector force to a fixed location of the constrained body; in this case, we have f = 3 degrees of freedom.

5. Use steps 2 and 3 to derive G and H for each body. These convert the f values in the "const raint force" Fe into the actual forces and torques on the bodies.

Often, some of the quantities r, A, G, ;�.nd H, which are nominally ma�ri ces , turn out to be scalar. Scalars can be handled as a special case in the implementation, or scaled identity matrices can be used.

We give examples of the constraint derivations for the constraints illus trated in Sec. 2.2. C.2 "Point-To-Point" Constraint This constraint is met if the two constrained points "P1 " and "P•" are at the same location (Fig. 18a). We thus define D to be the vector separating the two points. The derivation proceeds analogously to that of the "point ii Notice that A is a 1 x 3 matrix, and iJ i s a scalar.

We apply an arbitrary pure torque, Fe , to the body, and no fo rce. We therefore have f = 3, and G =0; H= I Notice that this constraint is "non-square" -we are applying 3 degrees of freedom to affect a scalar constraint "deviation."

Figure 1 :

 1 Figure 1: All primitive bodies obey Newton's laws. e.g. (a) A ball released in gra vity falls; (b) A ball thrown in gravity moves in an arc.

Figure 2 :

 2 Figure 2: "Point-to-nail" constraint. A user creates a pendulum by fix ing an endpoint of a rod at some location in space. The constraint causes the rod to "fly" into pl ace, assembling the pendulum. See videotape [Cal tech '87 Demo Reel].

Figure 3 :

 3 Figure 3: A "point-to-point" constraint. A users adds a second rod to the pendulum of Fig. 2, to create a compound pendulum. See videotape [Caitech '87 Demo Reel].

Figure 4 :

 4 Figure 4: Other constraint examples. (a) "Point-to-path" constraint. This constraint pulls objects along user-specified paths. (b) "Orientation" con straint. We rotate the rod to make its axis parallel to the slot axis.

Figure 5 :

 5 Figure 5: Constraint forces holding together compound pendulum of Fig. 3d. The constraint forces model the internal forces of a real-world pendulum. (a) shows forces on the upper rod, (b) shows forces on the lower rod: F9 is gravity pulling down on rods. F1 is the "point-to-nail" constraint force on the upper

  al. 88,Golub and Van Loan 83,Ralston and Rabinowitz 78, Boyce and Deprima 77]. We also have embedded into lisp an exten sion to "Einstein Summation Notation" for mathematical expressions [Barr 83,Misner, Thorne and Wheeler 73]; this makes it quite simple to create lisp functions by merely typing in the mathematical formulae using the sa.me notation with which we derive them.

  Such that constraint stays met, despit-e motion and other forces. gravity 2\1 {,] II \stays iJ fixed � ____ .,

Figure 6 :

 6 Figure 6: The inverse dynamics problem for dynamic constraints.

  Assembling a model: a.. Given a constraint to be met: nail •---� .f1 b. Introduce for . ce : constrai nt �-force ,;?nail c. Force pulls object: (((•� d. Force l!llows object: e. Constraint is met: F igure 7: Meeting a. Constraint. The constraint force pulls the ball towards the nail, then brings the ball to rest a.t the nail.

Figure 8 :

 8 Figure 8: Maintaining a Constraint. (a-d) The constraint force adapts to hold constraint even as object moves and other forces act on it. The con straint force pulls up, to counteract gravity, and sideways, to keep the pen dulum's inertia from flinging it sideways of£ the nail.

Figure 9 :

 9 Figure 9: "Point-to-nail" constraint. The "deviation" measure is D(Y) = XP(Y)-X0, the constraint force is Fe, and the constraint torque is b>< Fe. 4 Calculating Constraint Forces Note: We suggest skimming Appendix A, to become fa miliar with our notation and formulation of rigid-body me chanics, before reading this sectio n.

  Fig.11gives the definition of the quantities which must be supplied for each constraint. The derivations of these quantities for various types of constraints are given in Appendix C.

  jj A measure of "deviation" for the constraint: D(Y, t) = 0 <==> constraint is met D is a d-dimensiona.l vector. d The number of dimensions of D. .Dl'l The rate of change of .D: t)(l>(Y(t), t):: �D(Y(t)) zj(l) is a d-dimensional vector. l)l2 ) The acceleration of .D: rP D<2>(Y(t),F(t), T(t),t):: dt • D(Y(t),t) D(2) will depend linearly on F and T; thus: n<2J = L <r'<Ylfl'' + A'<Y)T') + iJ<Yl 3 matrix corresponding to the net force on body i i we have one 5Uch for each body in the constraint.

Figure 11 :

 11 Figure 11: Quantities associated with a. dynamic constraint. Fe Is the un known "constraint force" for the constraint. Y, :F, and T are the state, net force, and net torque in the model, as defined in Eqn. 9 (Appendix A). See discussion in Sec. 4.1. Derivations of these quantities for various constraints are given in Appendix C

1

  Analytically, if .Do f;. 0 the solution to the equation in Fig.12asymptotically approaches 0, but doesn't ever rea.ch D = 0. Nllltlerically, however, we soon reach D = 0 within error tolerances.

Figure 10 :

 10 Figure 10: Constraint-force calculation for a "point-to-nail" constraint (details in Sec. 4.2): Constraint force has components opposing gravity ( -F0), opposing motion (-F.), and pulling towards nail (-Fv). (a) User specifies constraint at center of mass. (b) Constraint force initially pulls towards nail. (c-e) Once ball is moving towards nail, constraint force turns around. (f) Constraint force slows ball. (g) Steady -state: Fv = Fv = 0, Fe = -F •. Net force on ball is 0; by Newton's first law, the ball remains at rest.

  Figure 12: D evolving over time. We ha.ve picked a second-order differential equation to describe D as a. function of time. The solution yields the behav ior required by the inverse dynamics problem-the constraint "deviation" decays down to 0, assembling the model, then remains at 0, maintaining the constraint. The• ra.te of assembly is controlled by the time constant r. See [Boyce a.nd Deprima. 77] for a discussion of second-order differential eq ua. tions.

  F9 -; mv-;:2 m(X-Xo)Thus we see the constraint force ha.s three components: One opposing the force of gra.vity, one opposing the ball's velocity, and one pulling the ball towards the nail. Fig.10illustrates the constraint force adapting to pull the ball to the nail, and bring it to rest. Once the ball is at rest at the na.il, we have X-Xo = 0 and v = 0; so the constraint force becomes Fe = -F9, yielding a. net force on the ball of F = Fe + F9 = 0 con•tre.iat.• j bodi•• i + 2:: cr; ffJ,; + A;Te)

  (r�G} + AtH;) bsdl.e• i F e; =Fe; B. = L (r�f1 + A �Te) bodie• i

  nail

Figure 13 :Figure 14 :

 1314 Figure 13: Multiple constraints: Each constraint contributes one line to the equation. The collection of constraints together yields a set of simultaneous linear equations, expressible a.s a linear matrix equation. n 'It,

[

  Isaacs and Cohen 87,Armstrong and Green 85]. Decrea;; ing the number of constraints in the model speeds up the constraint-force calculation.

namic

  Analysis To Animate Articulated Bodies Such As Humans and Robots, GraphicsInterface, 1985.

Fig. 16

 16 Fig. 16 summarizes the equations of motion of a rigid body. A full discussion of rigid-body dynamics can be found in [Fox 67,Goldstein 83].A.l Notes On The Equations Of Motion• The Orientation Matrix R: R transforms tensors from body coordinates to world coordinates (see Fig.17). As we numerically

  17

  ). As we numerically integrate R, numerical noise tends to cause R to drift away from a pure rotation, yielding noticeable skewing. This can be alleviated by using a feedback technique, as in [Barr 83). Alternatively, we can rep resent the orientation as a qua.ternion Q (see[Shoemake 85] for an introduction to quaternions). The equation of motion for Q is (see[Misner, Thorne and W heeler 73]):

Figure 17 :

 17 Figure 17: A rigid body. Orientation matrix R transforms vectors from body coordinates to world coordinates.we use 1-1 rather than I. I;;;,' dy can be precomputed for each body; we convert to world coordinates using R:

  b

  ;r' = (f,R)Ib;, �.R T + RI;,;.. <_fR T) (10) = w'RI;" ,;..a T + RI;,;..a w' T = w'I-1 + x-1w' T We have substituted w'R. for f,R, according to the equations of motion (Appendix A). We will also need the derivative of w: w =l-1 Z f. w = ( ;ji; I-1 ) f + r-1 (j,l) = Ct;I-1 )L + I-1T = w'I-' l + I-1w' T l + 1-11' (ll) = w'w + r-1(w'T l + T) =1-'(l xw+T) We have substituted T for ft.. l according to the equations of motion. We can now differentiate b: b = Rbbod.l', and f,'b = f,(R.bbcdy) = ( f, R)bood• = w'Rbbody =w•b =wxb f,.. 'b =( f.w) x �+wx(f,b) � = (f,w) x b+w x (w x b) = W1 (l x w + 1')) x h + w x (w x b) = (b'T x-1)T + (b'T x-'(l x w) +w x (w x b) )

  Figure 18:

  (a.) "Point-to-point" constraint. We apply equal-and opposites to the two constrained points, to ca.use the deviation fu nction D(Y) = Xp,(Y) -Xp1 (Y) to go to 0. (b) "Orientation" constraint. We rotate the body to cause the deviation fu nction D(Y) = b • N -I to go to 0. The constraint torque is given directly by Fe ; there is no force due to the constraint. C.l "Point-To-Path" Constraint This constraint is met at a time t if the constrained point "P" is on the path, at gpath(t); it is the same as the "point-to-nail" constraint (Example 1 in Sec. 4) but with a nail that moves. Thus several terms have a dependency on gp at h , D(Y, t) = XP(Y)gpath(t) iJ(li(Y, t) = vp (Y) -f,.XP•1 h (t) iJ<�>cy, f, t , t) = ap (Y, f, f) -.£,-.XP•th (t) i3 =b•Tr-'(L xw)+wx(wxi> ) -.£,-.fpath(t)

  keeping with Newton's third law, the two bodies must exert equal and opposite fo rces on each other. We apply an arbitrary force, Fe , to one of the constrained points, and the negation of that forc e, -Fe, to the other. We thus have f = 3, and defineG1 =1; G2 = -1; H1 =b� H2 = -b; C.3 "Orientation" ConstraintThis constraint is met if a specified unit vector b fixed in the body lines up with a unit vector N fi xed in the world (see Fig.18b). We could defin e D to be the angle between the vectors; it easier, however, if we define D to be 0 when the cosine of the angle (i.e. the dot-product of the vectors) is 1. Thus we have d = 1, and:

Figure 19 :

 19 Figure 19: Linking Chain between Two To wers. The chain swings naturally aft er assembly.

Figure 20 :Figure 22 :

 2022 Figure 20: Linking Chain between Two Towe rs , continued

Figure 23 :

 23 Figure 23: Spa.Ce Station Assembly, continued .

Figure 24 :

 24 Figure24: Cardbouse Assembly. We use "point-to-point," and "pt�to, plane" constraints to assemble a cardhouse.

Figure 25 :

 25 Figure 25: Cardhouse Assembly, completed.

Figure 26 :

 26 Figure 26: Pandora.'s Chain. The chain links are instructed to connect togetherand hook to a. trap door. To rsion springs keep the links roughly perpendicnla.r to each other. Gravity a.nd viscous damping are applied to all bodies. The chains and trap door swing naturally once they are assembled.

Acknowledgements

Thanks to: Jed Lengyel, for being a "guinea pig'' user of the modeling system, and for putting together the video system that was used in making animations; John Snyder, for the rendering of the animations; John Platt, for discussions and for numerical software; and Brian Von Herzen, for the motivat ion and modeling for the Harwood space-station self-assembly animations and for the polygonal playing card models. Figure 21: Linking Chain between Two Towers, continued

where: