SimMechanics Release Notes
The SimMechanics 2.1.1 Release Notes on page 1-1 describe the changes in this product since Version 2.1.
Note Version 2.1.1 is available as part of the R13 SP1 release.
If you are upgrading from a version earlier than Version 2.1.1, you should also see: SimMechanics 2.1 Release Notes on page 2-1 SimMechanics 2.0.1 Release Notes on page 3-1 SimMechanics 2.0 Release Notes on page 4-1
Printing the Release Notes
If you would like to print the Release Notes, you can link to a PDF version.

Contents

SimMechanics 2.1.1 Release Notes
Major Bug Fixes. 1-2 Upgrading from an Earlier Release. 1-3 SimMechanics Demos. 1-3
SimMechanics 2.1 Release Notes
New Features. Computer-Aided Design Translator for SolidWorks. import_physmod Command for Generating Models. Special Documentation and Examples for CAD Translator.

2-2 2-2 2-2 2-2

Major Bug Fixes. 2-4 Upgrading from an Earlier Release. 2-5 Installing the CAD Translator for SolidWorks. 2-5 Known Software and Documentation Limitations. SolidWorks Versions Compatibility for CAD Translator. Uninstalling and Reinstalling the CAD Translator. Platform Limitations. SimMechanics Users Guide Not Updated.

2-6 2-6 2-6 2-6 2-6

SimMechanics 2.0.1 Release Notes
Major Bug Fixes. 3-2 Code Compilation with lcc Compiler. 3-2
SimMechanics 2.0 Release Notes
New Features. Code Generation with Real-Time Workshop. Accelerator Mode and Discrete Events. Expanded Handle Graphics Visualization Features. Variable Body Mass and Inertia Tensor. Vectorized Joint Initial Condition Actuation. Damped Linear Force Element Blocks. New Mechanical Branching Bar Block.
4-2 4-2 4-2 4-2 4-2 4-2 4-2 4-3
Major Bug Fixes. 4-4 Upgrading from an Earlier Release. 4-5 Reconfiguring Virtual Reality Visualization. 4-5 Known Software and Documentation Limitations. Virtual Reality Toolbox Viewers. Configurable Subsystems Limitation. Unsupported Nonvirtual Subsystems. Currently Unsupported Simulink Features and Tools. Avoiding Singular Initial Configurations. Virtual Reality Visualization.
4-6 4-6 4-6 4-6 4-6 4-6 4-7

Major Bug Fixes

SimMechanics 2.1.1 includes these important bug fixes made since Version 2.1. If you are viewing these Release Notes in PDF form, please refer to the HTML form of the Release Notes, using either the Help browser or the MathWorks Web site, and use the link provided.
Upgrading from an Earlier Release
SimMechanics 2.1.1 includes the following changes since Version 2.1.

SimMechanics Demos

The some of the demo models of SimMechanics, located in toolbox/physmod/ mech/mechdemos, have been updated.
New Features. Computer-Aided Design Translator for SolidWorks. import_physmod Command for Generating Models. Special Documentation and Examples for CAD Translator. Major Bug Fixes. 2-2 2-2 2-2 2-2 2-4
Upgrading from an Earlier Release. 2-5 Installing the CAD Translator for SolidWorks. 2-5 Known Software and Documentation Limitations SolidWorks Versions Compatibility for CAD Translator Uninstalling and Reinstalling the CAD Translator. Platform Limitations. SimMechanics Users Guide Not Updated. 2-6 2-6 2-6 2-6 2-6

New Features

Two important features have been introduced in SimMechanics 2.1 since Version 2.0.1.
Computer-Aided Design Translator for SolidWorks
SimMechanics 2.1 contains a standalone computer-aided design (CAD) translator for the SolidWorks CAD platform. You install and use the CAD translator with SolidWorks but independently of MATLAB. The translator converts a SolidWorks CAD machine assembly into an XML file in a special Physical Modeling format. The XML file represents the assembly in a form that you can use with SimMechanics to generate a block diagram model dynamically equivalent to the original CAD assembly. To generate models, use the new import_physmod command. See import_physmod Command for Generating Models below. Relative to the main directory of your SolidWorks installation, the CAD translator installs to Applications/TheMathWorks/SimMechanics/.
import_physmod Command for Generating Models
Using the CAD translator (see above), you export a CAD assembly into an XML file in the Physical Modeling format. This XML file represents the assembly in a form that you can use to generate a SimMechanics block diagram model dynamically equivalent to the original CAD assembly. You generate the model from the XML file with the import_physmod command. The generated model contains Bodies and Joints that represent the original parts and degrees of freedom of the assembly. You can edit the generated model to add additional SimMechanics and Simulink blocks.

Special Documentation and Examples for CAD Translator
Special help files accompany the CAD translator. They include two HTML pages and a special PDF book.
These files are independent of the MATLAB help system. Relative to the main directory of your SolidWorks installation, the special CAD translator documentation installs to Applications/TheMathWorks/SimMechanics/help/. The SimMechanics Users Guide has not been updated for this release. Some examples of SolidWorks CAD assemblies are installed to Applications/ TheMathWorks/SimMechanics/examples/. One of the examples is a robot arm assembly. The corresponding XML file in Physical Modeling format, robot.xml, is located in the toolbox/physmod/mech/mechdemos directory, relative to your MATLAB root.
SimMechanics 2.1 includes these important bug fixes made since Version 2.0.1. If you are viewing these Release Notes in PDF form, please refer to the HTML form of the Release Notes, using either the Help browser or the MathWorks Web site, and use the link provided.
This section presents issues involved in upgrading from SimMechanics 2.0 to SimMechanics 2.1.
Installing the CAD Translator for SolidWorks
The computer-aided design (CAD) translator for SolidWorks is independent of MATLAB. You must install it separately, then interface it with SolidWorks.

Installing

