TSi Dynamics Version 1.1

Introduction

TSi Dynamics is a Mathematica package containing functions that support the assembly of mathematical models and simulation models for mechanical systems.  In its current form, TSi Dynamics accommodates tree structures composed of combinations of rigid and flexible bodies and simple and compound kinematic joints, and which may also have supplemental algebraic and/or differential constraints. Thus, systems with closed loops and nonholonomic joints can be treated. 

TSi Dynamics creates fully nonlinear, explicit models. The model building process requires that the user provide defining data for individual joints and bodies, as well as the system interconnection structure. The joints are defined in terms of their primitive action parameters from which all the required kinematic relations are derived. Thus, a user can contrive unusual joint configurations and is not restricted to a predefined set of standard joints. The equations are formulated in Poincar‚'s form of Lagrange's equations that admits the standard Lagrange equations as a special case. However, Poincar‚'s form allows the exploitation of quasi-velocities which can significantly simplify the equations of motion.

TSi Dynamics assembles the equations of motion and it contains functions that facilitate simulation. Simulation may done within Mathematica  or externally using some other simulation software. TSi Dynamics includes functions which assist in performing simulations using the NDSolve function in Mathematica and other functions to construct a C-code subroutine that compiles as a "MEX-file" for use with MATLAB/SIMULINK.

Getting Started

This notebook illustrates the tools provided in TSi Dynamics including a few new functions contained in Version 1.1. TSi Dynamics is contained in the directory TSiDyn_1. That directory also contains several examples. TSiDyn_1 can be placed anywhwere you like, but a search path must be defined for it so that it will be loaded properly.  For example, if you place it in the mathmatica packages directory with full path name "c:\math\packages\tsidyn_1", then the following command does this.

$Path=Join[$Path,{"c:/math/packages/tsidyn_1/"}]

If you don't want to do this every time you want to use the package, then the above command can be placed in the init.m notebook found in the Packages directory. That way the path will be set up every time Mathematica is initialized.

Once the search path is defined, TSi Dynamics can be loaded with the command

