Multibody Syst Dyn (2013) 30:153-171 DOI 10.1007/s11044-013-9352-0

Solvability of reactions in rigid multibody systems with redundant nonholonomic constraints

Marek Wojtyra • Janusz Fraczek

Received: 25 July 2012 / Accepted: 30 January 2013 / Published online: 5 March 2013 © The Author(s) 2013. This article is published with open access at Springerlink.com

Abstract The problem of calculating joint reaction forces in rigid body mechanisms with redundant constraints, both geometric and nonholonomic, is discussed. When constraint equations are dependent, some of the constraint reactions are unsolvable, i.e., cannot be uniquely determined using a rigid body model, whereas some others may be solvable. In this paper, analytic conditions, which must be fulfilled to obtain unique values of selected reaction forces in the presence of dependent nonholonomic constraints, are presented and proven. The concept of direct sum, known from linear algebra, is exploited. These purely mathematical conditions are followed by numerical methods that enable detection of constraints with uniquely solvable reactions. Similar conditions and methods were proposed earlier for holonomic systems. In this contribution, they are generalized to the case of linear nonholonomic constraints. An example of constraint reactions solvability analysis, for a mechanism subjected to redundant nonholonomic constraints, is presented.

Keywords Redundant constraints • Nonholonomic systems • Constraint reactions

1 Introduction

Dynamic analysis of a multibody system consists in calculating motion resulting from loads and driving constraints imposed on the system [10, 16, 28]. The time history of kinematic quantities, i.e., positions, velocities, and accelerations of bodies, is determined during the analysis. It is not necessary to calculate constraint reaction forces when the system is fric-tionless and its motion is the only object of interests [10]. Nevertheless, constraint reactions are frequently calculated along with the kinematic parameters describing motion, since they are needed, e.g., for design purposes or structural analyses.

M. Wojtyra (El) • J. Fraczek

Institute of Aeronautics and Applied Mechanics, Warsaw University of Technology, Nowowiejska 24,

00-665 Warsaw, Poland

e-mail: mwojtyra@meil.pw.edu.pl

J. Fraczek

e-mail: jfraczek@meil.pw.edu.pl

In multibody modeling, the constraint reactions not always can be uniquely determined. The problem of joint reactions indeterminacy, in engineering simulations of rigid multibody systems, is most often caused by redundant constraints presence. Redundant constraints are defined as constraints that can be removed without changing the kinematics of the system, i.e., they appear when the same degree of freedom is restricted more than once. Real-world mechanisms quite often are purposely overconstrained, mainly in order to strengthen or simplify their construction. Redundant constraints can be present in both holonomic and nonholonomic systems.

Multibody modeling of overconstrained systems is impeded. When redundant constraints are present in the mathematical model of holonomic multibody system, the constraint equations are dependent [10, 16]; hence, the constraint Jacobian matrix is rank-deficient. Therefore, special approaches are needed to handle equations of motion, e.g., [15, 20, 23, 27].

In the context of this article subject, special attention should be paid to the joint reactions calculation. In multibody modeling, sometimes it is assumed that constraint equations are independent (e.g., [4, 19]), hence redundant constraints, and thus reactions solvability problems, are simply ignored. Most often, however, redundant constraints are handled using one of two essentially different approaches. The first one consists in modifying the set of equations in order to exclude dependent equations [16, 28] (redundant constraints elimination is an example of this approach). The second one consists in leaving the set of equations unchanged and solving it using algorithms capable to deal with dependent equations (the minimum norm solution, sometimes based on the pseudoinverse technique [32], and the augmented Lagrangian formulation [2, 3] can serve as examples of this approach). In either case, purely mathematical operations, not supported by laws of physics, are performed to find the reactions [36]. Hence, the obtained solution usually reflects properties of the redundant constraints handling method rather than physical properties of the investigated multibody system. It should be noted, however, that contribution [11] indicates that in some cases combination of additional weighting factors (which are tuned to reflect elasticity of links and joints) with minimum norm solution leads to calculation of realistic reactions. So far, no general method of factors selection is proposed. Similarly, in the case of penalty methods, the authors of [13, 14, 30] formulated an idea that a flexibility aspect can be incorporated by relating the penalty factors to the real mechanism stiffness, which allows us to solve the indeterminacy problem without adding any major complication to equations of motion. However, so far no general rules on how to relate the penalty parameters with flexibility properties of multibody system parts and joints are provided.

A direct consequence of constraints redundancy is that some or all constraint reactions cannot be uniquely determined using a rigid body model [33]. Moreover, in the case of Coulomb friction in joints, the simulated motion might be nonunique [8]. It can be proven, however, that in the case of an overconstrained rigid body mechanism subjected to holo-nomic constraints, despite the fact that all constraint reactions cannot be uniquely determined, selected reactions or selected groups of reactions can be specified uniquely [33, 34]. In [34], methods of constraint reaction solvability analysis were proposed. An algebraic criterion, allowing detection of the joints for which reactions are uniquely solvable, was formulated. This criterion was followed by numerical methods for finding such joints.

In this contribution, systems subjected to both holonomic and linear nonholonomic constraints are investigated. It is shown that when redundancy problems and reaction forces are concerned, both types of constraints can be treated jointly and uniformly. The methods of constraint reactions solvability analysis are generalized to systems subjected to linear nonholonomic constraints. Firstly, the algebraic criterion is revisited and its applicability to

nonholonomic systems is considered. Then numerical methods of reactions solvability analysis are discussed. Finally, an exemplary system with redundant nonholonomic constraints is analyzed.

2 Equations of motion for constrained systems

2.1 Generalized coordinates and constraint equations

Consider a rigid multibody system described by n dependent generalized coordinates q. Holonomic constraints imposed on coordinates are represented by a set of m nonlinear scalar algebraic equations:

$ (q,t) = 0m x 1 • (1)

We assume that constraints (1) are consistent, i.e., no contradictory conditions are imposed on coordinates q.

The holonomic constraint equations (1) can be differentiated with respect to time to obtain constraints for velocities:

$ qq+ $t = 0mx1, (2)

where $ q is the constraint Jacobian matrix which may be rank-deficient and $ t is the vector of partial derivatives with respect to time.

Linear nonholonomic and scleronomic constraints (Pfaffian constraints) are represented by a set of h scalar equations, linear in velocities q [26]:

* (q)q = 0hx1. (3)

Note that in the case of redundant constraints presence, matrix * can be rank-deficient. It should be emphasized that the existence of a single nonintegrable velocity constraint does not necessarily mean that a system is nonholonomic, since this constraint may prove to be integrable by virtue of the remaining constraint equations. To check the holonomy of the system, the Frobenius theorem may be employed [5, 26]. It is important that our further considerations are valid for both integrable and nonintegrable Pfaffian constraints.

Equation (3) can be appended to Eq. (2) to obtain velocity level constraints for the whole system:

cq + b = o, (4)

where respective definitions of coefficient matrix C (often called a constraint matrix) and vector b are the following:

C(q,t) =

, b(q,t) =

(m+h)xn

We assume that constraints (4) are consistent, i.e., no contradictory conditions are imposed on velocities q. For our further considerations, it is important that when the mechanism is overconstrained, rows of matrix C are linearly dependent, hence

r = rank(C) < m + h. (6)

It is worth noting that constraint matrix C may be rank-deficient even when both 0 q and 9 are full-rank matrices.

The degree of redundancy p is defined as

p — m + h — r.

The velocity constraints (4) can be differentiated with respect to time in order to obtain constraints for accelerations:

with vector r defined as

r = —

Cq — r — 0(m+h)x1>

qq )qq+ t qq+ &tt

(* q )qq

In our further considerations, we will utilize absolute Cartesian coordinates [16, 28]. The main reason for choosing this type of coordinates is the fact that, in this approach, all kinematic pairs are treated in the same way, i.e., each pair is represented as a subset of Eq. (4). In the case of other types of coordinates, e.g., relative joint coordinates, the constraint equations are formulated only for selected, loop closing joints. The uniform treatment of all kinematic pairs simplifies calculation of joint reaction forces, however, it should be noted that the problem of joint reaction solvability—crucial for our considerations—is an issue of mechanism's structure, and from this point of view the choice of coordinates is irrelevant.

2.2 Equations of motion

Equations of motion for the considered multibody system can be written in the following form [10, 16]:

Mq + CTX — F(q, q,t),

where M denotes the mass matrix (which is invertible), F is the vector containing the external forces as well as the velocity dependent inertia terms, and X is the vector of Lagrange multipliers.

The Lagrange equations of the first kind (Eq. (10)) and the acceleration level constraints (Eq. (8)) can be written jointly in a matrix form, to obtain an index-1 set of differential-algebraic equations:

Equation (11) can be solved for the accelerations and the Lagrange multipliers (it can be shown that the Lagrange multipliers X are uniquely solvable only when matrix CM-1Cr is nonsingular, i.e., when matrix C has full row rank). In the case of constraints redundancy, however, matrix C is rank deficient, thus the leading matrix in Eq. (11) becomes singular, and Eq. (11) has a nonunique solution. More precisely, the solution for accelerations q is unique, and the solution for multipliers X is non-unique. To justify this statement, it is helpful to define matrix VnX(n-r), which is an orthogonal complement to C (which is now assumed

M CT" q F

C 0 X r

to be rank deficient, as indicated in Eq. (6)), so that

CV — 0 and rank

(n+p)xn

Premultiplying the first equation of Eq. (11) by V eliminates the Lagrange multipliers (it is very important that F does not depend on x, hence X does not appear on the right-hand side of Eq. (11) and can be eliminated) and yields the following set of equations (sometimes called a null space formulation [18, 29]):

Matrix M is invertible, thus—due to Eq. (12)—the rank of the leading matrix in Eq. (13) equals n. Constraints (8) are consistent; hence Eq. (13) has a unique solution for q. Obviously, it is possible to find another matrix V* = V, spanning the null space of the constraint matrix C, and to obtain a different version of Eq. (13). In such a case, the following holds:

"Vr M" "Vr F"

C q = r

V* = VA,

(n-r)x(n-r),

where A is an invertible transformation matrix. Thus, when V* appears instead of V in Eq. (13), it is sufficient to premultiply the first equation by A-r, to obtain the original form of Eq. (13). This indicates that the solution for q does not depend on the choice of matrix V, thus it is a unique solution.

It should be noted that, since only constraints at acceleration level are represented in Eq. (11), some stabilization techniques are usually adopted to account for constraints given by Eqs. (1) and (4).

3 Constraint reactions

3.1 Methods of handling redundant constraints

The presence of redundant constraints induces rank deficiency of constraint matrix C, thus special methods are needed to handle equations of motion. It is possible to use algorithms capable to deal with dependent equations, e.g., the pseudoinverse-based methods [32] or the augmented Lagrangian method [2, 3]. The other possibility is to modify the set of equations in order to exclude dependent equations [16, 28]; and this approach—usually referred to as redundant constraints elimination method—was adopted here to simulate the exemplary mechanism, presented in the subsequent section. It should be emphasized that, regardless of the method used to handle redundant constraints, the problems with Lagrange multipliers uniqueness are circumvented, not solved [36].

