Solvability of reactions in rigid multibody systems with redundant nonholonomic constraints

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.


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 frictionless 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.
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 holonomic 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.

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: 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: 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 velocitiesq [26]: 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: where respective definitions of coefficient matrix C (often called a constraint matrix) and vector b are the following: We assume that constraints (4) are consistent, i.e., no contradictory conditions are imposed on velocitiesq. For our further considerations, it is important that when the mechanism is overconstrained, rows of matrix C are linearly dependent, hence It is worth noting that constraint matrix C may be rank-deficient even when both Φ q and Ψ are full-rank matrices.
The degree of redundancy p is defined as The velocity constraints (4) can be differentiated with respect to time in order to obtain constraints for accelerations: with vector Γ defined as 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.

Equations of motion
Equations of motion for the considered multibody system can be written in the following form [10,16]: 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 λ 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 differentialalgebraic equations: Equation (11) can be solved for the accelerations and the Lagrange multipliers (it can be shown that the Lagrange multipliers λ are uniquely solvable only when matrix CM −1 C T 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 accelerationsq is unique, and the solution for multipliers λ is non-unique. To justify this statement, it is helpful to define matrix V n×(n−r) , which is an orthogonal complement to C (which is now assumed to be rank deficient, as indicated in Eq. (6)), so that Premultiplying the first equation of Eq. (11) by V T eliminates the Lagrange multipliers (it is very important that F does not depend on λ, hence λ 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 forq. 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: 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 −T , to obtain the original form of Eq. (13). This indicates that the solution forq 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).

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.

Uniquely determined constraint reactions
The generalized constraint reactions, for the whole multibody system, can be calculated as [10,16,28]: Similarly, the generalized reactions in a specified kinematic pair can be calculated as where C X denotes a matrix built of these rows of C that correspond to constraint equations representing the kinematic pair, and λ X 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 C Y 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 (f X ) 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 f X . This is an immediate consequence of the nonuniqueness of multipliers λ, 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 R n 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 ⊕ 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 ∈ X and y ∈ Y . 2. If x 1 + y 1 = x 2 + y 2 , provided that x 1 , x 2 ∈ X, and y 1 , y 2 ∈ Y , then x 1 = x 2 and y 1 = y 2 .
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 C T that correspond to the constraints representing this kinematic pair (X = span(C T X )), and let Y be a linear space spanned by the other columns of matrix C T (Y = span(C T Y )). Moreover, let Z be a linear space spanned by all columns of matrix C T (Z = span(C T )). Note that in the case of constraint redundancy, the following inequality is valid (compare Eq. (6)): thus the same vector of generalized reactions f can be obtained for different vectors of Lagrange multipliers λ (see Eq. (15)). It can be shown that: If Z is a direct sum of X and Y (Z = X ⊕ Y ), then the generalized constraint reaction f X 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 where λ X corresponds to the vectors spanning the linear space X and λ Y corresponds to the vectors spanning the linear space Y . Assume that λ * and λ * * are two different vectors of Lagrange multipliers. The generalized reactions corresponding to these multipliers are the following: and Since f * = f * * , the 2nd condition of the direct sum definition yields f * X = f * * X (and similarly, . The presented considerations show that, when the direct sum condition is fulfilled, the term f X = C T X λ X (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 λ X may not have a unique solution. The Lagrange multipliers λ X can be uniquely determined only if C X 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.

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 C X , C Y , 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]: where Q is an orthogonal ((m+h)×(m+h)) matrix, E is an orthogonal (n×n) permutation matrix (it consists of zeros and ones only) and R is a rectangular ((m + h) × 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 ∈ Z can be decomposed into sum of two vectors w X ∈ X and w Y ∈ 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 ∈ Z and write it as a linear combination of matrix C T columns: where n (m+h)×1 = [ n T X n T Y ] T is a vector of linear combination coefficients. An auxiliary vector d, defined below, is more suitable for further considerations: 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: 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: where elements d 1 , . . . , d r are uniquely defined by (24) and elements d r+1 , . . . , d m+h can be arbitrarily chosen; Q X denotes a matrix consisting of only these rows of Q that correspond to the considered kinematic pair. Hence, vector w X can be written in the form: The first r elements of vectord are zeros, the remaining p can be chosen arbitrarily. Thus, we can state that vector w X 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 B X : where matrix Q col X consists of last p columns of Q X . In conclusion, we can say that vector w X can be determined uniquely when matrix B X 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 B X 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.

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 fourwheel 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 W 1 and W 2 can change their inclinations with respect to the robot platform 1, whereas edges W 3 and W 4 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 W 1 and W 2 in a synchronized manner. Due  A global inertial reference frame π 0 = x 0 y 0 is established on the ground 0, and bodyfixed local reference frames π i = x i y i are embedded in the moving bodies (Fig. 1). The absolute Cartesian coordinates that describe the mechanism form vector q: where r i = [ x i y i ] T represents the position of the local reference frame π i origin with respect to the global frame π 0 , and ϕ i is the angle of the local frame π i rotation with respect to the global frame π 0 . 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 x 0 axis, firstly turns right, then travels forward, and finally turns left and starts to move along a line parallel to x 0 .
The holonomic constraints imposed on the modeled multibody system can be written jointly as a set of m = 17 nonlinear algebraic equations: Formulas for the constraint Jacobian matrix Φ q and for vector Φ t (see Eq. (2)) as well as for Γ 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: where v (0) K is the linear velocity at contact point K, and n (0) K is a unit vector perpendicular to the allowed direction of motion (the components of both vectors are expressed in the global frame π 0 ). The other symbols (i, Ω, R i , s (i) K ) are explained in the Appendix (note that R T i ΩR i = Ω). 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), Γ is a scalar given by the following formula:

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 Φ 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: Appending the Pfaffian constraints matrix to the Jacobian matrix gives the constraint matrix C: Let q 0 denote the vector of coordinates representing the mechanism configuration shown in Fig. 1: where a = 0.05 (m), d x = c x /2, d y = c y /2 − b y , and quantities c x , c y , b y are defined in the Appendix.
It can be calculated that rank Φ q q 0 = 17, rank Ψ q 0 = 4, r= rank C q 0 = 20 < 17 + 4. (35) Matrix Φ q is a full-rank matrix, whereas Ψ 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 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: . . .
Next, Eq. (27) was utilized to calculate eleven matrices B X (note that degree of redundancy p = 2, thus only two last columns of Q X matrices were used to build corresponding Q col X matrices): Matrices corresponding to knife-edge kinematic pairs W 1 -W 4 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.

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 M 21×21 = diag(m 1 , m 1 , J 1 , . . . , m 7 , m 7 , J 7 ), and with terms C and Γ 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 q 0 (see Eq. (34)) and initial velocityq 0 : 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 overconstrained mechanism with frictionless kinematic pairs, the individual reactions are nonunique 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.

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 Γ 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 W 4 , thus reaction in this pair was set to zero as well. The alternatively eliminated constraint equations are associated with kinematic pairs C, W 1 , W 2 , and W 4 . It is worth noting that in each case elimination affected, among others, reactions in pairs B and W 3 . 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.

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- 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 overconstrained 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 de-spite the redundancy of constraints. The reactions solvability analysis methods, previously proposed for holonomic systems, are generalized to mechanisms subjected to linear nonholonomic 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 P i = Q j (for example, in Fig. 1 joint D is formed by bodies 2 and 4, P 2 = Q 4 = D). The constraint equations describing this joint can be derived by requiring that the point P i on body i coincides with point Q j on body j :  where s (k) U is the position vector of point U k in the local reference frame π k and R k is the direction cosine matrix transforming quantities from π k to π 0 : Let us consider a translational joint, formed by bodies i and j , and assume that point P i belongs to body i, and point Q j as well as vector v ij belong to body j . Points P i and Q j lie on a line parallel to the axis of relative joint motion, and vector v ij is perpendicular to this axis. The translational joint can be described by two scalar constraint equations. The first one represents the fact that vector P i Q j is perpendicular to vector v ij , and the second equation represents the fact that body i does not change its orientation with respect to body j : where ψ ij 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 i and Q j , at every time instant, is given by a scalar, time-dependent function d(t): where u ij 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 The nonzero Jacobian matrix entries corresponding to Eq. (42) are the following: 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: Vectors Γ (see Eq. (9)), representing revolute and translational joints as well as the driving constraints, can be respectively expressed as