You can obtain the self-extracting executable archive SimMechanics_Translator_for_SolidWorks.exe in two ways. When you install SimMechanics 2.1, the archive is installed to toolbox/ physmod/import/standalone/solidworks/win32/, relative to your MATLAB root. Alternatively, you can download the archive independently of MATLAB at www.mathworks.com/products/simmechanics/. Once you have the archive, run it. Click Setup. The setup utility searches your hard drive for the SolidWorks installation location. If you wish to specify the location yourself, you can interrupt the search at any time by pressing the spacebar and then entering a full or partial path to your SolidWorks directory. Once the installer finds SolidWorks, it prompts you to confirm. After you confirm, the installation runs to completion. A README file, with additional instructions, appears after installation.

Linking to SolidWorks

After you have completed the preceding steps to install the translator, start SolidWorks. Choose Add-Ins from the Tools menu. Then select SimMechanics.

Known Software and Documentation Limitations
This section describes known software and documentation problems in Version 2.1.
SolidWorks Versions Compatibility for CAD Translator
The SimMechanics computer-aided design (CAD) translator for Solid Works is designed to work with SolidWorks 2001Plus SP2.0 or later versions.
Uninstalling and Reinstalling the CAD Translator
If you remove the CAD translator for SolidWorks and want to reinstall it, you must reboot your system first after uninstalling. Then you can reinstall the translator.

Platform Limitations

The SolidWorks CAD platform and the SimMechanics CAD translator for SolidWorks are available only for the Microsoft Windows operating system.
SimMechanics Users Guide Not Updated
The SimMechanics Users Guide has not been updated for this release. It is identical to the Version 2.0 Users Guide. The computer-aided design (CAD) translator, installed separately from MATLAB, contains its own special, standalone documentation. This documentation is independent of the MATLAB help system. See Special Documentation and Examples for CAD Translator on page 2-2.
SimMechanics 2.0.1 includes the following bug fix.
Code Compilation with lcc Compiler
Attempting to compile code generated from SimMechanics models using lcc led to errors in Version 2.0. This bug is now fixed.
New Features. Code Generation with Real-Time Workshop. Accelerator Mode and Discrete Events. Expanded Handle Graphics Visualization Features Variable Body Mass and Inertia Tensor. Vectorized Joint Initial Condition Actuation. Damped Linear Force Element Blocks. New Mechanical Branching Bar Block. 2-1 2-1 2-2 2-2 2-2 2-2 2-3 2-3
Major Bug Fixes. 2-4 Upgrading from an Earlier Release. 2-5 Reconfiguring Virtual Reality Visualization. 2-5 Known Software and Documentation Limitations Virtual Reality Toolbox Viewers. Configurable Subsystems Limitation. Unsupported Nonvirtual Subsystems. Currently Unsupported Simulink Features and Tools. Avoiding Singular Initial Configurations. Virtual Reality Visualization. 2-6 2-6 2-6 2-6 2-6 2-6 2-7
A number of important features have been introduced in SimMechanics 2.0 since Version 1.1 (Release 13). This section summarizes them.
Code Generation with Real-Time Workshop
You can now generate stand-alone C code from SimMechanics models using Real-Time Workshop. You can also export it to a PC kernel using xPC Target. SimMechanics 2.0 requires Real-Time Workshop 5.0.1 or later to generate code for the Rapid Simulation (RSIM) target.
Accelerator Mode and Discrete Events
Accelerator mode now supports discrete events associated with the Joint Stiction Actuator and the Point-Curve Constraint.
Expanded Handle Graphics Visualization Features
The Handle Graphics visualization tool of SimMechanics has a large number of new features, including recording and playback of simulation animations.

Variable Body Mass and Inertia Tensor
The Sensors & Actuators library features a new Variable Mass & Inertia Actuator block. This block accepts Simulink signals that vary the mass and/or inertia tensor of the Body block to which it is connected. The block only simulates the changes in the bodys mass properties and does not simulate the effect of the thrust forces or torques resulting from the accretion or ejection of mass.
Vectorized Joint Initial Condition Actuation
The Joint Initial Condition Actuator block now lists all the primitives in the Joint block to which it is connected. You can set initial positions/angles and initial velocities separately for each prismatic and revolute primitive.
Damped Linear Force Element Blocks
This release introduces a new block library, Force Elements, with two blocks representing damped linear harmonic oscillator forces. The Body
Spring-Damper block simulates a damped spring-like force between two bodies. The Joint Spring-Damper block actuates a joint primitive with a damped spring-like force or torque.
New Mechanical Branching Bar Block
The Utilities library features the new Mechanical Branching Bar block. You can use this block to map multiple sensor and actuator lines to one sensor/ actuator port on a Joint, Constraint, or Driver or to one Body coordinate system on a Body.
SimMechanics 2.0 includes these important bug fixes made since Version 1.1 (Release 13).
This section presents issues involved in upgrading from SimMechanics 1.1 to SimMechanics 2.0.
Reconfiguring Virtual Reality Visualization
The SimMechanics documentation assumes that you are using the default Virtual Reality Toolbox viewer. You are free to install and use the blaxxun Contact 4.4 plug-in viewer with a Web browser. The plug-in is supported only on Microsoft Windows platforms.
Installing the Browser Plug-in Viewer
Important tips for installing the blaxxun Contact browser plug-in viewer: Use Version 4.4 of the Contact viewer that ships with the Virtual Reality Toolbox and which is available for Web download from the MathWorks. Do not upgrade to a later version of the blaxxun Contact plug-in. Before installing the blaxxun Contact plug-in, turn off the MATLAB Web server. This allows the blaxxun viewer to register itself over the Web and complete its installation. You can turn the MATLAB Web server back on after installing. You can also avoid installing the Web server altogether. The Microsoft Internet Explorer browser, Version 5.5 or higher, requires manual reconfiguration for use with the blaxxun viewer. You must reset the Java-related network security settings, as explained in the Virtual Reality Toolbox Users Guide.