The redundant constraints elimination method consists in removing dependent equations from the mathematical description of multibody system, hence only the subset of independent equations is analyzed. The redundant constraints elimination may be done "manually"—selected constraints, which are unnecessary when kinematics is concerned, are simply not included into the model [28]. The other possibility is to eliminate redundant constraints automatically. The automatic elimination may be performed prior to constraint equations formulation [23] or during the preprocessing phase of simulation [25]. However, the most popular and the simplest method of automatic redundant constraints elimination (implemented in widely used multibody software) is based on the constraint matrix C investigation [16, 28]. In this method, all constraint equations are formulated and then they are divided into subsets of independent and dependent equations. The division is usually made according to the results of the constraint matrix decomposition. The equations classified as dependent are excluded from the mathematical model.

It must be stressed that, regardless of redundant constraints selection method, the reactions associated with eliminated constraints are set to zero. It is quite obvious that in real mechanisms reactions corresponding to constraints neglected during analysis are not likely to be constantly equal zero. Moreover, it is worth noting that setting the reactions of eliminated constraints to zero, transfers their loads to the constraints that remain in the mathematical model. Consequently, redundant constraints elimination affects not only reactions of excluded constraints but also reactions of remaining constraints.

3.2 Uniquely determined constraint reactions

The generalized constraint reactions, for the whole multibody system, can be calculated as [10, 16, 28]:

f = CTX = F - Mq. (15)

Similarly, the generalized reactions in a specified kinematic pair can be calculated as

fx = CTX Xx, (16)

where Cx denotes a matrix built of these rows of C that correspond to constraint equations representing the kinematic pair, and Xx is a vector of Lagrange multipliers associated with these constraints. Note that Eq. (4) can always be reordered, so that

where none of the rows of matrix CY corresponds to the considered kinematic pair.

In the case of redundant constraints presence, the generalized reactions for the whole system (f) are uniquely determined, however, some or all reactions between pairs of interacting bodies (fx) cannot be uniquely determined using a rigid body model. In other words, it is impossible to tell how the total reaction f is distributed between the individual constraint reactions fx. This is an immediate consequence of the nonuniqueness of multipliers x, discussed in Sect. 2.2.

Fortunately, in many cases, some of the constraint reactions are solvable (i.e., they can be uniquely determined using a rigid body model) despite the presence of redundant constraints in the system taken as a whole. It is possible to generalize the methods of reactions solvability analysis, presented in [34] (for holonomic systems only), to the case of systems subjected to both holonomic and linear nonholonomic constraints.

The concept of direct sum, known from linear algebra [31], can be exploited to check whether the particular examined reactions can be uniquely calculated. To recall the definition of direct sum, assume that Z is a linear vector space in Rn whereas X and Y are the subspaces of Z. We say that Z is a direct sum of subspaces X and Y, and we denote it as Z = X 0 Y, when the following conditions are fulfilled:

1. Z is a sum of subspaces X and Y (we denote it as Z = X + Y), which means that any

vector z € Z can be represented as z = x + y, where x e X and y e Y.

2. If xi + y1 = x2 + y2, provided that xi, x2 e X, and y:, y2 e Y, then xi = x2 and = y2.

To formulate analytical condition for checking constraint reactions solvability, let us consider a kinematic pair represented by holonomic or nonholonomic constraint equations. Let X be a linear space spanned by these columns of matrix CT that correspond to the constraints representing this kinematic pair (X = span(CX)), and let Y be a linear space spanned by the

other columns of matrix CT (Y = span(CY)). Moreover, let Z be a linear space spanned by all columns of matrix CT (Z = span(CT)). Note that in the case of constraint redundancy, the following inequality is valid (compare Eq. (6)):

dim(Z) = rank(C) = r<m + h, (18)

thus the same vector of generalized reactions f can be obtained for different vectors of Lagrange multipliers X (see Eq. (15)). It can be shown that:

If Z is a direct sum of X and Y (Z = X 0 Y), then the generalized constraint reaction fX in the considered kinematic pair is uniquely determined.

To prove the above statement, the properties of direct sum can be utilized. Let us write the Lagrange multipliers vector as X = [ XTX XY ]T, where XX corresponds to the vectors spanning the linear space X and xy corresponds to the vectors spanning the linear space Y. Assume that x* and x** are two different vectors of Lagrange multipliers. The generalized reactions corresponding to these multipliers are the following:

f* = ct x* = cX xX + cY xY = fX + fY (19)

f** = ct x** = cX x*X* + cT xy* = fX* + fY*- (20)

Since f* = f**, the 2nd condition of the direct sum definition yields fX = fX* (and similarly,

fY = fY*).

The presented considerations show that, when the direct sum condition is fulfilled, the term fX = CXXX (generalized reaction force in the considered kinematic pair—see Eq. (16)) can be uniquely determined. It should be stressed, however, that the problem of finding Lagrange multipliers XX may not have a unique solution. The Lagrange multipliers XX can be uniquely determined only if CX is a full rank matrix.

The above results can be utilized to check whether constraint reactions in the selected kinematic pair can be uniquely determined, regardless of possibility that the problem of calculating Lagrange multipliers may have infinitely many solutions.

3.3 Numerical detection of uniquely determined constraint reactions

Three numerical methods of examination whether the direct sum condition is fulfilled, and thus whether reactions in the particular kinematic pair are uniquely solvable, were proposed in [34]. The first method is based on ranks comparison of matrices CX, CY, and C, the second is based on singular value decomposition of matrix C; these methods will not be discussed here. This paper presents an improved version of the third method (the most effective one), which is based on QR decomposition of the constraint matrix C. The QR method decomposes matrix C into matrices [12, 17]:

CE = QR or CT = ERT QT, (21)

where Q is an orthogonal ((m + h) x (m + h)) matrix, E is an orthogonal (n x n) permutation matrix (it consists of zeros and ones only) and R is a rectangular ((m + h) x n) upper triangular matrix.