<<TSiDynam`

The Basics of TSi Dynamics

What's in the Package?

To find out what functions are available in TSi Dynamics simply use the command:

?TSiDynamics

Information about any of these functions can be obtained in the usual way with ?Function. Most users will find that only a few of these functions are necessary for creating models and simulations. The most important are:

   Joints: generates all kinematic functions for list of joints
   
   CreateModelSim: creates the dynamic model parameters for tree structures from system data
   
   EndEffector & EndEffectorVelocity: creates the position, orientation and velocity of a body frame located at specified nodes as functions of configuration coordinates
   
   SpringPotential & LeafPotential: the former creates the potential energy function of a spring (defined by a local energy function) connected between specified nodes, the latter performs a similar function for energy storage elements that act between a system node and the space frame
   
   GeneralizedForce & DamperForce: the former is a general function for forming the generalized force vector from a local force acting at a specified node, and the latter is specific for a damper (defined by a local disspation function) connected between specified system nodes
   
   AlgConstrainedSys & DiffConstrainedSys: these account for algebraic and differential constraints added to a tree structure
   
   MakeODEs & CreateMEXFile: the former assembles the ordinary differential equations from the model parmeters generated by other functions into a form for numerical computation within Mathematica, the latter builds a C source function for computation in MATLAB/SIMULINK.
   
After describing the basic data structures used by TSi Dynamics we will illustrate the use of these and some other functions.

Data Structures

Tree Structures

TSi Dynamics builds models for mechanical systems which have a tree topology. Chain structures are a special case. Every system contains a base reference frame which is designated body '0'. Otherwise, bodies and joints can be numbered arbitrarily.

A tree is composed of a set of defining chains. For instance consider a tree composed of the following sequences of bodies:

                                       0,1,2,4      0,1,2,3,5      0,1,2,3,6

All definingchains of any tree will start with body 0, so we need not list it. However, the body sequences alone do not adequately define a tree. For instance bodies 5 and 6 both connect to body 3, but they will do so through different joints. This information can be provided by defining each chain as an ordered list of pairs Ä each pair consisting of a body and its inboard joint: {inboard joint, body}. For example:

            {{1,1},{2,2},{5,4}}      {{1,1},{2,2},{3,3},{4,5}}      {{1,1},{2,2},{3,3},{6,6}}
            
  indicates that body 5 connects to body 3 at joint 4, and body 6 connects to body 3 at joint 6.
  
A tree is defined by the data structure:

    Tree = {list of chains}
    chain = ordered list of pairs {inboard joint, body}
              =  {{first inboard joint, first body},...,{last inboard joint, last body}}

Joints

Joints constrain the relative motion between two bodies, or more precisely, between two reference frames Ä one fixed in each body on either side of the joint. In TSi Dynamics, joints are of two types: simple and compound. Simple joints directly define relative motion between the two reference frames, whereas compound joints require intermediate reference frames to define the overall motion.

Simple joints are characterized by the number of degrees of freedom, r, an r-vector of joint quasi-velocities, p, and a 6xr joint map matrix, H. Across the joint the relative velocity vector (Dw, DV) is Hp. Moreover, H is a constant (independent of the joint configuration) and the columns represent the joint action axes in the outboard frame (by convention). A compound joint is equivalent to a sequence of simple joints. Thus, it is neccessary to define a set of numbers that represent the degrees of freedom associated with each intermediate frame and a corresponding set of (constant) joint map matrices. When defining a joint in TSi Dynamics, it is necessary to also assign names for both the joint quasi-velocities and joint configuration variables.

A k-frame compound joint with n degrees of freedom is defined by the data structure:

                                                            {r, H,q,p} 

where 
   r = k-vector whose elements define the number of degrees of freedom for each simple joint, with n = r1 + úú + rk. 
   H = [H1 .. Hk], a matrix composed of the k joint map matrices of the simple joints.
   q = n-vector of joint coordinate names.
   p = n-vector of joint quasi-velocity names.  

Let us consider some elementary examples.

A two degree of freedom revolute joint as a simple joint (spherical joint):

r1={2};H1={{1,0},{0,0},{0,1},{0,0},{0,0},{0,0}};
q1={a1x,a1z};p1={w1x,w1z};

A two degree of freedom revolute joint as a compound joint (universal joint):

r2={1,1};H2={{1,0},{0,0},{0,1},{0,0},{0,0},{0,0}};
q2={a2x,a2z};p2={w2x,w2z};

The function Joints computes all of the necessary quantities for modeling a list of joints. These include the kinematic matrix V(q), the Euclidean configuration matrix X(q), and the overall joint map matrix H(q). These are reported as lists.

JointLst={{r1,H1,q1,p1},{r2,H2,q2,p2}};

{V,X,H}=Joints[JointLst];

Print[V]

Print[X]

Print[H]

Notice that joints 1 and 2 produce the same Euclidean configuration matrix. Thus, the configuration coordinates have the same physical meaning. Other observations are: (i) the overall joint map matrix for a compound joint depends on the joint configuration parameters, and (ii) the kinematic matrix is simpler for a compound joint representation than for a simple joint representation. A free or unconstrainedbody would be modeled using a six degree of freedom joint, such as

r3={6}; H3=IdentityMatrix[6];
q3={ax,ay,az,x,y,z}; p3={wx,wy,wz,ux,uy,uz};

JointLst={{r3,H3,q3,p3}};

{V,X,H}=Joints[JointLst];

Print[V]

Print[X]

Print[H]

Rigid Bodies

A rigid body is defined by its mass, inertia matrix and the location of distinguished points or nodes where joints or sensors may be located. We assume the following:

1) There is a distinguished point that corresponds to the inboard joint of the body. The body frame has its origin located there. 
2) The center of mass and all other points of interest (nodes) including outboard joint locations are defined in the body frame.
3) The inertia matrix is defined in the body frame and it is the inertia matrix about the center of mass.

The data for a rigid body is organized in a list as follows> A rigid body with k outboard nodes is defined by the data structure:
 
                                    {com,{out1,..,outk},m,Inertia}
where 
   com is the center of mass location, 
   outi = {node number, location} for the ith outboard node, 
   m is the mass, and 
   Inertia is the inertia tensor (about the center of mass).

Flexible Bodies

TSi Dynamics accommodates flexible bodies which satisfy the following assumptions:

    • Local deformations are small, so that linear stiffness (quadratic strain energy) and dissipation (quadratic dissipation function) relationships apply.
    • Body deformations can be characterized by a finite set of deformation coordinates, denoted by x. 
    • Body frame center of mass location and joint locations and orientations can be defined as (affine) linear functions of the deformation coordinates.
     
Any flexible body model in which a modal representation of flexure is valid satisfies these assumptions. Even with these assumptions, large global deformations are possible in which case the body inertia matrix as represented in the primary body frame may be a function of the deformation coordinates.

The primary body reference frame is a frame fixed in an infinitesimal body element with origin at the inboard joint location. We use the following notation for specifying the location of points in the primary frame under deformation. Let x be the n-vector of deformation coordinates. From the linearity assumption the position y of any point P in the body under deformation can be characterized by a matrix C and the relation

                                                   y = C.{1, x}

where C is a 3x(n+1) matrix. This is the type of information required to locate the center of gravity in the deformed body.

More information is typically required for general nodes at which sensors or joints are located. Here it is necessary to determine both the position of the node and the angular orientation of a local reference frame fixed in an infinitesimal element at the node. Under the small deformation assumption the orientation can be specified by any set of local parameters. We use Euler angles (3-2-1 convention). The orientation and position of the frame is defined by a 6-vector Y, the first three elements being the angles and the last three the position vector all relative to and specified in the primary body frame. As above, Y can be characterized by a matrix C

                                                  Y = C.{1,x}

where C is a 6x(n+1) matrix.

A flexible body with k outboard nodes and n deformation coordinates is defined by the data structure:

                       {Ccom,{out1,..,outk},m,{M(x),B,K},x,v}

where Ccom is a 3x(n+1) matrix that defines the center of mass location, 
   outi = {node number, Couti}, where Couti is a 6x(n+1) matrix that defines the orientation and location of the ith outboard node
   m is the mass 
   M(x) is the inertia matrix, B is the dissipation matrix, K is the stiffness matrix
   x is an n-vector of deformation coordinate names
   v is an n-vector of deformation velocity names

A First Example: The Double Pendulum

Building the Equations of Motion

Building a model in TSi Dynamics generally proceeds in the following way:

    Define joint data.
    Define body data.
    Define interconnection structure.
    Define potential energy function and generalized forces.
    Assemble the model.

The following is a simple example.

Define Joint Data

r1={1};H1={{1},{0},{0},{0},{0},{0}};
q1={a1x};p1={w1x};

r2={1};H2={{1},{0},{0},{0},{0},{0}};
q2={a2x};p2={w2x};

JointLst={{r1,H1,q1,p1},{r2,H2,q2,p2}};

Define Body Data

com1={0,0,-l1}; mass1=m1; out1={2,{0,0,-l1}};
Inertia1=DiagonalMatrix[{0,0,0}];

com2={0,0,-l2}; mass2=m2; out2={3,{0,0,-l2}};
Inertia2=DiagonalMatrix[{0,0,0}];

BodyLst={{com1,{out1},mass1,Inertia1},{com2,{out2},mass2,Inertia2}};

Define Interconnection Structure

TreeLst={{{1,1},{2,2}}};

Define Potential Energy

In this case only gravity contributes to the potential energy. The only generalized forces are external torques acting at the two joints.

g=gc; PE=0; Q={T1,T2};

{V,X,H,M,Cp,Fp,p,q}=CreateModelSim[JointLst,BodyLst,TreeLst,g,PE,Q];

We can look at some results.

Print[MatrixForm[V]]

Recall that V returns as a list of kinematic matrices Ä one for each joint. Hence in this case we get two 1x1 matrices. This can be assembled into a single block diagonal matrix but it is more efficient to retain the list form.

Print[MatrixForm[M]]

Print[MatrixForm[Fp]]

When doing time consuming computations you may want to save results periodically and then reload them for later computations. Since, as an analysis proceeds, different data sets may be saved under the same file name, it is a good idea to delete any older file before the new one is saved.

DeleteFile["dbl_pend.dat"];
Save["dbl_pend.dat",p,q,V,X,H,Cp,Fp,M];

Of course, the data is reloaded by the standard:

<<dbl_pend.dat;

Simulation: Numerical Solution

We can easily construct a simulation within Mathematica. First assemble the system parameter matrices computed above into a system of ordinary differential equations using the TSi Dynamics function MakeODEs.

Equations=MakeODEs[p,q,V,M,Cp,Fp,t]

Now, replace parameter symbols by numbers and set up initial conditions.

DataReplacements={m1->1,l1->2,m2->1,l2->1,gc->1,T1->0,T2->0}

Equations1=Equations/.DataReplacements

vars=Join[p,q]

InitialConditions=Table[ToExpression[StringJoin[ToString[vars[[i]]],"[0]==0"]],{i,Length[vars]}]

InitialConditions={w1x[0] == 0, w2x[0] == 0, a1x[0] == .1, a2x[0] == .2}

Equations2=Join[Equations1,InitialConditions]

Finally, integrate the equations using the Mathematica function NDSolve and plot the results.

sols=NDSolve[Equations2,vars,{t,0,45}]

Plot[Evaluate[a1x[t]/.sols],{t,0,45}]

Creating MEX-Files

It may be desireable to use the system model in an external program. TSi Dynamics provides functions to create C-Code which compiles as a MEX-File for use as an S-function in SIMULINK. An S-function may have inputs, outputs so that it can be interconnected with other subsystems within SIMULINK and it can have parameters which can be defined from within SIMULINK. In the following example, we define the joint torques as inputs and the (y,z) position of mass 2 as the outputs. Parameters include the system two masses, the two lengths and the gravitational constant.

<<dbl_pend.dat;

We can use the TSi Dynamics function EndEffector to define the functions.

ChainLst={{1,1},{2,2}};
TerminalNode=3;
Xout=EndEffector[ChainLst,TerminalNode,BodyLst,X];
yout={Xout[[2,4]],Xout[[3,4]]}

Inputs = Q;
Outputs = yout;
MEXFilename = "dbl_pend.c";
paramvec={m1,l1,m2,l2,gc};
PassedParams = {"X0","m1","l1","m2","l2","gc"};
PassedParamsDimensions = {{2*Length[p],1},{1,1},{1,1},{1,1},{1,1},{1,1}};
DeleteFile["dbl_pend.dat"];

Save["dbl_pend.dat",p,q,Inputs,Outputs,PassedParams,
                 PassedParamsDimensions,V,X,H,Cp,Fp,M,MEXFilename];

Building More Complex Models

Introductory Remarks

TSi Dynamics provides several routes to developing models. The most direct is to use the functions CreateModel or CreateModelSim. These two functions work essentially the same way but produce models in slightly different form. The latter is almost always preferable for creating simulations. One example of its use has already been given above for creating a model of the double pendulum. The following more elaborate examples illustrate some other functions in TSi Dynamics.

Partial Suspension System

We consider a simple split axle suspension in planar motion. The system consists of a chassis, two axles, two tires, two springs Ä each connecting nodes on one axle and the chassis, and two torsional dampers. The modeling process begins in the standard way with joint and body definitions. The coordinate system is chosen with z down and motion takes place in the y-z plane (y points to the right).

Joint Data

(* Joint 1 - planar motion *)
r1={3};
q1={theta1,y,z};p1={w1,vy,vz};
H1={{1,0,0},{0,0,0},{0,0,0},{0,0,0},{0,1,0},{0,0,1}};
(* Joint 2 - suspension *)
r2={1};
q2={theta2};p2={w2};
H2={{1},{0},{0},{0},{0},{0}};
(* Joint 3 - suspension *)
r3={1};
H3={{1},{0},{0},{0},{0},{0}};
q3={theta3};p3={w3};

JointLst={{r1,H1,q1,p1},{r2,H2,q2,p2},{r3,H3,q3,p3}};
{V,X,H}=Joints[JointLst];

Body Data

(* Body 1 - Chassis *)
cm1={0,0,0};
out1={{2,{0,b/2,a/2}},{3,{0,-b/2,a/2}},{6,{0,b/2+a/2,-a/2}},{7,{0,-b/2-a/2,-a/2}}};
I1=DiagonalMatrix[{Ixx,Iyy,Izz}];

(* Body 2 - Wheel & Axle *)
cm2={0,d,0};
out2={{4,{0,d,R}},{8,{0,Sqrt[2]*a,0}}};
I2=DiagonalMatrix[{Wxx,Wyy,Wzz}];

(* Body 3 - Wheel & Axle *)
cm3={0,-d,0};
out3={{5,{0,-d,R}},{9,{0,-Sqrt[2]*a,0}}};
I3=DiagonalMatrix[{Wxx,Wyy,Wzz}];

BodyLst={{cm1,out1,m1,I1},{cm2,out2,mw,I2},{cm3,out3,mw,I3}};

System Structure

Notice that this system has a tree structure with two chains.

TreeLst={{{1,1},{2,2}},{{1,1},{3,3}}};

Potential Energy Constructions

The potential energy comes from two sources in addition to the gravitational potential: (i) the tires, and (ii) the springs. The tires are leaves of the tree and they react with the space frame. Computing the potential energy of the tires illustrates the use of the function LeafPotential.

It is assumed that each tire potential energy function is defined in terms of the space coordinates of the tire contact point (the tire outboard node). In this example, we suppose that only the vertical displacements are significant.

(* Define Tire PE in terms of Space Coordinates *)

yy={{x2,y2,z2},{x3,y3,z3}};
Pot=(k/2)*(z2-z20)^2+(k/2)*(z3-z30)^2;

Note that z20 and z30 may be specified as externally defined functions of time to simulate riding over rough terrain. LeafPotential replaces the absolute coordinates by the generalized coordinates.

(* Obtain Tire PE in terms of generalized Coordinates *)

TirePE=LeafPotential[BodyLst,TreeLst,X,H,Pot,yy];

TirePE

The spring potential energy function SpringPotential is a specialized function for springs connecting two nodes in the system. It is assumed that the potential energy is a function of the length of the spring. Because only one variable is involved it is convenient to specify the potential energy as a function rather than as an expression in which the variables must be separately identified. SpringPotential constructs the potential energy function for the spring in terms of the generalized system coordinates.

SpringPot[x_]:=(kdel/2)*(x-ls)^2
q=Join[q1,q2,q3];
SpringPE1=SpringPotential[6,8,SpringPot,TreeLst,BodyLst,X,q];
SpringPE2=SpringPotential[7,9,SpringPot,TreeLst,BodyLst,X,q];

PE=TirePE+SpringPE1+SpringPE2;

SpringPE1

SpringPE2

Generalized Force

The only generalized forces acting are due to the dampers. We assume that the forces due to the torsional dampers are proportional to the relative angular velocities of the axles with respect to the chassis. As an alternative one could place a linear damper between two nodes and use the function DamperForce to compute the generalized force.

Q={0,0,0,-c*w2,-c*w3};

Gravity

Since, we have chosen the z-axis to point down, we change the sign of the gravitational acceleration

g=-gc;

Obtain Equations of Motion

Notice that since the joint kinematic data V, X, and H has already been computed, we include this data in the argument list of CreateModelSim to avoid recomputing it.

{V,X,H,M,Cmat,F,pp,qq}=CreateModelSim[JointLst,BodyLst,TreeLst,g,PE,Q,V,X,H];

M

Simplify[Cmat.pp]

F

A Flexible Two Link Manipulator

Consider a two link, planar mechanism in which the lower end of link 1 is attached to the space frame by a revolute joint and the upper end of link 1 and the lower end of link 2 are joined by a second revolute joint. Assume that the links have the same mass per unit length, rho, modulus of elasticity, El, cross sectional area moments (constant along the length) {Ixx,Iyy,Izz}, and the lengths L1 and L2, respectively.

Modeling a Flexible Link

There are many approaches to modeling individual flexible elements. One view fixes a "primary body reference frame" in an infinitesimal element at the inboard joint. The location of any point in the deformed body are specified in this frame. As the body undergoes motion, it is necessary to track the location and orientation of the primary body frame as well as the location in that frame of any interesting body points. Thus, the number of degrees of freedom associated with a free deformable body will be the six degrees of freedom associated with the primary frame plus some number of degrees of freedom associated with the body deformations. Usually, it is desired to characterize body deformations by a finite number of degrees of freedom and this is often accomplished by using finite element or (truncated) modal expansion methods.      

As an example take one of the links in the two link manipulator. We fix a reference frame in the link at its inboard joint. Suppose that the beam lies along the z-axis (points upward) with the inboard joint at the bottom. Moroever, assume that the link bends only in the y-z plane. Within this frame the link deforms as a centelevered beam and we can characterize the beam deformations in terms of its modes.  A cantelevered (Euler) beam of length L has the characteristic equation

																				Cos[b*L] * Cosh[b*L] = -1

which has a countable set of solutions br, r = 1,2,3...; the first three roots are approximately:

bL={0.600*Pi,1.49*Pi,2.50*Pi};

They define the natural frequencies by the relation

w[i_]:=(bL[[i]])^2*(El*Ixx/(rho*L^4))^(1/2);

Correspondingly, the mode shapes are:

W[i_,z_]:=((Cosh[bL[[i]]*z/L]-Cos[bL[[i]]*z/L])/(Cosh[bL[[i]]]+Cos[bL[[i]]])

              - (Sinh[bL[[i]]*z/L]-Sin[bL[[i]]*z/L])/(Sinh[bL[[i]]]+Sin[bL[[i]]]));

From the above information we can construct, for example, a one-mode model for a link of length L. First, the inertia, damping and stiffness matrices are:

M={{Ia,0,0,0,-m,0,0},
   {0,Ib,0,m,0,0,a},
   {0,0,Ic,0,0,0,0},
   {0,m,0, m,0,0,b},
   {-m,0,0,0,m,0,0},
   {0,0,0,0,0,m,0},
   {0,a,0,b,0,0,1}};

B=0*IdentityMatrix[7]; 
K=DiagonalMatrix[{0,0,0,0,0,0,(0.600*Pi)^2}]*(El*Ixx/(m*L^3));

where Ia, Ib, Ic are respectively Ixx m L^2/A, Iyy m L^2/A, Izz m L^2/A. Let y1 denote the modal coordinate, then the y position of the right end and anglular rotation about the x-axis (phiend) under deformation are given by:

yend=W[1,L]*y1; phiend=D[W[1,z],z]/.z->L;

and the center of mass (which corresponds to the midpoint) is

ycom=W[1,L/2]*y1;

Consequently, we have the matrices Cout and Ccom:

Cout={{0,D[W[1,z],z]/.z->L},
      {0,0},{0,0},
      {0,0},{0,W[1,L]},{L,0}};
Ccom={{0,0},{0,W[1,L/2]},{L,0}};

The link body data is then

Link1={N[Ccom],{{2,N[Cout]}},m,{M,B,K},{y11},{v11}}/.L->L1

Link2={N[Ccom],{{3,N[Cout]}},m,{M,B,K},{y21},{v21}}/.L->L2

BodyList={Link1,Link2};

r1={1};
H1=Transpose[{{1,0,0,0,0,0}}];
q1={theta1};p1={w1};

r2={1};
H2=Transpose[{{1,0,0,0,0,0}}];
q2={theta2};p2={w2};

JointList={{r1,H1,q1,p1},{r2,H2,q2,p2}};

TreeList={{{1,1},{2,2}}};
PE=0; Q={0,0,0,0};

{V,X,H,M,Cp,Fp,p,q}=CreateModelSim[JointList,BodyList,TreeList,g,PE,Q];

We can examine some of the results:

M

Simplify[Cp.p]

Fp

The absolute position and orientation of a frame at the end of the second link can be obtained with the function EndEffector.

EndPos=EndEffector[3,TreeList,BodyList,X];

Let's take a look at the position vector:

EndPos[[{1,2,3},4]]

Constraints

Introductory Remarks

Many systems of interest do not have a simple tree structure but can be viewed in terms of a tree with additional algebraic and/or differential constraints. A key advantage of the Lagrange formulation is that the explicit formulation of the relations that constrain the configuration of a system is usually avoided by appropriate choice of the generalized coordinates. However, with complex systems, a correct choice of generalized coordinates may be far from obvious. 

One systematic method for dealing with mechanical configurations containing closed loops is to break (free) a minimal number of joints so as to obtain a tree structure which, of course, has a natural set of configuration coordinates. The loop-breaking joint constraints are then added as suplemental algebraic relations. The number of degrees of freedom of the system is typically reduced by the number of independent relations. If these relations are (globally) solvable for a like number of coordinates then these coordinates can be explicitly eliminated. When there are k such constraints, the number of coordinates is reduced by k and the number of quasi-velocities is also reduced by k. It may not be possible to eliminate all k excess coordinates, but it is always possible to eliminate k quasi-velocities.

Some chain and loop structures give rise to differential constraints Ä such as when rolling of one body on another occurs. In general, when k independent differential constraints are present it is possible to eliminate k quasi-velocities. If the constraints are "partially integrable" then, in principle, it is also possible to eliminate some coordinates as well. When the constraints are "completely integrable" the number of configuration degrees of freedom are reduced by k and the system is said to be "holonomic." Otherwise it is "nonholonomic."

Even when it is not possible or practical to eliminate all k configuration coordinates it is fairly straightforward to derived the reduced set of equations in which the k quasi-velocities are removed. These equations form a valid, closed system of differential equations that completely describe the system behavior. Presently, TSi Dynamics provides functions to generate this latter set of equations unless the user specifies a set of configuration coordinates that are to be removed.

A Simple Closed Chain

Consider a planar two bar linkage in which the lower end of bar 1 is constrained by a revolute joint on the y-axis at y=0 and the upper end of bar two slides on the z-axis. The upper end of bar 1 and the lower end of bar 2 are connected by a revolute joint.

r1={1};
H1=Transpose[{{1,0,0,0,0,0}}];
q1={theta1};p1={w1};

r2={1};
H2=Transpose[{{1,0,0,0,0,0}}];
q2={theta2};p2={w2};

JointList={{r1,H1,q1,p1},{r2,H2,q2,p2}};

Mass1=m1; CenterOfMass1={0,0,L1/2}; OutboardNode1={2,{0,0,L1}};
InertiaMatrix1=DiagonalMatrix[{m1*L1^2/12,m1*L1^2/12,0}];

Mass2=m2; CenterOfMass2={0,0,L2/2}; OutboardNode2={3,{0,0,L2}};
InertiaMatrix2=DiagonalMatrix[{m2*L2^2/12,m2*L2^2/12,0}];

BodyList={{CenterOfMass1,{OutboardNode1},Mass1,InertiaMatrix1},
          {CenterOfMass2,{OutboardNode2},Mass2,InertiaMatrix2}};

TreeList={{{1,1},{2,2}}};
PE=0; Q={0,0};

{V,X,H,M,Cp,Fp,p,q}=CreateModelSim[JointList,BodyList,TreeList,g,PE,Q];

ChnBodyList={{CenterOfMass1,{0,0,L1},Mass1,InertiaMatrix1},
             {CenterOfMass2,{0,0,L2},Mass2,InertiaMatrix2}};

EndPos=EndEffector[ChnBodyList,X];

G=EndPos[[{2},4]]

The the relation G=0 is not solvable for either theta1 or theta2 for all values of L1 and L2 Ä although, as we will describe below, it can be solved for particular values. Thus, we will not try to eliminate any configuration coordinate. The constrained dynamics are obtained with:

{Mm,Cm,Fm,Vm,phat,qhat}=AlgConstrainedSys[M,Cp,Fp,V,G,p,q]

These parameters define the constrained system dynamics. Now, let's consider the special case: L1=L2=L. This is the only situation where G=0 can be solved for either theta1 or theta2, so that either angle can be used to parameterize the system configuration. For simplicity, also set m1=m2=m. In this case the system matrices are:

{Mm,Cm,Fm,Vm}=Simplify[{Mm,Cm,Fm,Vm}/.{L1->L,L2->L,m1->m,m2->m}]

It is convenient toassemble the corresponding differential equations.

Eqns=MakeODEs[phat,qhat,Vm,Mm,Cm,Fm,t]

Notice that there are two kinematic equations and one dynamic equation. Since L1=L2=L it is easy to show that theta2=-2 theta1. Hence we can make this replacement to obtain:

SpEqns=Simplify[Eqns/.theta2[t]->-2*theta1[t]]

Simple trigonometric identities can be used to verify that these equations are equivalent to those given by Ginsberg for this example.

Another Closed Chain

Consider a two bar linkage in the x-y plane. The left end of bar 1 is constrained to slide on the x-axis and the right end of bar two slides on the y-axis. The right end of bar 1 and the left end of bar 2 are connected by a revolute joint.

r1={1,1};
H1=Transpose[{{0,0,1,0,0,0},{0,0,0,1,0,0}}];
q1={theta1,x};p1={w1,vx};

r2={1};
H2=Transpose[{{0,0,1,0,0,0}}];
q2={theta2};p2={w2};

JointList={{r1,H1,q1,p1},{r2,H2,q2,p2}};

Mass1=m1; CenterOfMass1={L1/2,0,0}; OutboardNode1={2,{L1,0,0}};
InertiaMatrix1=DiagonalMatrix[{0,m1*L1^2/2,m1*L1^2/2}];

Mass2=m2; CenterOfMass2={L2/2,0,0}; OutboardNode2={3,{L2,0,0}};
InertiaMatrix2=DiagonalMatrix[{0,m2*L2^2/2,m2*L2^2/2}];

BodyList={{CenterOfMass1,{OutboardNode1},Mass1,InertiaMatrix1},
          {CenterOfMass2,{OutboardNode2},Mass2,InertiaMatrix2}};

TreeList={{{1,1},{2,2}}};
PE=0; Q={0,0,0};

{V,X,H,M,Cp,Fp,p,q}=CreateModelSim[JointList,BodyList,TreeList,g,PE,Q];

DeleteFile["Linkage.dat"];

Save["Linkage.dat",p,q,V,X,H,Cp,Fp,M];

<<Linkage.dat;

ChnBodyList={{CenterOfMass1,{L1,0,0},Mass1,InertiaMatrix1},
             {CenterOfMass2,{L2,0,0},Mass2,InertiaMatrix2}};

EndPos=EndEffector[ChnBodyList,X];

G=EndPos[[{1},4]]

Since G=0 always has a solution for x, we will eliminate x from the constrained equations.

{Mm,Cm,Fm,Vm,phat,qhat}=AlgConstrainedSys[M,Cp,Fp,V,G,p,q,{2}]

Simplify[Cm.phat]

The Rolling Disk: a Classic Example of Differential Constraints

As an illustration we take an example from Neimark and Fufaev. Consider a disk that rolls without slipping on the x-y plane (z is up). Assume that the disk is of mass m and radius R. One approach is to ignore the the rolling constraints and formulate the equations of motion for the free disk in space, and then add the required constraints. A body frame is established with origin at the center of the disk. The six degree of freedom, simple joint is defined by:

r={6};
H=IdentityMatrix[6]

q={psi,theta,phi,x,y,z};
p={wx,wy,wz,vx,vy,vz};
JointList={{r,H,q,p}};

In this setup the x,y,z coordinates locate the center of the disk in the space frame and psi, theta, phi are Euler angles in the z-x-y (or 3-2-1) convention. 

The body data is:

Mass=m; CenterOfMass={0,0,0}; OutboardNode={2,{0,-R*Sin[psi],-R*Cos[psi]}};
InertiaMatrix=DiagonalMatrix[{J,Iy,Iy}];
BodyList={{CenterOfMass,{OutboardNode},Mass,InertiaMatrix}};

and the remaining data:

TreeList={{{1,1}}};
PE=0; Q={0,0,0,0,0,0};

The unconstrained disk model is obtained with

{V,X,H,M,Cp,Fp,p,q}=CreateModelSim[JointList,BodyList,TreeList,g,PE,Q];

DeleteFile["disk2.dat"];

Save["disk2.dat",p,q,V,X,H,Cp,Fp,M];

Now, we formulate the differential constraints. Rolling without slipping implies that the velocity of the disk contact point must be zero. An expression for the (angular and translational) velocity at the outboard node as a function of the configuration variables is easily obtained using the function EndEffectorVelocity.

<<disk2.dat;

{ChainList}=TreeList;
VCont=EndEffectorVelocity[ChainList,BodyList,X,H,p]

The rolling constraint requires that the translation velocity of the contact point (the last three of the elements of the six-vector VCont) be zero. Thus, we obtain the constrained system:

{Mm,Cm,Fm,Vm,T,phat}=DiffConstrainedSys[M,Cp,Fp,V,VCont[[Range[4,6]]],p,q,{4,5,6}];

The dynamics are reduced to three dimensions and the original six quasi-velocities are reduced to three, in fact, we have:

phat

The set of configuration coordinates is not reduced by the function DiffConstrainedSys. In general, a set of differential constraints may not admit any such reduction. Such would be the case if the constraints were completely nonholonomic. In the present case, however, the constraints are "partially" integrable and from basic geometry one can see that (cg) hight = Rcos(theta)-R. Using this relationship, the coordinate z can be eliminated from the equations.

{Sols}=Solve[{X[[1]][[3,4]]==R*(Cos[theta]-1)},{z}]

Because the translation parameters x,y,z are, in fact, the space coordinates, z is precisely the hight of the center of gravity. Inspection shows that only Mm and Vm depend on z so we define:

Mmo=Simplify[Mm/.Sols];

Vmo=Simplify[Vm[[ {1,2,3,4,5} ]]/.Sols];

Now, we assemble the governing equations.

Eqns=MakeODEs[phat,q[[{1,2,3,4,5}]],Vmo,Mmo,Cm,Fm,t]

Save["disk2.dat",p,q,V,X,H,Cp,Fp,M,Eqns];

Since x,y,z locate the center of gravity it is of interest to compute the contact point.

ChnBodyList={{CenterOfMass,{0,-R*Sin[psi],-R*Cos[psi]},Mass,InertiaMatrix}};

{ContactPoint}=Simplify[Transpose[EndEffector[ChnBodyList,X][[{1,2,3},{4}]]]/.Sols]

Consider the equations of motion as given by the list Eqns. Careful inspection of the equations suggests that representation of the angular velocity in a fame that does not rotate about the body-x axis may simplify them. Thus, we transform the angular velocity coordinates, wy, wz, via the relations:

a1=wyy[t]==Cos[psi[t]]*wy[t]-Sin[psi[t]]*wz[t]
a2=wzz[t]==Sin[psi[t]]*wy[t]+Cos[psi[t]]*wz[t]

The inverse transformation rules are obtained using Mathematica's Solve function

Trans=Simplify[Flatten[Solve[{a1,a2},{wy[t],wz[t]}]]]

<<disk2.dat;

The transformed equations are obtained with the function StateTransformation:

StateVars={psi,theta,phi,x,y,wx,wy,wz}

TrEqns=StateTransformation[Eqns,StateVars,Trans,{wyy,wzz},t]

These equations are obviously equivalent to thedisk equations given by Neimark and Fufaev. Another derivation is given by Meirovitch. To compare our equations with his, it is necessary to reduce them to second order form by eliminating the quasi-velociites (thus, we get Lagrange's equations), and to perform a minor transformation of anglular coordinates:

Rules1=Flatten[Solve[TrEqns[[{1,2,8}]],{wx[t],wyy[t],wzz[t]}]]

Rules2={wx'[t]->D[wx[t]/.Rules1,t],
        wyy'[t]->D[wyy[t]/.Rules1,t],
        wzz'[t]->D[wzz[t]/.Rules1,t]}

Rules3={theta[t]->th[t]-Pi/2,
        theta'[t]->th'[t],
        theta''[t]->th''[t]}

Rules4={phi[t]->-pi/2+ph[t],
        phi'[t]->ph'[t],
        phi''[t]->ph''[t]}

LagEqns=Simplify[TrEqns[[{5,6,7}]]/.Rules1/.Rules2/.Rules3/.Rules4]

These equations are easily confirmed to be equivalent to those given by Meirovitch.

Timing Considerations

Introductory Remarks

The computer time required for building a model can vary greatly depending on the computer and its configuration. But it is also very much dependent on the system being analysed and the way in which the system parametric data is defined. Excessive execution time can sometimes arise from suprisingly minor factors in the formulation, whereas obvious indicators of system complexity such as the number of degrees of freedom may be of little consequence. The number of symbolic parameters used in the data definition can be significant, but even more important is how those parameters are introduced. 

The following paragraphs illustrate some of these considerations and show the value of simplifying the kinematic expressions by replacing recurrent expressions by temporary variables. All of the times shown below were obtained using a Macintosh Quadra 800 with 24 MB od Ram.

Replacements for Repeated Expressions

When a set of expressions contain many repetitions of of some subexpression then it makes sense to introduce a temporary variable for it and perform any required operations on that expression only once. This is, of course, a standard procedure for making numerical computation more efficient. TSi Dynamics includes a function, KinematicReplacements, which provides a limited extension of this idea to symbolic model assembly.  

KinematicReplacements analyzes the combined set of expressions in the lists V,X,H (after the joint kinematics are computed) for repeated expressions and replaces them with temporary variables. In addition to the new V,X,H lists, a set of inverse replacement rules are returned which allows the original V,X,H to be recovered whenever desired. Meanwhile, computations can be performed using the "simplified" kinematic expressions. The repeated expressions that are replaced are of two types: (i) expressions that depend on configuration coordinates as well as symbolic system data parameters, and (ii) expressions that depend only on symbolic system data parameters. Expressions that depend on system coordinates must be reintroduced before certain operations (such as differentiation with respect to a coordinate) are performed. To facilitate this, the inverse replacement rules are returned in two sets: rules1 and rules2. The application of rules1 restores all expressions containing the configuration coordinates. rules2 must be applied after rules1 to restore all parametric expressions.

Normally, the function KinematicReplacements would be employed after the computation of V,X,H using Joints. The simplified V,X,H are then used throughout the calculation of generalized forces and potential energy functions. The coordinated dependent expressions are then restored by applying rules1 before using the functions CreateModel or CreateModelSim. Then the parameter dependent expressions can be restored by applying rules2.

This process will have marginal affect on simple problems and is likely to produce slightly longer computation times. But it can have a dramatic affect on computation times for more complex problems.

Partial Suspension System Revisited

To see how this process works we will take another look at the Partial Suspension System example describe above. The computations are performed in two ways: without and then with the use of  KinematicReplacements as described in the preceeding paragraphs. Timing data is obtained in both cases and it will be seen that the use of the  KinematicReplacements strategy slightly slows things down.

Joint Data

(* Joint 1 - planar motion *)
r1={3};
q1={theta1,y,z};p1={w1,vy,vz};
H1={{1,0,0},{0,0,0},{0,0,0},{0,0,0},{0,1,0},{0,0,1}};
(* Joint 2 - suspension *)
r2={1};
q2={theta2};p2={w2};
H2={{1},{0},{0},{0},{0},{0}};
(* Joint 3 - suspension *)
r3={1};
H3={{1},{0},{0},{0},{0},{0}};
q3={theta3};p3={w3};

JointLst={{r1,H1,q1,p1},{r2,H2,q2,p2},{r3,H3,q3,p3}};
{V,X,H}=Joints[JointLst];

Body Data

(* Body 1 - Chassis *)
cm1={0,0,0};
out1={{2,{0,b/2,a/2}},{3,{0,-b/2,a/2}},{6,{0,b/2+a/2,-a/2}},{7,{0,-b/2-a/2,-a/2}}};
I1=DiagonalMatrix[{Ixx,Iyy,Izz}];

(* Body 2 - Wheel & Axle *)
cm2={0,d,0};
out2={{4,{0,d,R}},{8,{0,Sqrt[2]*a,0}}};
I2=DiagonalMatrix[{Wxx,Wyy,Wzz}];

(* Body 3 - Wheel & Axle *)
cm3={0,-d,0};
out3={{5,{0,-d,R}},{9,{0,-Sqrt[2]*a,0}}};
I3=DiagonalMatrix[{Wxx,Wyy,Wzz}];

BodyLst={{cm1,out1,m1,I1},{cm2,out2,mw,I2},{cm3,out3,mw,I3}};

System Structure

TreeLst={{{1,1},{2,2}},{{1,1},{3,3}}};

Potential Energy Constructions

(* Define Tire PE in terms of Space Coordinates *)

yy={{x2,y2,z2},{x3,y3,z3}};
Pot=(k/2)*(z2-z20)^2+(k/2)*(z3-z30)^2;

(* Obtain Tire PE in terms of generalized Coordinates *)

TirePE=LeafPotential[BodyLst,TreeLst,X,H,Pot,yy];

TirePE

SpringPot[x_]:=(kdel/2)*(x-ls)^2
q=Join[q1,q2,q3];
SpringPE1=SpringPotential[6,8,SpringPot,TreeLst,BodyLst,X,q];
SpringPE2=SpringPotential[7,9,SpringPot,TreeLst,BodyLst,X,q];

Timing[SpringPotential[6,8,SpringPot,TreeLst,BodyLst,X,q]]

PE=TirePE+SpringPE1+SpringPE2;

SpringPE1

SpringPE2

Generalized Force

Q={0,0,0,-c*w2,-c*w3};

Gravity

g=-gc;

Obtain Equations of Motion

Timing[CreateModelSim[JointLst,BodyLst,TreeLst,g,PE,Q,V,X,H]]

{V,X,H,rules1,rules2}=KinematicReplacements[V,X,H,q];

rules1

rules2

(* Define Tire PE in terms of Space Coordinates *)

yy={{x2,y2,z2},{x3,y3,z3}};
Pot=(k/2)*(z2-z20)^2+(k/2)*(z3-z30)^2;

(* Obtain Tire PE in terms of generalized Coordinates *)

TirePE=LeafPotential[BodyLst,TreeLst,X,H,Pot,yy];

SpringPot[x_]:=(kdel/2)*(x-ls)^2
q=Join[q1,q2,q3];
SpringPE1=SpringPotential[6,8,SpringPot,TreeLst,BodyLst,X,q];
SpringPE2=SpringPotential[7,9,SpringPot,TreeLst,BodyLst,X,q];

Timing[SpringPotential[6,8,SpringPot,TreeLst,BodyLst,X,q]]

PE=TirePE+SpringPE1+SpringPE2;

SpringPE1

SpringPE2

{PE,Q,V,X,H}={PE,Q,V,X,H}/.rules1;

Timing[CreateModelSim[JointLst,BodyLst,TreeLst,g,PE,Q,V,X,H]/.rules2]

Tractor Semi-Trailer

An example which benefits significantly from the use of replacement rules is a tractor semi-trailer. The goal of a recent study was to show how deformable tires, tire parameters, and front end alignment parameters affect turning radius of a tractor semi-trailer. This example uses a simple model for a deformable pneumatic tire, namely, a single deformation parameter (proportional to sideslip angle), so that the cornering force is given by

                             F=-kappa*arctan[vy/vx]

Only the computation of the right front tire cornoring force is given below.

Main body (tractor) joint

(* Joint 1  *)
r1={3}; q1={theta,x,y}; p1={wth,vx,vy};
H1={{0,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1},{0,0,0}};

Tractor front end joints

Define castor sr (right) and sl (left) and camber tr (right) and tl (left). Note that positive castor and camber implies spindle points down and forward and inward.

(* Joint 2 right front spindle *)
r2={1}; q2={delta2}; p2={wdel2}; 
H2=Transpose[{{-Sin[sr],-Sin[tr],1,0,0,0}}];

(* Joint 3 left front spindle *)
r3={1}; q3={delta3}; p3={wdel3}; 
H3=Transpose[{{-Sin[sl],Sin[tl],1,0,0,0}}];

Tractor-Semi joint

(* Joint 4 5th wheel *)
r4={1}; q4={delta4}; p4={wdel4};
H4=Transpose[{{0,0,1,0,0,0}}];

JointLst={{r1,H1,q1,p1},{r2,H2,q2,p2},{r3,H3,q3,p3},{r4,H4,q4,p4}};

Body Data

Tractor body frame is located midway between rear axles on body centerline and at axle hight. It is right handed with x-forward, y-left, and z-up. 

	           wb = wheelbase (186")
	           R=tire radius (12")
	           2b=track (72-11/16")
	           pin= pin location from rear bogie (8")
	           d=tractor rear axle location from rear bogie (28")
	           tb=trailor base (388")
	           kappa=tire coefficient (approx 14,000 N?)

wb=186*.0254; R=12*.0254; b=(72+11/16)*.0254;
pin=8*.0254; d=28*.0254; tb=388*.0254;
Jtracxx=1.414*10^3; Jtracyy=5.684*10^4; Jtraczz=5.698*10^4;Jtracxz=1.96*10^3;
Jtrailxx=2.202*10^5; Jtrailyy=4.852*10^5; Jtrailzz=5.065*10^5;
mtrac=1.15*10^4; mtrail=4.376*10^4; mtire=40;
Ixx=mtire*(R^2)/2; Iyy=mtire*(3*R^2+(11*.0254)^2)/12; Izz=Iyy;

(* Body 1 tractor *)
cm1={2.4466,0,.9052};
out1={{2,{wb,-b,0}},{3,{wb,b,0}},{4,{pin,0,0}},{5,{d,-b,-R}},{6,{d,b,-R}},{7,{-d,-b,-R}},{8,{-d,b,-R}}};
I1={{Jtracxx,0,Jtracxz},{0,Jtracyy,0},{Jtracxz,0,Jtraczz}};

(* Body 2 right front wheel *)
cm2={0,0,0}; out2={{9,{sr*R,tr*R,-R}}};
I2=DiagonalMatrix[{Ixx,Iyy,Izz}];

(* Body 3 left front wheel *)
cm3={0,0,0}; out3={{10,{sl*R,-tl*R,-R}}};
I3=DiagonalMatrix[{Ixx,Iyy,Izz}];

(* Body 4 trailor *)
cm4={-5.77,0,1.0592}; 
out4={{11,{-tb,-b,-R}},{12,{-tb,b,-R}}};
I4=DiagonalMatrix[{Jtrailxx,Jtrailyy,Jtrailzz}];

BodyLst={{cm1,out1,mtrac,I1},{cm2,out2,mtire,I2},{cm3,out3,mtire,I3},{cm4,out4,mtrail,I4}};
TreeLst={{{1,1},{2,2}},{{1,1},{3,3}},{{1,1},{4,4}}};

q={theta,x,y,delta2,delta3,delta4};
p={wth,vx,vy,wdel2,wdel3,wdel4};

Computations Without Camber

tr=0.0; tl=0.0;

Joint Calculations

{V,X,H}=Joints[JointLst];

Tire Force Computations
For comparison, we compute the front right tire force two ways: (i) without replacements, and (ii) with replacements. In this case, there is an improvement from 166 sec to 4 sec.

(* Front right tire *)
ChnLst={{1,1},{2,2}};
Force={0,0,0,0,-kappa*ArcTan[v9y/v9x],0};
VelNames={w9x,w9y,w9z,v9x,v9y,v9z};
TerminalNode=9;
Timing[GeneralizedForce[ChnLst,TerminalNode,BodyLst,X,H,q,p,Force,VelNames]]

{V,X,H,rules1,rules2}=KinematicReplacements[V,X,H,q];

rules1

(* Front right tire *)
ChnLst={{1,1},{2,2}};
Force={0,0,0,0,-kappa*ArcTan[v9y/v9x],0};
VelNames={w9x,w9y,w9z,v9x,v9y,v9z};
TerminalNode=9;

Timing[GeneralizedForce[ChnLst,TerminalNode,BodyLst,X,H,q,p,Force,VelNames]]

Computations With Small Numerical Camber

tr=0.01; tl=0.01;

Joint Calculations

{V,X,H}=Joints[JointLst];

Tire Force Computations
As before, compute the front right tire forces. In this case, there is an improvement from 597 sec to 18 sec.

(* Front right tire *)
ChnLst={{1,1},{2,2}};
Force={0,0,0,0,-kappa*ArcTan[v9y/v9x],0};
VelNames={w9x,w9y,w9z,v9x,v9y,v9z};
TerminalNode=9;
Timing[GeneralizedForce[ChnLst,TerminalNode,BodyLst,X,H,q,p,Force,VelNames]]

{V,X,H,rules1,rules2}=KinematicReplacements[V,X,H,q];

rules1

(* Front right tire *)
ChnLst={{1,1},{2,2}};
Force={0,0,0,0,-kappa*ArcTan[v9y/v9x],0};
VelNames={w9x,w9y,w9z,v9x,v9y,v9z};
TerminalNode=9;
Q1=GeneralizedForce[ChnLst,TerminalNode,BodyLst,X,H,q,p,Force,VelNames];

Timing[GeneralizedForce[ChnLst,TerminalNode,BodyLst,X,H,q,p,Force,VelNames]]

Bibliography

To use TSi Dynamics requires some knowledge about dynamics. For those users seeking more information on the concepts behind the TSi Dynamics modeling formalism, we list a few reference books. Some have been cited above in our examples.

V. I. Arnold, Mathematical Methods of Classical Mechanics, Springer-Verlag: New York, 1989 (2nd ed.)

S. H. Crandall, D. C. Karnopp, E. F. Kurtz, D. C. Pridmore-Brown, Dynamics of Mechanical and Electromechanial Systems, McGraw-Hill: New York, 1968.

F. Gantmacher, Lectures in Analytical Mechanics, MIR: Moscow, 1975.

J. H. Ginsberg, Advanced Engineering Dynamics, Harper & Row: New York, 1988.

H. Goldstein, Classical Mechanics, Addison-Wesley: Reading, 1981 (2nd ed.)

D. T. Greenwood, Classical Mechanics, Prentice-Hall: Englewood Cliffs, 1977.

T. R. Kane and D. A. Levinson, Dynamics: Theory and Application, McGraw-Hill: New York, 1985.

C. Lanczos, The Variational Principles of Mechanics, University of Toronto Press: Toronto, 1970 (4th ed.)

L.Meirovitch, Methods of Analytical Dynamics, McGraw-Hill: New York, 1970.

Ju. I. Neimark and N. A. Fufaev, Dynamics of Nonholonomic Systems, American Mathematical Society: Providence, 1972.

R. M. Rosenberg, Analytical Dynamics of Discrete Systems, Plenum: New York, 1977.