This section describes known software and documentation problems in Version 2.0.
Virtual Reality Toolbox Viewers
The default viewer functions on UNIX/Linux platforms, including the Mac OS X, but only with virtual reality code written in the VRML 2 standard. The blaxxun Contact plug-in viewer works only on Windows platforms. It can read virtual reality code in both VRML 1 and VRML 2 formats.
Configurable Subsystems Limitation
Simulink configurable subsystems work with SimMechanics blocks only if all of the block choices have consistent port signatures.
Unsupported Nonvirtual Subsystems
For Iterator, Function-Call, Triggered, and While Iterator subsystems cannot contain SimMechanics blocks.
Currently Unsupported Simulink Features and Tools
Certain Simulink features and tools are currently not functioning with SimMechanics blocks or SimMechanics models. Consult the section on SimMechanics limitations.
Avoiding Singular Initial Configurations
Avoid starting a machine in a singular configuration. Its subsequent motion violates assembly tolerances, as the simulation incorrectly removes one or more necessary constraints. A common singular configuration is one where the machine could move in two or three dimensions, but starts in exactly one or two dimensions. Work around an initial singularity by slightly misaligning the joint axes, within assembly tolerances, before starting the simulation. You can tame singularities encountered later in a simulation by selecting the Use robust singularity handling check box on the Constraints pane of the Mechanical Environment Settings dialog.
Virtual Reality Visualization
Consult the Virtual Reality Toolbox Version documentation for product limitations. The presentations of virtual reality visualization in the SimMechanics Users Guide assume that you are using the default viewer.
Rendering and Animation Bugs with the blaxxun Plug-in Viewer
The Virtual Reality Toolbox documentation contains instructions for using the optional blaxxun Contact 4.4 plug-in viewer for Web browsers. You must use the plug-in viewer if you want to view virtual reality scenes written in the VRML 1 format. These are sporadic minor problems with the blaxxun plug-in viewer: If the machine is not initially visible in the browser, try closing the browser and restarting, or try the browsers Refresh or Reload buttons. When you update a diagram, the virtual reality scene might not update as it should. Click the browsers Refresh or Reload buttons to update the scene manually. Microsoft Internet Explorer 6 sometimes fails if you use it as the default browser for virtual reality the blaxxun plug-in viewer. Restart the visualization from the Mechanical Environment Settings dialog if this happens.

1.2 Relative versus absolute coordinate formulations
for DAE analysis, has diculty solving such systems.
Relative versus absolute coordinate formulations
Another consideration is the need for formulating the equations in a way that is both transparent to analysis and that eciently computes the acceleration vector v. In this regard, the two most popular choices of coordinates are relative coordinates and absolute coordinates (or reference coordinates). In the case of absolute coordinates, each body in the mechanism is given six degrees of freedom by describing its position as a reference point on the body and its orientation with a set of four Euler parameters. In this approach, the mass matrix M (q) is particularly simple, being block diagonal. But every interaction between a pair of bodies must be represented by a set of constraint equations. (An example is two bodies connected by a joint.) Invariably this leads to a model that has a large number of conguration variables (the dimension of q can be quite large even for simple systems), and a correspondingly large number of constraint equations. Despite these disadvantages, absolute coordinate formulations are typically the method of choice for commercial packages, where software criteria such as uniformity and ease of constraint representation often dominate. By contrast, relative coordinate approaches minimize the number of coordinates necessary for representing the conguration by implicitly parameterizing certain constraints (for example, joint interactions) between bodies. This re-parameterization is accomplished by restricting the relative motion between bodies to an allowable subspace. This typically results in far fewer variables in the conguration vector q and a corresponding reduction in the number of constraint equations1 , as compared to the absolute coordinate formulation. While the dimension of q and the number of constraint equations is signicantly reduced, a drawback with this approach is that the mass matrix M (q) now becomes dense and the constraint equations more complicated to express. The computational cost of constructing and inverting the mass matrix contributes signicantly to the overall computational cost of the formulation, and so is an important aspect to consider. The evolution of recursive computational techniques makes it possible to factorize the mass matrix and to invert it very eciently, even when using a relative coordinate formulation [2, 7, 10, 11, 21, 22, 28]. For this reason, the SimMechanics solvers use a relative coordinate formulation together with recursive computational procedures. This is in keeping with the general approach to system simulation currently used in Simulink, which is based around ODEs. Clearly it is necessary that the formulation be amenable to integration by all of the existing numerical solvers in Simulink. The choice of relative coordinates allows an important class of mechanisms (for example, many robotic systems) to be simulated without constraint equations. Where constraints do arise, they will be few in number, making it much easier to solve them eciently. There are systems where the number of constraints is proportional to the number of bodies, even when relative coordinates are used. Such systems are pathological and do not arise often in practice.

2.1 Analysis of simple chain systems
Analysis of simple chain systems
Consider the n-link serial chain in Figure 1, where the links are numbered consecutively from the tip of the chain to the base, which also acts as the inertial reference frame. Each body is connected to two joints, an inboard joint (closer to the base), and an outboard joint (closer to + the tip). An arbitrary joint, the kth joint, has an inboard side labeled Ok attached to body k +1 and an outboard side labeled Ok attached to the kth body. On the kth joint, the vector from the outboard side of the joint Ok to the center of mass on the kth body is denoted by p(k), and the vector from Ok to Ok1 is denoted by r(k, k 1). This vector plays an important role in shifting forces and velocities between bodies. The spatial velocity V (k) R6 at the outboard side of the kth joint is dened by V (k) = [(k)T , v(k)T ]T , where (k) R3 is the angular velocity of the kth body and v(k) R3 is the linear velocity of the kth body at the point Ok. Here we assume that both vectors are expressed in inertial coordinates. The spatial force f (k) R6 applied by the kth joint to the kth body at the point O(k) is dened by f (k) = [T (k)T , F (k)T ]T , where T (k) R3 is the applied torque and F (k) R3 is the applied force, again expressed in inertial coordinates. The spatial inertia matrix M (k) is dened as M (k)= J(k) m(k)(k) p m(k)(k) m(k)I3 p R66 ,
where m(k) is the mass of the kth body, and J(k) R33 is the inertia tensor of the kth body about the center of mass, computed in inertial coordinates. Given a vector p R3 we let p R33 denote the cross product matrix generated from the vector p. Thus p(k) denotes the cross product matrix generated from p(k). Denote the time derivative of V (k) by (k). The following kinematics equations describe the motion of the chain: V (k) = T (k + 1, k)V (k + 1) + H T (k)q(k), V (n + 1) = 0,

k = n, n 1,. , 1

(k) = (k + 1, k)(k + 1) + H (k)(k) + a(k), (n + 1) = 0, q

k = n, n 1,. , 1 (5).