In the literature, one can easily find a QR decomposition without the use of permutation matrix E, however, this variant of decomposition is not convenient for us. The column permutation E is chosen in such a way that absolute values of diagonal elements in R are decreasing. This fact is important for our considerations, since it assures that, for decomposed matrix C with rank r, the last p (see Eq. (7)) rows of matrix R consist of zeros only.

The QR decomposition of constraint matrix C is useful when formulating a numerical method of checking whether an arbitrarily chosen vector w e Z can be decomposed into sum of two vectors wX e X and wY e Y, both of them given uniquely. The proposed criterion is based on the 2nd condition of direct sum of spaces definition.

Let us consider vector w e Z and write it as a linear combination of matrix CT columns:

w = WX + WY = CX nX + CY nY = CT n = ERT QT n, (22)

where n(m+h)x1 = [nX nY ]T is a vector of linear combination coefficients.

An auxiliary vector d, defined below, is more suitable for further considerations:

d(m+h)x1 = QT n. (23)

The transformation (23) maps d to n and n to d uniquely, since Q is a square and orthogonal matrix. This transformation can be used to rewrite Eq. (22) in the following form:

w = CT n = ERT d. (24)

We can state that Eq. (24) uniquely defines only the first r elements of vector d, since only the first r rows of matrix R are nonzero. The remaining p elements of d can be arbitrarily chosen. Vector d satisfying Eq. (24), and thus vector n, can be written in the following form:

d — d + d, d— [di

d — [o •

n — Qd —

dm+h ] ,

nx Qx d— Qx

.nr. .Qr. Qr

(d + d),

where elements d1,...,dr are uniquely defined by (24) and elements dr+1 ,...,dm+h can be arbitrarily chosen; QX denotes a matrix consisting of only these rows of Q that correspond to the considered kinematic pair.

Hence, vector wX can be written in the form:

wx — nx — Qx d + Qx d.

The first r elements of vector d are zeros, the remaining p can be chosen arbitrarily. Thus, we can state that vector wX can be determined uniquely only when the arbitrary elements of d vanish due to multiplication by zero. This observation implies the following definition of an auxiliary matrix BX:

[Bx ]nxp — CX Qx>1'

where matrix QX consists of last p columns of QX.

In conclusion, we can say that vector wX can be determined uniquely when matrix BX consists of zeros only. This result can be easily used for detecting kinematic pairs for which

reactions can be uniquely determined. Firstly, the QR decomposition of the constraints matrix C should be calculated. Then, rows of matrices Q and C, corresponding to the kinematic pair being investigated, ought to be extracted. Finally, matrix Bx should be calculated and examined if it contains only zero elements.

It should be mentioned that the method presented here, as well as two other methods described in [34], merely require operations on the constraint matrix C, therefore, it is relatively easy to implement them in a multibody code.

4 Example

4.1 Structure and kinematics of the mechanism

Since applications of multibody methods to robotics have always been interesting to us [7, 9, 22], a mobile robot was a natural candidate to serve as an example in this paper. The constraint reaction solvability method was utilized to analyze a simplified model of a four-wheel mobile robot connected with a one-wheel trolley (Fig. 1). The system is modeled as a planar mechanism, in which wheels are replaced by semicircular knife edges. Edges Wi and W2 can change their inclinations with respect to the robot platform 1, whereas edges W3 and W4 are rigidly attached to the platform. The linkage consisting of bodies 2, 3, 4, 5, and 6, driven in joint H, is used to rotate edges Wi and W2 in a synchronized manner. Due

Fig. 1 Simple model of the mobile robot and its kinematic scheme

8 0.02

8 0 CL

Q -0.02

- dw \

Time (s)

Fig. 2 Relative displacement in the translational joint vs. time

to appropriately chosen linkage dimensions, it is possible to change the wheel inclination by ±90°. The edge W5 is rigidly attached to the trolley platform 7. The trolley is connected to the robot via revolute joint A.

A global inertial reference frame n0 = x0y0 is established on the ground 0, and body-fixed local reference frames n = x,yi are embedded in the moving bodies (Fig. 1). The absolute Cartesian coordinates that describe the mechanism form vector q:

q = [qf qi qi q! qi qff, qi = [rf wf,

where r, = [xt y, ]T represents the position of the local reference frame n, origin with respect to the global frame n0, and ^ is the angle of the local frame n, rotation with respect to the global frame n0.

Seven revolute joints (A, B, C, D, E, F, G) and one translational joint (H) are described by holonomic (geometric) constraints. Appropriate equations are presented in the Appendix. Driving constraints (represented by holonomic and rheonomic constraints, described in the Appendix) and are imposed on relative motion in joint H. The time history of the displacement d(t) of body 6 with respect to body 1 is presented in Fig. 2. The steering is chosen so that the robot, which initially is moving along x0 axis, firstly turns right, then travels forward, and finally turns left and starts to move along a line parallel to x0.

The holonomic constraints imposed on the modeled multibody system can be written jointly as a set of m = 17 nonlinear algebraic equations:

* (q,t) = [

= i * A

*h =01

Formulas for the constraint Jacobian matrix $q and for vector $t (see Eq. (2)) as well as for r vector (see Eq. (9)) are presented in the Appendix.

The ground-wheel interactions are represented by knife-edge Pfaffian constraints. The knife-edge constraint equation can be derived by requiring that linear velocity at the point of contact is parallel to the edge. The following scalar equation, linear in velocities, is obtained:

K o SK

(0)'-(0) = (R. nff (r, + o R sK^O = [ (R, nf)T n

(0) 1 O tll^i 1l Vl^iO-f /^l/"»! J n t />i"\nto /■»+• r\A1nt Of Ofli^ (0)

= vKq, = 0, (30)

where vK) is the linear velocity at contact point K, and nK) is a unit vector perpendicular to the allowed direction of motion (the components of both vectors are expressed in the

Table 1 Geometric parameters of knife-edge constraints

Edge Body i Point K s(i) sK (i) nK

W1 2 B [0 of [0 1]T

W2 3 C [0 of [0 1]T

W3 1 L [0 - 4a]T [0 1]T

W4 1 M [0 4a]T [0 1]T

W5 7 N [0 0]T [0 1]T

a = 0.05 (m)

global frame n0). The other symbols (i, i2, R;, sK) are explained in the Appendix (note that Rf i2 R; = i2). The components of vectors are presented in Table 1 (dimensions of the robot are presented in the Appendix).

In the case of knife-edge constraint defined by Eq. (30), r is a scalar given by the following formula:

rK = -(QR;4>y ii<Pi. 4.2 Constraint reactions solvability analysis

Detailed formulas for constraint Jacobian entries are provided in the Appendix. There are m = 17 holonomic constraint equations and n = 21 coordinates, thus the Jacobian matrix

0 q has 17 rows and 21 columns.

The nonholonomic constraints are present due to h = 5 knife-edge kinematic pairs, thus Eq. (30) can be used to calculate the nonzero entries of the Pfaffian constraints matrix for the whole system:

0 9 B 0 0 0 0 0

0 0 9 C 0 0 0 0

9 L 0 0 0 0 0 0

9 M 0 0 0 0 0 0

0 0 0 0 0 0 9 N

Appending the Pfaffian constraints matrix to the Jacobian matrix gives the constraint matrix C:

Let q0 denote the vector of coordinates representing the mechanism configuration shown in Fig. 1:

q0 = [20a 0 0 | 24a 0 0 | 16a 0 0 •••

20a + dx dy 0 | 20a - dx dy 0 j 20a 0 0 j 5a 0 0]T,

where a = 0.05 (m), dx = cx/2, dy = cy/2 — by, and quantities cx, cy, by are defined in the Appendix.

It can be calculated that

rank($q(q0)) = 17, rank(9(q0)) = 4, r = rank(C(q0)) = 20 < 17 + 4. (35)

Matrix $ q is a full-rank matrix, whereas 9 is rank-deficient. Obviously, the constraint matrix C is rank-deficient as well. The same results can be obtained for any other nonsingular configuration q of the mechanism. This indicates that the considered multibody system is overconstrained. The degree of redundancy (see Eq. (7)) can be calculated as

p = m + h - r = 17 + 5 - 20 = 2.

Dependency of constraints was detected, thus the constraint reaction solvability analysis was performed prior to dynamic simulations. Firstly, QR decomposition of the constraint matrix C was calculated (see Eq. (21)). Then matrix C was divided into fourteen submatrices, corresponding to seven revolute joints, one translational joint, one driving constraint equation, and five knife-edge constraints, respectively. Similarly, matrix Q was divided into fourteen submatrices:

[CrevA]2x21

[CrevG]2x21 [Ctraff ]2x21 [Cdrvng]1x21 [Ced^1 ]1x18

. [CedW5]1x21

[QrevA]2x22

[QrevG]2x22 [Qtraff ]2x22 [Qdrvng]1x22 [QedW1 ]1x22

. [QedW5]1x22 .

Next, Eq. (27) was utilized to calculate eleven matrices BX (note that degree of redundancy p = 2, thus only two last columns of QX matrices were used to build corresponding QX matrices):

0 0.408 0.082 0 -0.408 01x16 0 0 0 0 0 01x16

0 0.408 -0.082 01x4 -0.408 01x13

01x4 0.408 01x4 0

01x17 01x17

01 x 4 T

edW2 ■

01x7 01x7

0.408 01x14 0 01x14

0 -0.408 01x20 0 0.707 01x20

0 -0.408 01x20 0 -0.707 01x20

Brev^ = BrevD = BrevE = BrevF = BrevG = Btraff = Bdrvng = BedW< = 0:

Matrices corresponding to knife-edge kinematic pairs W1-W4 as well as matrices corresponding to revolute joints B and C have nonzero elements, thus constraint reactions in these pairs cannot be uniquely determined using the rigid body model. Matrices corresponding to revolute joints A, D, E, F, and G, to translational joint H, and to the driving constraint equation consist of zeros only, thus related reactions can be uniquely determined. This concludes the constraint reactions solvability analysis.

Table 2 Masses and central moments of inertia

Body i 1 2 3 4 5 6 7

mi (kg) Ji (kg m2) 10 0.25 1 0.005 1 0.005 1 0.005 1 0.005 1 0.005 5 0.05

4.3 Simulated motion of the mechanism

Masses of mechanism bodies and their central moments of inertia are presented in Table 2. In each moving part, the center of its mass coincides with the local reference frame origin.

Equations of motion for the investigated mechanism can be written in the form of Eq. (11), with no external forces applied to the system (note that the edge-ground forces are modeled as constraint reactions), with M21x21 = diag(m1,m1,/1 ,...,m7,m7,J7), and with terms C and r discussed previously.

The mechanism is overconstrained, therefore, the leading matrix in the equations of motion is singular. The constraints elimination method (see Sect. 3.1) was chosen to handle the redundancy problem. Redundant constraints may be selected in many ways, thus several variants of dynamic simulations were carried out. In each simulation, the equations of motion were integrated (the fourth-fifth order Runge-Kutta formula [6], implemented in MATLAB® as ode45 function, was used) with the initial configuration given by q0 (see Eq. (34)) and initial velocity q0:

q0 = [ v o o ! v o o ! v o o | v 0 0 •••

v 0 0 | v 0 0 | v 0 0]T, v = 0.1 (m/s). (39)

Diverse joint reactions were calculated for different selections of redundant constraints, however, in each case the same mechanism motion (to within numerical precision) was obtained. This illustrates that—as it was discussed in Sect. 2.2—in the case of an over-constrained mechanism with frictionless kinematic pairs, the individual reactions are non-unique but their resultant effect (when motion is concerned) is unique.

The trajectory of robot platform center of mass as well as the trajectories of knife-edge contact points, observed during simulations, are presented in Fig. 3.

4.4 Calculated constraint reactions

Three variants (out of several other possible) of redundant equations elimination were studied and presented here. Since the degree of redundancy equals two, in each possible variant two constraint equations need to be eliminated.

In the first variant, scalar constraint equations no. 6 and 21 were eliminated from the set of constraints at velocity level, i.e., from Eq. (4). Appropriate elimination of matrix C (see Eq. (33)) rows no. 6 and 21, as well as vector r elements no. 6 and 21, was performed. Consequently, the Lagrange multipliers associated with the neglected equations were canceled, and a reduced set of equations of motion (Eq. (11)), with nonsingular leading matrix, was obtained. Note that the eliminated constraint equation no. 6 represents revolute joint C, thus—as a result of the elimination process—the y component of joint C reaction was arbitrarily set to zero. Constraint equation no. 21 represents knife-edge kinematic pair W4, thus reaction in this pair was set to zero as well.

■0.2

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Position x (m)

Fig. 3 Trajectories of selected points of the robot

In the second variant of redundant equations elimination, constraints no. 18 and 21, representing knife-edges W1 and W4, respectively, were eliminated. In the third variant, constraint equations no. 19 and 21, associated with knife-edges W2 and W4, respectively, were deleted. Again, appropriate reactions were set to zero due to purely mathematical operations, not supported by the physics of the system.

For all three variants of redundant constraints elimination, dynamic simulations were carried out and constraint reactions were observed (as it was mentioned earlier, simulated motion was the same in all cases). Reactions observed in joints A, D, E, F, G, H, the steering force associated with driving constraints, and knife-edge W5 reaction had the same time histories in all simulations, whereas reactions in joints B and C, as well as knife-edges W1-W4 reactions, were different in each case. The results of simulations corroborate correctness of constraint reactions solvability analysis performed in Sect. 4.2.

Global x and y components of revolute joint B reaction force and perpendicular to the edge component of the knife-edge W3 reaction are presented in Fig. 4 as examples of nonunique reactions. In the same Fig. 4, components of revolute joint D reaction force and the knife-edge W5 reaction are shown as examples of uniquely solvable reactions.

The alternatively eliminated constraint equations are associated with kinematic pairs C, W1, W2, and W4. It is worth noting that in each case elimination affected, among others, reactions in pairs B and W3. This result shows that the process of elimination influences not only reactions of neglected constraints but also reactions of constraints remaining in the mathematical model.

It should be mentioned that the redundant constraints elimination method can be seen as a substitution of the original mechanism by a kinematically equivalent mechanism without redundant constraints. The problem is that—when reactions are concerned—there is no equivalent model without redundant constraints. Elimination of a constraint results in setting its reaction to zero, which is unlikely to be observed in the original mechanism.

5 Conclusions and comments

The presence of redundant constraints leads to nonuniqueness of some or all reactions, as long as the mechanism parts are treated as rigid bodies. Mathematically, redundant con-

Fig. 4 Nonunique (left) and unique (right) constraint reactions

straints existence results in system of equations with infinite number of solutions. Thus, the rigid body equations of motion are insufficient to calculate realistic joint reactions, observed in real-world overconstrained multibody systems.

Although the problem of reactions uniqueness—crucial for results credibility—is often encountered in practical calculations, surprisingly small attention is paid to it. In the great majority of general purpose multibody simulation packages detailed information whether particular joint reaction is solvable or unsolvable is unavailable, and the problem of reactions solvability is simply neglected. The software users are frequently advised to replace over-constrained models with kinematically equivalent models without redundant constraints. The problem that the nonredundant model does not reflect the real system is usually not discussed. If the software user fails to follow the advice, the redundant constraints detected in the model are handled automatically, using one of several available methods. This procedure may easily lead to erroneous results. For example, quite often the multibody software user is not aware that specified joint reaction force is nonunique and his attempt to introduce friction into this joint is pointless.

This paper presents a theoretical background and numerical methods for detecting kinematic pairs of rigid body mechanism for which reactions can be uniquely determined deSpringer

spite the redundancy of constraints. The reactions solvability analysis methods, previously proposed for holonomic systems, are generalized to mechanisms subjected to linear non-holonomic constraints. It should be emphasized that the obtained results are valid for both integrable and nonintegrable Pfaffian constraints. An example of a mobile robot analysis is provided to illustrate the discussed problems.

The proposed reactions solvability methods are based on constraint matrix investigation, thus they can be easily implemented in existing multibody software. Required computational effort is small, since the constraint matrix must be computed anyway. Moreover, in most cases, reactions solvability analysis can be performed only once, during the preprocessing phase of simulation (provided that the initial mechanism configuration is nonsingular), along with redundant constraint elimination process.

The uniqueness of reactions depends only upon the kinematic structure of the mechanism, thus it is irrelevant which type of coordinates is utilized to model the multibody system. The methods presented in the paper employ absolute coordinates' formulation. Adopting them to use other types of coordinates, for example, relative ones, is not straightforward and requires additional research.

Finally, going beyond the scope of this paper, we should point out that in order to find unique values of all reactions, rigid body assumption must be rejected. Engineering experiences suggest that flexibility of bodies, elasticity of bearings, assembly stresses, thermal loads, joint clearances as well as geometric imperfections (in [24], it was shown that, in the case of redundantly constrained mechanisms, the rigid body assumption must be accompanied by a postulate that links are geometrically perfect) should be taken into account, when solving for realistic reactions. In many models, elasticity of bodies is considered, however, some effects associated with constraint redundancy are still observed, even when mechanisms are modeled as entirely flexible systems, e.g., [1, 21]. Quite often, in a multibody model, only selected rigid bodies are replaced by their deformable substitutes, e.g., [37], which not always leads to correct results. It can be shown that problems with finding realistic reactions, typical for overconstrained systems, may be observed when analyzing partially flexible models without redundant constraints [35]. Hence, problems with finding reactions in overconstrained mechanisms are not limited to purely rigid body models.

Acknowledgements Research supported by the Institute of Aeronautics and Applied Mechanics statutory funds.

Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

Appendix: Holonomic constraints and their derivatives

The investigated mechanism, presented in Fig. 1, consists of seven moving bodies interconnected by seven revolute joints and one translational joint. The joints are represented by geometric constraint equations. Moreover, driving constraints are imposed on relative motion in the translational joint.

Let us consider a revolute joint formed by bodies i and j at point Pi = Qj (for example, in Fig. 1 joint D is formed by bodies 2 and 4, P2 = Q4 = D). The constraint equations describing this joint can be derived by requiring that the point Pi on body i coincides with point Qj on body j :

*r = ri + Risf - rj - RjsQ = 02xi, (40)

Table 3 Geometric parameters of joints

i j sP s(j) sQ vij fij

A 7 1 [7 a 0]T [-8a 0]T — -

B 2 1 [0 0]T [4a 0]T - -

C 3 1 [0 0]T [-4a 0]T - -

D 4 2 [Cx/2 - Cy/2]T [-bx - by]T - -

E 5 3 [-Cx/2 - Cy/2]T [bx - by]T - -

F 6 4 [0 Cy - by]T [-Cx/2 Cy/2]T - -

G 6 5 [0 Cy - by]T [Cx/2 Cy/2]T - -

H 6 1 [0 10a]T [0 - 10a]T [10]T 0

a = 0.05 (m), ç = arctan( 1/4), bx = a sin ç, by = a cos ç, cx = 4a - bx, Cy =y J(5a)2 - C2

where s® is the position vector of point Uk in the local reference frame nk and Rk is the direction cosine matrix transforming quantities from nk to n0:

Rk = Rk(çk) =

cos Çk sin Çk

- sin Çk cos Çk

Let us consider a translational joint, formed by bodies i and j, and assume that point Pi belongs to body i, and point Qj as well as vector vij belong to body j. Points Pi and Qj lie on a line parallel to the axis of relative joint motion, and vector vij is perpendicular to this axis. The translational joint can be described by two scalar constraint equations. The first one represents the fact that vector Pi Qj is perpendicular to vector vij, and the second equation represents the fact that body i does not change its orientation with respect to body j :

(rj + RjsQ - r - R;s«)TRjvj

Ç; - Çj - fij

where ^¡j is a constant angular value.

The driving constraint (single scalar equation) for a translational joint can be derived by requiring that the distance between points p and Qj, at every time instant, is given by a scalar, time-dependent function d(t):

= (rj + RjsQ - r; - R;sff Rjuj) - d(t) = 0

where u,j is a unit vector parallel to the line of relative motion.

Geometric parameters characterizing joints are presented in Table 3. Note that a = 0.05 (m), and BD = CE = a, DF = EG = 5a, BC = LM = 8a, AC = 4a, AN = 7a (see Fig. 1).

The nonzero Jacobian matrix entries (see Eq. (2)) corresponding to Eq. (40) can be calculated as

= k = [I2X2 ßR;s«], = [ * rj ] = [-I2X2 -a Rj sQ ], a =

0 -1 1 0

The nonzero Jacobian matrix entries corresponding to Eq. (42) are the following:

$q — r*r _t-(r,vj))T -(RvjVQR,s

vq' L r' ni L 0ix2

]—[(R0 v

(Rjv(f)T -(Rjv(f)TQ(rj - r, - R,s«)'

j — 1

The nonzero Jacobian matrix entries corresponding to Eq. (43), and the nonzero elements of the driving constraints partial derivative with respect to time (see Eq. (2)), are the following:

< < — <] — [-(Rj -(Rj ujy QR, s«]

) ]■

< _ <] _ [(Rjuj))T - (Rj u<(Jj))T Q (rj - r, - R,s

Vectors r (see Eq. (9)), representing revolute and translational joints as well as the driving constraints, can be respectively expressed as

rr — R, spV,2 - Rj sftf,

(RjviJj )T(2Q(rj - ri),pj + (rj - r,)^2 - R,sf^j - ¿)2)'

rd — (Rjuj))T(2q(rj - r,)<pj + (rj - r,)^2 - R,sf&j - <p,)2) + dtt.

References

1. Aarts, R.G.K.M., Meijaard, J.P., Jonker, J.B.: Flexible multibody modelling for exact constraint design of compliant mechanisms. Multibody Syst. Dyn. 27(1), 119-133 (2012)

2. Bayo, E., Ledesma, R.: Augmented Lagrangian and mass-orthogonal projection methods for constrained multibody dynamics. Nonlinear Dyn. 9(1-2), 113-130 (1996)

3. Bayo, E., Garcia de Jalon, J., Serna, M.A.: A modified Lagrangian formulation for the dynamic analysis of constrained mechanical systems. Comput. Methods Appl. Mech. Eng. 71(2), 183-195 (1988)

4. Blajer, W.: On the determination of joint reactions in multibody mechanisms. J. Mech. Des. 126(2), 341-350 (2004)

5. Bloch, A.: Nonholonomic Mechanics and Control. Interdisciplinary Applied Mathematics, vol. 24. Springer, New York (2003)

6. Dormand, J.R., Prince, P. J.: A family of embedded Runge-Kutta formulae. J. Comput. Appl. Math. 6(1), 19-26 (1980)

7. Fraczek, J., Wojtyra, M.: Teaching multibody dynamics at Warsaw University of Technology. Multibody Syst. Dyn. 13(3), 353-361 (2005)

8. Fraczek, J., Wojtyra, M.: On the unique solvability of a direct dynamics problem for mechanisms with redundant constraints and Coulomb friction in joints. Mech. Mach. Theory 46(3), 312-334 (2011)

9. Fraczek, J., Wojtyra, M.: Application of general multibody methods to robotics. In: Arczewski, K., et al. (eds.) Multibody Dynamics—Computational Methods and Applications, pp. 131-152. Springer, Berlin (2011)

10. Garcia de Jalon, J., Bayo, E.: Kinematic and Dynamic Simulation of Multibody Systems: The Real-Time Challenge. Springer, New York (1994)

11. Garcia de Jalon, J., Gutierrez-Lopez, M.D.: Multibody dynamics with redundant constraints and singular mass matrix. In: Proc. of the 2nd Joint International Conference on Multibody System Dynamics, Stuttgart, Germany (2012)

12. Golub, G.H., Van Loan, C.F.: Matrix Computations, 3rd edn. Johns Hopkins University Press, Baltimore (1996)

13. González, F., Kövecses, J.: Use of penalty formulations in the dynamic simulation of redundantly constrained multibody systems. In: Proc. of the 2nd Joint International Conference on Multibody System Dynamics, Stuttgart, Germany (2012)

14. González, F., Kövecses, J.: Use of penalty formulations in dynamic simulation and analysis of redundantly constrained multibody systems. Multibody Syst. Dyn. 29(1), 57-76 (2013)

15. González, F., Kövecses, J., Teichmann, M., Courchesne, M.: Techniques for the simulation of redundantly constrained multibody systems. In: Proc. of the ECCOMAS Thematic Conference Multibody Dynamics 2011, Brussels, Belgium (2011)

16. Haug, E.J.: Computer Aided Kinematics and Dynamics of Mechanical Systems. Allyn & Bacon, Boston (1989)

17. Kim, S.S., Vanderploeg, M.J.: QR decomposition for state space representation of constrained mechanical dynamic systems. J. Mech. Transm. Autom. Des. 108, 176-182 (1986)

18. Laulusa, A., Bauchau, O.A.: Review of classical approaches for constraint enforcement in multibody systems. J. Comput. Nonlinear Dyn. 3(1), 011005 (2008)

19. Malczyk, P., Fraczek, J.: Cluster computing of mechanisms dynamics using recursive formulation. Multi-body Syst. Dyn. 20(2), 177-196 (2008)

20. Malczyk, P., Fraczek, J.: A divide and conquer algorithm for constrained multibody system dynamics based on augmented Lagrangian method with projections-based error correction. Nonlinear Dyn. 70(1), 871-889 (2012)

21. Meijaard, J.P., Brouwer, D.M., Jonker, J.B.: Analytical and experimental investigation of a parallel leaf spring guidance. Multibody Syst. Dyn. 23(1), 77-97 (2010)

22. Mianowski, K., Wojtyra, M.: Virtual prototype of a 6-DOF parallel robot. In: Proc. of Eleventh World Congress in Mechanism and Machine Science, Tianjin, P.R. China, pp. 1604-1608 (2004)

23. Müller, A.: A conservative elimination procedure for permanently redundant closure constraints in MBS-models with relative coordinates. Multibody Syst. Dyn. 16(4), 309-330 (2006)

24. Müller, A.: Generic mobility of rigid body mechanisms. Mech. Mach. Theory 44(6), 1240-1255 (2009)

25. Müller, A.: Semialgebraic regularization of kinematic loop constraints in multibody system models. J. Comput. Nonlinear Dyn. 6(4), 041010 (2011)

26. Neimark, J.I., Fufaev, N.A.: Dynamics of Nonholonomic Systems. Am. Math. Soc., Providence (1972)

27. Neto, M.A., Ambrosio, J.: Stabilization methods for the integration of DAE in the presence of redundant constraints. Multibody Syst. Dyn. 10(1), 81-105 (2003)

28. Nikravesh, P.E.: Computer-Aided Analysis of Mechanical Systems. Prentice Hall, Englewood Cliffs (1988)

29. Phong, D.V., Hoang, N.Q.: Singularity-free simulation of closed loop multibody systems by using null space of Jacobian matrix. Multibody Syst. Dyn. 27(4), 487-503 (2012)

30. Ruzzeh, B., Kövecses, J.: A penalty formulation for dynamics analysis of redundant mechanical systems. J. Comput. Nonlinear Dyn. 6(2), 021008 (2011)

31. Strang, G.: Linear Algebra and Its Applications, 4th edn. Brooks/Cole, Pacific Grove (2005)

32. Udwadia, F.E., Kalaba, R.E.: Analytical Dynamics: A New Approach. Cambridge University Press, New York (1996)

33. Wojtyra, M.: Joint reaction forces in multibody systems with redundant constraints. Multibody Syst. Dyn. 14(1), 23-46 (2005)

34. Wojtyra, M.: Joint reactions in rigid body mechanisms with dependent constraints. Mech. Mach. Theory 44(12), 2265-2278 (2009)

35. Wojtyra, M., Fraczek, J.: Joint reactions in rigid or flexible body mechanisms with redundant constraints. Bull. Pol. Acad. Sci., Tech. Sci. 60(3), 617-626 (2012)

36. Wojtyra, M., Fraczek, J.: Comparison of selected methods of handling redundant constraints in multi-body systems simulations. J. Comput. Nonlinear Dyn. 8(2), 021007 (2013)

37. Zahariev, E., Cuadrado, J.: Dynamics of over-constrained rigid and flexible multibody systems. In: Proc. of the 12th IFToMM. World Congress, Besançon, France (2007)