Here q(k) Rnqk is the vector of conguration variables for the kth joint, the columns of H T (k) R6nqk span the relative velocity space between the kth body and the (k + 1)th body, and T (k + 1, k) is the transpose of the matrix (k + 1, k) = I3 r(k + 1, k) 0 I3 ,
where the vector r(k + 1, k) is the vector from the outboard size of the (k + 1)th joint to the outboard side of the kth joint. The vector a(k) R6 represents the Coriolis acceleration of the kth body and is given by a(k) = (k + 1)(k) (k + 1)(v(k) v(k + 1)) .
The role of various terms can be claried with an example. Suppose the kth joint is a simple

kth hinge + Ok p(k+1) r(k+1,k)
+ Ok-1 Ok+1 (k+1)th hinge + Ok+1 r(k,k-1) Ok-1 k-1th hinge

Toward base

Toward tip
Figure 1: Bodies and joints in a serial chain.
revolute joint. This joint only allows rotation between the (k + 1)th body and the kth body. In this case, the variable q(k) represents the rotation angle between the two bodies and the matrix H(k) is given by H(k) = [uT , 0T ], with uk R3 being the unit vector along the instantaneous 31 k axis of the revolute joint (in inertial coordinates). The kinematic relationships for more complex joints can also be expressed in this form. The matrix (k + 1, k) shifts spatial forces (according to the usual rules for shifting forces between points on a rigid body) from the point O(k) to the point O(k + 1) while the transpose (or adjoint) shifts spatial velocities from O(k + 1) to O(k). Given the joint velocity degrees of freedom (DoFs) q(k) and acceleration DoFs q (k), these equations provide a simple recursive procedure for determining the velocities and accelerations of the bodies constituting the chain. The equations of motion, formulated about the reference points O(k), are f (k) = (k, k 1)f (k 1) + M (k)(k) + b(k) , T (k) = H(k)f (k). The vector b(k) R6 is the gyroscopic force vector at O(k): b(k) = (k)I(k)(k) m(k) (k) (k)p(k) . f (0) = 061 , k = 1,. , n (6)
The second equation in 6 is obtained by performing a virtual work balance across the kth joint (a massless entity). With the externally applied generalized torque vector2 for the joint given
These could be torques or forces applied by motors around various joint axes.
by T (k), we obtain T T (k)qk = f T (k)H T (k)qk , qk Rnqk.

2.2 Recursive solution

Since qk is arbitrary, this gives the desired relationship between the externally applied generalized torque T (k) and the reaction force f (k). Readers familiar with optimal control theory will recognize equation 5 and equation 6 as being similar to those that arise in the optimal control equations of discrete dynamic systems [6]. After substituting f (k) = P (k)(k) + z(k), sweep methods lead to a discrete Riccati equation for P (k) (which can be solved recursively) and then directly to the desired joint accelerations. The solution, when expressed in matrix form, leads to a square factorization of the mass matrix and its inverse.

Recursive solution

We can express these recursive relationships in matrix form as follows. Sum the velocity recursion and use the fact that (i, i) = I6 and (i, k)(k, j) = (i, j). The second relationship follows from the fact that shifting a spatial force from Oj to Ok and then from Ok to Oi is equivalent to shifting the force directly from Oj to Oi. Then

V (k) =

T (i, k)H T (i)qi
A natural denition of the matrix operators follows: H T = diag[H T (1),. , H T (n)]. The operator is dened by: = I66 (2, 1). 0 I66. 0 0. .

This procedure yields a linear system of equations for q (k) and (k). A back substitution gives the ABM, with the joint accelerations q (k) computed using the circled pivot blocks. Since all op erations are done block-wise, the computation count is clearly O(n). These recursive techniques exploit the structure and sparsity pattern of the matrices that arise in the multibody dynamics problem to eciently simulate systems represented by acyclic graphs. Recursive techniques are easily extended to arbitrary tree topologies and provide the underpinning of an ecient computational approach to the general case where constraints arise.

Constrained systems

A general purpose solver must be capable of dealing with constrained mechanical systems where the underlying graph is cyclic. Systems that have a cyclic graph are reduced to an open spanning tree by the removal of a subset of the joints. The joints in this cut-set are replaced by a set of constraint equations to ensure that the spanning tree undergoes precisely the same motion as the original mechanism with the cyclic topology. The choice of joints in the cut-set is often not unique and has a direct eect on the number of generated constraint equations and the sparsity pattern of the resulting constraint Jacobian. The problem of making a suitable choice is related to the problems of ordering columns in sparse matrices and can be addressed with the techniques of graph theory [9]. One such choice is to choose the cut-set so as to maximize the number of independent loops, thereby transforming the constraint equations into smaller sets of independent equations which are easier to solve. SimMechanics automatically selects the
3 CONSTRAINED SYSTEMS 3.1 Ecient computation of the Lagrange multipliers
cut-set but does allow users to determine this set manually. Astute analysts should be able to intuitively select a cut-set that exploits the properties of the specic mechanism being analyzed, if necessary.
Ecient computation of the Lagrange multipliers
The equations of motion of the spanning tree can be expressed using the recursive relationships outlined in section 2. This step results in the descriptor form given by equation 1, equation 2, and equation 3. Since the descriptor form is an index-3 DAE, it is necessary to reduce the index by explicitly dierentiating the constraints. If we dierentiate the constraint g(q, t) = 0 in equation 3 twice, we get g t 2 (GHv) Hv + 2 G Hv + g GH v + q t t2 GHv + = 0 = 0. (9) (10)
Substituting the equation of motion in equation 2 into equation 10 gives M v = f + H T GT GHM 1 H T GT = f The vector f is given by f = GHM 1 f G 2g (GHv) Hv 2 Hv 2 q t t. (13). (11) (12)

Equation 11 and equation 12 can be treated as an index-1 DAE [29], or we can solve for directly, and treat the system as an ODE on a manifold. In the index-1 DAE approach, it is possible to exploit the structure of the equations of motion, together with carefully designed multi-step integration schemes, to implement very ecient domain specic solvers for multibody systems [29]. Such an approach still suers from numerical drift, however, and is not an option for SimMechanics, which must work seamlessly with the current suite of ODE solvers in Simulink. For this reason, we are forced to consider the ODE approach. There are two problems with this approach: First, we need to compute the Lagrange multipliers eciently. Second, because we have reduced the DAE to an ODE on a manifold, blindly integrating the resulting ODE numerically gives a solution that drifts away from the desired manifold given by g(q, t) = 0 [5, 27]. The same is true of the index-1 DAE approach [29]. An alternative procedure is to eliminate the Lagrange multipliers from the analysis and to make use of coordinate partitioning to split the variables into independent and dependent sets [19, 32]. The existence of such a partitioning, at least locally, follows directly from the implicit function theorem and is equivalent to an explicit local parameterization of the manifold. Our principal objection to this approach is that, if the subset of independent variables is integrated, it results in a projection of the mass matrix that eliminates the structure that we exploit in the O(n) computations. In practical applications, it is necessary to carefully monitor and recompute the parameterization as the simulation proceeds
3.2 Handling singularities
since the choice of independent and dependent variables is somewhat arbitrary and is unlikely to be globally applicable. The approach taken in SimMechanics is a direct approach in some ways similar to the direct approach used to solve the index-1 DAEs arising in models that have algebraic loops. The Lagrange multipliers are computed using ecient recursive computations to evaluate the coecient matrix = GHM 1 H T GT (14) and the vector f. To demonstrate this, consider the problem of forming the coecient matrix. The velocity Jacobian can be factorized to give G = AT H T , (15) where the conguration-dependent matrix A(q, t) is determined by the nature of the constraints. Using this factorization, it is possible to generate the coecient matrix in O(nc n) operations, where nc is the number of independent velocity constraint equations. Similarly the vector f can be computed very eciently using the factorized form of the velocity Jacobian in equation 15. This allows f to be computed as the sum of a free acceleration solution (where the constraint forces are set to zero) and the product of the matrix A and its time derivative with the spatial acceleration vector determined from knowledge of q and v (but not v). Again this computation can be implemented using the recursive techniques outlined in section 2. Solving for the Lagrange multipliers, still requires the coecient matrix to be factorized. This is an O(n3 ) operation. Here we see the reasoning behind the use of a relative coordinate c formulation. Since the number of constraint equations nc is typically small, it is still possible to eciently compute the Lagrange multipliers. We note here that since the matrix is symmetric and positive-denite, it is possible to solve for the Lagrange multipliers using an iterative solution technique like conjugate gradient iteration, all implemented using the matrix factorizations already discussed. In the absence of roundo eects, this approach could generate a solution in O(nc n) operations. However, in the general case where the coecient matrix might be ill-conditioned, and in the absence of a suitable preconditioner, we choose to use the direct approach.

Handling singularities

An important issue arises when implementing the solvers to determine the Lagrange multipliers. Mechanical systems can pass through singular congurations during their motion. Singular congurations exist when the number of independent constraint equations is reduced on a nite set of points in the conguration space. At such points, the coecient matrix becomes rank decient, although it is still possible in principle to determine a vector of Lagrange multipliers (no longer unique) that satises the determining equation. Unfortunately, in the vicinity of these singular points the ill-conditioning of makes the computation of extremely sensitive to roundo errors in f [31, 17]. The most popular approach to deal with singular congurations is the use of damped least squares. This approach minimizes the sum, f
3.3 Compensating for numerical drift
where is the damping factor. Here controls the tradeo between keeping the residual small and keeping the norm of small in the vicinity of a singularity. Unfortunately, there is no simple or generic method for determining a suitable damping factor. So we use a slightly dierent approach based on a truncated QR decomposition of. This is similar, in some respects, to the truncated singular value decomposition approach [17] but is computationally more ecient. The approach is based upon QR decomposition of with full column pivoting [15] E = QR , (17)
where E is a permutation matrix determined by the pivot selection. The pivots are chosen so that the absolute values of the terms on the diagonal of R appear in numerically decreasing order: |R(1, 1)| > |R(2, 2)| >. |R(nc , nc )|. It is well known that this algorithm can be used to estimate the numerical rank of a matrix [15], and in our computations, we estimate the rank r according to r = max : R(k, k) > nc M R(1, 1) , (18)
where M is the machine epsilon (distance from one to the next oating point number)3. If r < nc the decomposition can be written as E = Q Q R11 R0 , (19)
where Q Rnc r , Q Rnc (nc r) , R11 Rrr , and R12 Rr(nc r). The solution is then =E

1 R11 QT f 0(nc r)1

which is the minimum norm solution (in the Euclidean metric) in the projected subspace. SimMechanics oers users the choice between a solver based on Cholesky decomposition of the coecient matrix, which is numerically more ecient, and a solution based on the truncated QR decomposition, which is more robust in the presence of singularities.
Compensating for numerical drift
In the absence of roundo, index reduction and integration of the resulting ODE would be perfectly adequate. Unfortunately, numerical drift of the solution from the invariant manifold g(q, t) requires alternative schemes. One such scheme, based on coordinate partitioning, has already been mentioned. Here the coordinates are split into independent and dependent sets. Finding a robust detection method for changing the parameterization is the most challenging aspect of the scheme. This scheme is not suitable for real-time HIL simulations, an important consideration for SimMechanics as it supports code generation for real-time deployment. Other schemes are based on stabilization and coordinate projection [3, 4, 5]. Stabilization involves the addition of extra terms to the equation of motion. These terms vanish on the manifold

3 The only truly reliable manner of determining the numerical rank of a matrix is the singular value decomposition [12], although the QR decomposition is known to work very well in practice.
3.4 Coordinate projection
g(q, t) = 0 but have the eect of making the solution asymptotically attractive to the manifold. If the solution does drift o the manifold, it is ultimately attracted back onto it, although there are no predened bounds on the extent of the drift. The most popular stabilization technique is Baumgarte stabilization. Its simplicity has made it a popular choice in engineering applications. But Baumgarte stabilization requires parameterization, and there is no known generic procedure for choosing these parameters to make the stabilization robust (see [14] for attempts in this regard). Indeed the choice of suitable parameters depends on the discretization scheme used to integrate the equations of motion [5], a serious practical limitation.

Coordinate projection

Coordinate projection involves the numerical discretization of the ODE. At the end of the discretization step, the solution is projected onto the invariant manifold. In [27], various practical aspects of these projection schemes are studied in some detail. For example, the eect of coordinate projection on step size selection, event location, and order of convergence are all considered. It is shown how the codes in the present Simulink ODE suite, ode23, ode45, ode113 and ode15s, can all be adapted to allow for coordinate projection without compromising accuracy or eciency [25]. A fundamental change to Simulink to enable mechanical simulation with SimMechanics has been the addition of a projection method to the ODE suite that can be called once the discretized solution has been updated following the acceptance of a successful step based on the error estimates.The scheme works for the one-step Runge-Kutta formulas ode23 and ode45, as well as for variable-step variable-order codes like ode113 [27]. Finally, the standard theory for convergence in BDF codes like ode15s is still applicable when projection is carried out in this way. An important property of the Simulink solver suite is its ability to localize and detect events using discontinuity locking (to ensure that the integrators see a continuous vector eld) and its ability to output solutions at any time in the interval of integration using continuous extension formulae [25]. Both of these features are aected by projection. Simulink allows users to rene the output from the solvers using highly ecient interpolation schemes, and these interpolated outputs typically satisfy the invariants to an accuracy comparable with the accuracy of the numerical solution. In SimMechanics, the interpolated values are further projected to ensure that the invariants are satised. It could be argued that projection of the interpolated values reduces eciency. Our approach avoids surprising users who might notice that the interpolants do not satisfy the invariants. More critical is the use of projection in event location. The event location capabilities of the current Simulink solvers are also built around the continuous extension formula and the use of switching functions to determine the exact location of events in time. To ensure that events are located correctly, the outputs of the interpolants are projected before sampling the switching functions. This makes event location more expensive, but is necessary for robust event detection. For example, the projection approach is appropriate for a one-step method used to compute an approximate solution at time tn+1 from a solution at tn. The step size is h, and tn+1 = tn + h.

The solution takes the form
3.5 Stabilization for code generation

qn+1 vn+1

= h(tn , qn , vn ).
If we consider the nonlinear position constraint g(q, t) = 0, the predicted position variables qn+1 are projected onto the closest point on the manifold4 , denoted by qn+1. Linearizing the constraints about the predicted value qn+1 gives
0 = g(tn+1 , qn+1 ) = g(tn+1 , qn+1 ) + G(tn+1 , qn+1 )(qn+1 qn+1 ) + O( qn+1 qn+). (22)
To leading order, the projected solution qn+1 is obtained from solving G = g qn+1 = If the Euclidean norm is used, = GT (GGT )1 g. (25) In [29], the projection process is repeated, and so is referred to as sequential projection, until a suitable level of convergence is attained. This is much more ecient than forming the full set of Karush-Kuhn-Tucker (KKT) equations to determine the exact minimization on the manifold and is easily motivated by the fact that error control has made the predicted solution qn+1 close to the manifold already. In practical applications, it is advantageous to perform the projection in a weighted norm. The projection weights reect the weights used in the error control of the ODE solver. Suitable convergence is somewhat ambiguous. In SimMechanics, we oer users the ability to terminate the projection process based on relative and absolute tolerances for. We also allow users to carry on the projection process to completion. This implies that the iteration is continued until the new iterates fail to reduce the residual due to rounding error. This is similar to the current approach taken in solving algebraic loops using direct methods in Simulink. The latter approach could be slightly more expensive.

(23) . (24)

Stabilization for code generation
Numerical projection aects an important feature of SimMechanics, the generation of standalone code for faster simulation and deployment on rapid prototyping systems with Real-Time Workshop R. For systems that have cyclic graphs, code implementation requires projection to be carried out in real time, and neither of the projection schemes outlined thus far are deterministic, in the sense that the number of iterations needed cannot be predicted. Stabilization codes thus provide an appealing alternative to the projection approach for real-time applications. The generalizations of the Baumgarte technique in [5] are appealing since they can be implemented without the need for selecting problem-specic parameters. For notational convenience, suppressing explicit time dependence, we can write the reduced ODE representation

4 How close depends on the metric and potentially on the weighting. Usually, if weights are used in ODE error estimation, the same weights should be used in the projection scheme.
(with the Lagrange multipliers eliminated) of equation 1, equation 2 and equation 3 as: z = f (z) 0 = c(z) = g(q) G(q)H(q)v. (26) (27)
Here the variable z is dened as z T = [q T , v T ]T. Reference [5] considers the family of stabilization methods given by z = f (z) F (z)c(z) , (28) where > 0 is a parameter, and F = D(CD)1 and C(z) = The matrix C(z) is given by C(z) = c z. (29)

(GHv) q

and D(z) is chosen such that CD is nonsingular; for instance D= which gives, F = H T GT (GH H T GT ) 0 H T GT (GH H T GT )1 I
(GHv) HGT (GH H T GT )1 q

H T GT 0

0 H T GT
(32) By discretizing the stabilization term independently of the dierential equation, one can show that an optimal choice of is proportional to the inverse of the step size [5]. The stabilization scheme is implemented as follows. Start with (qn , vn ) at time tn and integrate with some ODE numerical method to advance the solution to (qn+1 , vn+1 ) at t = tn+1. Stabilize, by modifying the solution qn+1 qn+1 F (qn+1 , vn+1 )c(qn+1 , vn+1 ). (33) = vn+1 vn+1 Since the term (GHv)/q is expensive to compute, a cheaper alternative is to use F = M 1 H T GT (GHM 1 H T GT ) 0 M 1 H T GT (GHM 1 H T GT )1 , (34)
obtained by neglecting the (GHv)/q term and setting D= M 1 H T GT 0 M 1 H T GT 0. (35)
3.6 Redundant constraints
This is similar to projection. However, this projection is carried out in a metric dened by the positive-denite mass matrix M. This scheme is particularly ecient because the term = GHM 1 H T GT has already been computed and factorized when the Lagrange multipliers were computed. The stabilization step is computationally almost free and ideal for real-time application. Finally, an even better choice is to use the cheaper variant of F , but to apply the scheme twice [3]. This adds negligible cost but signicantly increases the accuracy. This is the scheme that users can choose when simulating systems in SimMechanics with stabilization and is the scheme used to generate code for constrained systems.

Redundant constraints

When using a mechanical simulation tool it is easy to generate redundant constraints. In this case the constraint Jacobian G fails to be surjective, and the redundant constraints need to be identied and removed. An additional consequence is that the constraint forces are no longer unique. This may not be important if the primary interest is the motion of the system, but it could be very important if the constraint forces themselves are of interest. Also SimMechanics allows users to move certain joints and constraints with known time trajectories (motion drivers). This is an important feature, often used in computed-torque control schemes and biomechanics applications. Here it is possible for the user to overspecify the motions or to impose inconsistent motions between the driven joints. SimMechanics uses a relative coordinate formulation, and the motion drivers are treated differently according to the nature of the graph associated with the system. For systems that have an acyclic graph, the motions induced by motion drivers can never be inconsistent and they can be dealt with very eciently using the recursive methods of section 2 without the need for any constraint equations. It is only systems that have a cyclic topology where inconsistencies can arise, and for these systems we need to be careful to distinguish between motion drivers applied to the joints in the spanning tree and motion drivers applied to joints in the cut-set. Drivers that are applied to joints in the spanning tree can still be dealt with without adding additional constraint equations using recursive methods. They must now be checked for consistency however. Drivers applied to joints in the cut-set do result in additional constraint equations, and if these are redundant, they can result in inconsistent motions. In this section we discuss how the redundant constraints are resolved and how potential inconsistencies are detected in SimMechanics. These issues can also arise during simulation when joints lock or unlock due to stiction behavior. We delay discussion of this topic until section 6. The primary tool used in this analysis is the QR decomposition with full column pivoting [15]. We begin by reordering the rows and columns of G G(q, t) = G11 (q) G12 (q) G21 (q) G22 (q) q1 q2 , (36)
T T where q T = [q1 , q2 ] has been partitioned into a set of unknown conguration variables q1 : R n1 and a set of known, time-dependent variables, q : R Rn2 , n + n = n. The columns of R 2 q G have been partitioned accordingly. Similarly the rows of G have been ordered such that the
rows associated with the homogeneous constraint equations appear rst and the rows associated with motion drivers applied to joints in the cut-set appear below them. The velocity constraints5 then take the form 0 G11 (q) G12 (q) q1 =. (37) g21 G21 (q) G22 (q) q2 t Redundancy and consistency are detected before the simulation is started. The rst task is to detect potential inconsistencies in the motion drivers applied to the joints in the spanning tree. To this end, we determine the eective rank of the following submatrix QR = G11 G21 E , (38)

where E is a permutation matrix determined by the pivoting strategy in the QR decomposition. Using the technique discussed in section 3.1, we determine the numerical rank of this matrix. To determine potential inconsistencies in the motion drivers, and to provide feedback to identify which drivers cause the diculty, the columns of the submatrix [GT GT ]T are appended, one by one, and the QR decomposition recomputed. After k such stages, we obtain the following decomposition G11 G12 I 0 E. (39) QR = G21 GFk Here Fk is the matrix that selects the rst k columns of [GT GT ]T. If, in this process, the rank of the augmented matrix increases, we have identied (1) that a potential inconsistency exists, and (2) the driver causing the inconsistency. Increasing the rank of the augmented matrix, and the fact that the elements of q2 are arbitrary quantities determined by the user, implies that a solution to the velocity constraint equation may not exist. Of course it is possible for the user to dene the motions of the drivers in such a way that they are consistent. It is impossible to determine whether the user will do this, and consequently an error is reported if inconsistencies are detected. Having veried that the motion drivers in the spanning tree are consistent, we move to the task of eliminating redundancies among the homogeneous constraints. To achieve this we compute the QR decomposition of the submatrix QR = GT E 11. (40)
This allows us to identify and eliminate the dependent rows. Finally we append the constraints generated by motion drivers applied to joints in the cut-set, one by one. After k stages, we obtain the following decomposition QR =

T GT GT Fk 11 21

where Fk selects out the rst k rows of G21. If the rank of the augmented matrix increases, we have again identied a possible inconsistency between the motion drivers and the source of the
For simplicity, we ignore nonholonomic constraints here, but they are easily included in the analysis.

inconsistency.

3.7 Consistent initial conditions
Typically mechanical systems have few explicit motion drivers, and this computation need only be performed once (except when discrete topology changes occur due to stiction). One diculty with this approach arises when a simulation is started from a singular conguration. In this case, the analysis eliminates a constraint or constraints that are instantaneously redundant. Consequently the system likely violates the constraint as the simulation progresses. SimMechanics detects the constraint violation, but at present does not automatically detect whether the simulation is being started from such a singular conguration.

Linearization

A large class of dynamical systems can be satisfactorily controlled using the highly developed tools of linear control theory. To facilitate this it is necessary to obtain linear models of the system dynamics that accurately reect the behavior of the system in the neighborhood of a specied nominal trajectory. For mechanical systems, we start with the reduced ODE form q = H(q)v v = F (t, u, q, v) , where F : R Rnu Rnq Rnv Rnv is given by: G 2g (GHv) Hv 2 Hv 2 )]. (46) F = M 1 [f + H T GT (GHM 1 H T GT )1 (GHM 1 f q t t We have included the dependence of the dynamics on exogenous inputs (external forces and motor torques) through the vector u : R Rnu. Linearization about a nominal trajectory (, q , v ) proceeds by introducing perturbations (u, q, v) and expanding the equations about u the nominal trajectory to get: F (t, u, q , q )u+ F (t, u, q , v )q+ F (t, u, q , v )v+O( u 2 , q 2 , v 2 ). v + v = F (t, u, q , v )+ u q v To rst order, this gives a time-varying representation of the perturbative dynamics: v = F F F (t, u, q , v )u + (t, u, q , v )q + (t, u, q , v )v u q v. (44) (45)
Simulink can numerically linearize dynamic systems to generate a linear representation. Recently, Simulink was extended to incorporate the analytic derivatives, where these exist, either determined internally or provided by the user, of blocks or subsystems by using linear fractional transformations. This is important since, without this extension it would not be possible to generate linear models that accurately reect the dynamics of cyclic systems. Very simply,
any perturbation (numerical or otherwise) introduced to determine the linear dynamics, must also satisfy the constraints to the same order of accuracy as the linearized dynamics. Since Simulink is essentially a general purpose solver, the perturbation is introduced to linearize the reduced dierential equations of motion without knowledge of the constraint manifold. Since cyclic systems have constraints, only a subset of the perturbed motions satisfy equation 44 and equation 45 to rst order and also satisfy the constraints to the same order. Acyclic systems can be linearized in the standard manner, so we concentrate on the cyclic case. Here a linearization of the mechanical subsytem is provided to Simulink by SimMechanics. Simulink can then linearize a system containing other blocks and subsystems that are connected to the mechanical subsystem. Since the linearization of an interconnection can be obtained by the interconnection of the linearizations, Simulink can construct the linearized model for the whole system (mechanical and other) very eciently. The task then falls on SimMechanics to register a linearization method with Simulink and to construct the desired linearization of the mechanical components. For cyclic systems this means linearizing the ODE on the constraint manifold. The perturbation must therefore satisfy g(, t) + G(, t)q = 0 q q g = 0. G(, t) q + G( + q, t)q + q q t (47) (48)

REFERENCES

ode15s.

Conclusions

In this paper we have addressed some of the interesting issues that arise when introducing mechanical simulation capability into an existing ODE-based simulation environment. We addressed the issue of computational eciency by demonstrating how a relative coordinate formulation can be implemented using recursive computations to produce ecient O(n) algorithms. The problems of Lagrange multiplier determination and preventing numerical drift were also addressed. A number of strategies were outlined. Each has dierent properties which make it applicable in dierent situations. Issues associated with redundant constraints and consistent motion drivers were discussed briey. A strategy based on QR decomposition was put forward as an ecient and reliable way to eliminate redundant constraints and detect consistency problems. The important topic of linearization on a manifold and how this is currently achieved in Simulink was also examined. Finally, we mentioned some of the issues that arise when simulating mechanical systems that undergo discrete topology changes and the implications for real-time simulation.

References

[1] E. Allgower, K. Georg, Numerical Continuation Methods: An Introduction, Springer-Verlag, 1990 [2] W. W. Armstrong, Recursive Solution to the Equations of Motion of an n-link Manipulator, Proceedings of the 5th World Congress on Theory of Machines and Mechanisms, Vol. 2, 1979, pp. 1343-1346 [3] U. Ascher, H. Chin, L. Petzold, S. Reich, Stabilization of Constrained Mechanical Systems with DAEs and Invariant Manifolds, J. Mech. Struct. Machines, 23, 1993, pp. 135-158 [4] U. Ascher, L. Petzold, Stability of Computational Methods for Constrained Dynamics Systems, SIAM J. SISI, 14, 1993, pp. 95-120 [5] U. Ascher, H. Chin, S. Reich, Stabilization of DAEs and invariant manifolds, Numer. Math., 67, 1994, pp. 131-149 [6] A. E. Bryson, Y. Ho, Applied Optimal Control, Blaisdell Publishing, 1969 [7] B. Cloutier, D. Pai, U. M. Ascher, The Formulation Stiness of Forward Dynamics Algorithms and Implications for Robot Simulation, Proceedings of the IEEE Conference on Robotics and Automation, 1995 [8] J. Dennis, R. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, SIAM, 1996
[9] I. Du, A. Erisman, J. Raid, Direct Methods for Sparse Matrices, Clarendon Press, Oxford, 1986 [10] R. Featherstone, The Calculation of Robot Dynamics Using Articulated-Body Inertias, International Journal of Robotics Research, 2(1), 1983, pp. 13-30 [11] R. Featherstone, Robot Dynamics Algorithms, Kluwer Academic Publishers, 1987 [12] G. Golub, C. Van Loan, Matrix Computations, Johns Hopkins University Press, 1996 [13] E. J. Haug, Computer Aided Kinematics and Dynamics of Mechanical Systems, Volume 1, Basic Methods, Allyn and Bacon, 1989 [14] E. Haug, R. Deyo, Real-Time Integration Methods for Mechanical System Simulation, NATO ASI Series, Springer, 1991 [15] N. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, 1996 [16] T. R. Kane, D. A. Levinson, Multibody Dynamics, Journal of Applied Mechanics, 50, 1983, pp. 1071-1078 [17] A. A. Maciejewski, C. A. Klein, Numerical Filtering for the Operation of Robotic Manipulators through Kinematically Singular Congurations, Journal of Robotic Systems, 5(6), 1988, pp. 527-552 [18] D. K. Pai, U. M. Ascher, Paul G. Kry, Forward Dynamic Algorithms for Multibody Chains and Contact, Proceedings of the IEEE Conference on Robotics and Automation, 2000 [19] T. Park, E. Haug, A Hybrid Constraint Stabilization Generalized Coordinate Partitioning Method for Machine Dynamic Simulation, J. Mech. Trans. Auto. Des., 108 (2), 1986, pp. 211-216 [20] R. E. Roberson, R. Schwertassek, Dynamics of Multibody Systems, Springer-Verlag, 1988 [21] G. Rodriguez, Kalman Filtering, Smoothing and Recursive Robot Arm Forward and Inverse Dynamics, IEEE Journal of Robotics and Automation, 3(6), 1987, pp. 624-639 [22] G. Rodriguez, A. Jain, K. Kreutz-Delgado, A Spatial Operator Algebra for Manipulator Modeling and Control, The International Journal of Robotics Research, 10(4), 1991, pp. 371-381 [23] W. O. Schiehlen, Multibody System Handbook, Springer-Verlag, 1990 [24] A. A. Shabana, Dynamics of Multibody Systems, Wiley, 1989 [25] Lawrence F. Shampine, Mark W. Reichelt, The MATLAB ODE Suite, SIAM J. Sci. Comp., 18, 1997, pp. 1-22 [26] Lawrence F. Shampine, Mark W. Reichelt, Jacek A. Kierzenka, Solving index-1 DAEs in MATLAB and Simulink, SIAM Review, 41, 1999, pp. 538-552