Explicit equations of motion for constrained mechanical systems with singular mass matrices and applications to multi-body dynamics

We present the new, general, explicit form of the equations of motion for constrained mechanical systems applicable to systems with singular mass matrices. The systems may have holonomic and/or non-holonomic constraints, which may or may not satisfy D'Alembert's principle at each instant of time. The equation provides new insights into the behaviour of constrained motion and opens up new ways of modelling complex multi-body systems. Examples are provided and applications of the equation to such systems are illustrated.


Introduction
One of the central problems in analytical dynamics is the determination of equations of motion for constrained mechanical systems.This problem was formulated more than 200 years ago and has been actively and continuously worked on since by many engineers, mathematicians and scientists.The initial description of constrained motion was established by Lagrange (1787).He invented the method of Lagrange multipliers specifically to handle constrained motion.Gauss (1829) introduced a general, new principle of mechanics for handling constrained motion, which is commonly called today as Gauss's Principle.Gibbs (1879) and Appell (1899) independently discovered the so-called Gibbs Appell equations of motion using the concept of quasi-coordinates.Dirac (1964) developed, using Poisson brackets, a recursive scheme for determining the Lagrange multipliers for singular, Hamiltonian systems.Udwadia & Kalaba (1992) discovered a simple, explicit set of equations of motion for general constrained mechanical systems.Their equations can deal with holonomic and/ or non-holonomic constraints that are not necessarily independent.All these investigators have used D'Alembert's principle as their starting point.The principle states that the total work done by the forces of constraints under virtual displacements is always zero.This assumption seems to work well in many situations and can be presumed as being at the core of classical analytical dynamics.Udwadia & Kalaba (2001, 2002) generalized their equations to constrained mechanical systems that may or may not satisfy D'Alembert's principle.
Though one may perhaps consider these equations to be the most general, explicit equations obtained to date that describe the motion of constrained mechanical systems, they are limited by the fact that, unlike Dirac's formulation, they cannot deal with singular mass matrices.Systems with singular mass matrices are not common in classical dynamics when dealing with unconstrained motion.As known (Pars 1979), when the minimum number of coordinates is employed for describing the (unconstrained) motion of mechanical systems, the corresponding set of Lagrange equations usually yield mass matrices that are non-singular; in fact they are symmetric, and positive definite (Pars 1979).However, singular mass matrices can, and do, arise in the modelling of complex, multi-body mechanical systems.Such occurrences are most frequent when describing mechanical systems with more than the minimum number of required generalized coordinates, so that the coordinates are then not in fact independent of one another, and are subjected to constraints.Often, such descriptions of mechanical systems that use more than the minimum number of required coordinates are helpful in setting up the equations of motion of complex systems in a less labour-intensive fashion (see Udwadia & Kalaba 1996 for several examples).At other times, one may be interested in decomposing a complex multi-body system into its constituent parts for each of which the equations of motion are known; one then wants to use these equations for each of the constituents to obtain the equations of motion of the composite system.Singular mass matrices can arise in such situations too.Thus, in general, singular mass matrices can arise when one wants greater flexibility in modelling complex mechanical systems.
In this paper, we develop general and explicit equations of motion to handle systems whether or not their mass matrices are singular.The equations are applicable to systems with holonomic constraints, non-holonomic constraints, or a combination of the two types, as well as systems where the constraint forces may or may not be ideal.We show that, in general, the motion of such systems with singular mass matrices may not always be uniquely defined.Further, we present the necessary and sufficient conditions for the equations of motions to uniquely determine the accelerations of the system.We then show that these general equations are identical to the familiar, explicit equations previously obtained by Udwadia & Kalaba (2001) when the mass matrices are restricted to being positive definite.We provide examples to show how systems with singular mass matrices can arise in the modelling of mechanical systems in classical mechanics and demonstrate how to use our equations for describing such systems.The new equations open up new ways of modelling complex multi-body systems.For, they permit one to decompose such systems into sub-systems each of whose equations of motion are known, and then combine these sub-system equations to obtain the equations of motion of the composite system in a straightforward and simple manner.

Explicit equations of motion for general constrained mechanical systems applicable to systems with singular mass matrices
Consider an unconstrained mechanical system.The motion of the system at any time t can be described, using Lagrange's equation, by M ðq; tÞ€ q Z Qðq; _ q; tÞ; ð2:1Þ with the initial conditions qð0Þ Z q 0 ; _ qð0Þ Z _ q 0 ; ð2:2Þ where q is the generalized coordinate n-vector; Q is an n-vector which is a known function of q, _ q and t; and n is the number of generalized coordinates.In this paper we shall assume, contrary to common practice, that the matrix M that describes the unconstrained motion of the system is a symmetric n by n matrix, which is, in general, semi-positive definite.By 'unconstrained' we mean here that the n coordinates, q, are independent of one another, or are to be treated as being independent of each other.
Suppose further that the unconstrained system is now subjected to the m constraints 4 i ðq; _ q; tÞ Z 0; i Z 1; 2; .; m; ð2:3Þ where k%m equations in the equation set (2.3) are functionally independent.We shall assume that the initial conditions (2.2) satisfy these m constraints.The constraints described by (2.3) include all the usual varieties of holonomic and non-holonomic constraints, and combinations thereof.
Our aim is to obtain explicit expressions for the acceleration, € qðtÞ, of the dynamical system at the time t, given: (i) its unconstrained equation of motion, as given by equation (2.1); (ii) the constraints, as given by the equation set (2.3); and (iii) the state of the system, _ qðtÞ and q(t), at time t.Assuming that the set of constraints (2.3) are smooth enough, we can differentiate them with respect to t to obtain the relation Aðq; _ q; tÞ€ q Z bðq; _ q; tÞ; ð2:4Þ where A is an m by n matrix whose rank is k, and, b is an m-vector.When the unconstrained system gets constrained, additional forces constraint forces, Q c will arise to ensure that the constraints are satisfied.Therefore, the equation of motion for the constrained system becomes M € q Z Qðq; _ q; tÞ C Q c ðq; _ q; tÞ: ð2:5Þ We note, again, that the matrix M in equation (2.5) can, in general, be singular.As stated by Udwadia & Kalaba (2001), the work done by the forces of constraints under virtual displacements at any instant of time t can be expressed as w T Q c ðq; _ q; tÞ Z w T C ðq; _ q; tÞ; ð2:6Þ where C ðq; _ q; tÞ is an n-vector describing the nature of the non-ideal constraints, which is determined by the mechanician and could be obtained by experimentation and/or observation (see Udwadia & Kalaba 2001 for details).The virtual displacement vector, w(t), is any non-zero n-vector that satisfies (Udwadia & Kalaba 1996;Udwadia et al. 1997) Aðq; _ q; tÞw Z 0: ð2:7Þ where h is an arbitrary n-vector.Equation (2.14) is the general, explicit equation of motion for constrained mechanical systems with non-ideal constraints.Here, we do not have any restriction regarding the positive definiteness of the symmetric mass matrix M. It is allowed to be singular.We state this as our first result.
Result 1.The general equation of motion of a constrained mechanical system described by relations (2.1) (2.3), whether or not the matrix M that arises in the description of the unconstrained motion of the system is singular, is given by where M is given by relation (2.13), and h(t) is an arbitrary n-vector.
We note that, in general, because of the second member on the right-hand side of equation (2.14), the acceleration of a system with a singular mass matrix is not necessarily unique.However, when the (mCn) by n matrix M has full rank (rankZn), this second member in equation (2.14) vanishes for then and so the equation of motion becomes unique.Hence, we have the following important result.
Result 2. When the matrix M has full rank, the equation of motion of the constrained system is unique and is given by the equation Next, we provide a necessary and sufficient condition for the matrix M to have full rank.

The necessary and sufficient condition for M to have full rank
When M has full rank, the equation of motion of the constrained system is unique and is given by relation (2.15).We present the necessary and sufficient condition for M to have full rank in lemma 3.1.
has full rank, n.Equation (2.15) gives the uniquely determined acceleration of the constrained system when M has full rank.
Corollary 3.2.A simple corollary of the above result is that when the matrix M is non-singular, the equation of motion of the constrained system is unique and is given by equation (2.15).This follows from the fact that the matrix M is then of full rank.

Remarks on Result 3
(i) In general, when the mass matrix is singular, the acceleration vector of the constrained system may not be unique.The acceleration of the constrained system is unique if and only if the matrix M has full rank.Since for most practical systems in classical mechanics a unique acceleration vector is experimentally observed, and therefore expected from our theoretical models, the rank of M may be used as a check to assess the reasonableness of a given method of modelling, or a model that is obtained, for a complex multi-body system.
(ii) Even when M has full rank, the presence of semi-definite, symmetric (singular) mass matrices in the description of the unconstrained motion of systems causes conceptual differences from the situation when the mass matrices are invertible.When the mass matrix is invertible, once Q c is determined (see equation (2.5)), the acceleration of the constrained system can be explicitly found from equation (2.5), since M is invertible.Not so, when the mass matrix is singular.Even if Q c is known, the acceleration of the constrained system cannot be directly determined from equation (2.5).
The determination of the acceleration of the constrained system requires both equations (2.4) and (2.5) to be used, as in (2.12).

Equivalence with previous results for symmetric, positive definite mass matrices
In this section, we will show that when the mass matrix M is symmetric and positive definite (and hence, non-singular), the equation of motion given by relation (2.15) reduces to that given in Udwadia & Kalaba (2001).Since M is symmetric and positive definite, M has full rank and the second member on the right in equation ( 2 has an inverse.Therefore, we have Postmultiplying equation (4.9) by B C , we have And we know that Thus, equation (4.10) yields 12Þ & As a known property of generalized inverses (Udwadia & Kalaba 1996), we have ð4:13Þ where we have denoted which can be rewritten as 3), (4.4) and (4.8), we get where the last equality follows because the matrices M 1/2 and J 1 are both invertible.Thus, equation (4.13) yields Substituting equation (4.17) in equation (2.15), we get Expanding equation (4.18) gives the relation ð4:19Þ which can be also rewritten as ð4:20Þ By equations (4.4) and ( 4.3), we obtain ð4:21Þ Using relation (4.9) in the first member on the right in equation (4.21), we then have ð4:23Þ Denoting the acceleration of the unconstrained system by aZM 1 Q, we get ð4:24Þ which is the same as the result obtained in Udwadia & Kalaba (2001).We note that when the matrix M is singular, M 1/2 does not exist, and the explicit equation of motion obtained earlier (Udwadia & Kalaba 2001), and given by equation (4.24), becomes inapplicable.

Examples
In this section, three examples shall be provided to show how systems with singular mass matrices can occur in the formulation of problems in classical mechanics and how they can be handled by equation (2.14).

(a ) Example 1
Consider a wheel of mass m and radius R rolling on an inclined surface without slipping, as shown in figure 1, with the gravitational acceleration g downwards.The angle of the inclined surface is a, where 0! a! p=2.
The system clearly has just one degree of freedom.Its kinetic energy is where I c is the moment of inertia around the centre of the wheel.
If y is the vertical displacement of the centre of the wheel as it rolls down the inclined plane, the wheel's potential energy can be simply expressed as V ZKmgy: ð5:2Þ Were we to take q and y as the independent generalized coordinates and use the Lagrangian to obtain Lagrange's equations of motion for the unconstrained system (since we are assuming the coordinates are independent), we would get the relations Note that the mass matrix M describing the unconstrained motion is singular now.This singularity is a consequence of the fact that in reality the system has only one degree of freedom and we are using more than the minimum number of generalized coordinates to describe the system, and pretending that these coordinates are independent.The advantage of doing this is that the unconstrained equations of motion, (5.5), can be trivially written down.However, the two coordinates y and q are in reality not independent of one another, since the wheel rolls down without slipping.They are related through the equation of constraint y Z Rq sin a ð5:8Þ that must be added to the formulation in order to model the dynamics of the physical situation properly.Differentiating equation (5.8) twice, we obtain ð5:9Þ so that A Z ½KR sin a 1 ð 5:10Þ ð5:11Þ Since in this problem the system is ideal and the constraint forces do no work under virtual displacements, C Z 0: ð5:12Þ By equation (5.10) we have Since the rank of is 2, M has full rank, and the equation of motion of the constrained system is unique and is given by relation (2.15), so that we have Z which is the correct equation of motion for the constrained system.Example 1 illustrates the following general, practical principle related to the modelling of many complex, multi-body engineering systems.
It may be appropriate to model a mechanical system in this case, for illustration, a single degree of freedom system in terms of more than the minimum number of coordinates required to describe its configuration in this case, in terms of the two coordinates (y, q).We proceed by pretending that these coordinates are independent, because doing this will generally lead to a simpler procedure for writing out Lagrange's equations when we are dealing with complex, mechanical systems.However, to represent the physical situation properly (the fact that coordinates are indeed not independent), we must add to the model the constraints between the redundant coordinates in this case, relation (5.8) , thereby treating it as a problem of constrained motion.The system's explicit equation of motion is then obtained through the use of equation (2.15).Whether the modelling procedure is appropriate or not that is, whether the choice of the additional coordinates used in the formulation of the problem beyond the minimum number required to describe the configuration of the system, and the requisite constraints used, are suitable to describe the dynamics correctly , is dictated by whether the matrix M has full rank (see Remarks on Result 3).When it does, the unique equation of motion for the system is correctly obtained.
For the choice that has been made in this example of the coordinates (y, q), and of the constraint (5.8), this is so, and we obtain the correct equations of motion, as shown in (5.17).
It should be noted that the process of obtaining a Lagrangian could be the most difficult, tricky task in having the correct equation of motion of a constrained (especially complicated) system.By employing more than the minimum number of required coordinates properly, one can significantly simplify the process of deriving the Lagrangian and the equation of motion for the unconstrained system.Then, one can put the information of the unconstrained system and constraints in equation (2.14) and let a computer do the rest for us, determining the equation of motion for the constrained system.

(b ) Example 2
Were we to add a zero mass particle (a virtual particle) moving in the x-direction to the system given in example 1, we would have The last equation simply says that there can be no force exerted in classical mechanics on a particle whose mass is zero.
We note now that the matrix M no longer has full rank.The equation of motion is therefore now given by relation (2.14); it is no longer unique.We next where h is arbitrary.Equation (5.23) is essentially the same as the equation of motion (5.17) for the system without the virtual particle.The last equation in the set (5.23) simply states that the acceleration of a particle of mass zero (to which only a zero force can be applied) is indeterminate.
(c ) Example 3 A system of two masses, m 1 and m 2 , connected with springs, k 1 and k 2 , is shown in figure 2. We shall model this system by decomposing it into two separate sub-systems that is, we consider it as a multi-body system as shown in figure 3. The two sub-systems are then connected together by the 'connection constraint', q 1 Zx 1 Cd, where d is the length of mass m 1 .We use the coordinates x 1 , q 1 and q 2 to describe the configuration of the two sub-systems, and first treat these coordinates as being independent in order to get the unconstrained equations of motion.We then 'connect' the two sub-systems by imposing the constraint q 1 Zx 1 Cd to obtain the equation of motion of the composite system shown in figure 2.
Consider figure 3.For sub-system 1, the kinetic energy and the potential energy are given by where l 10 is the unstretched length of the spring k 1 .For sub-system 2, the kinetic energy and the potential energy are given by where l 20 is the unstretched length of the spring k 2 .Thus, for the two sub-systems, the total kinetic energy becomes and the total potential energy is Defining x 1 Z x 1 Kl 10 , and q 2 Z q 2 Kl 20 , by equations (5.24) and (5.25) we get Using Lagrange's equation with the Lagrangian Lðx 1 ; q 1 ; q 2 ; _ x 1 ; _ q 1 ; _ q 2 ÞZ T KV , the equations of motion for the unconstrained system can be written as . Decomposition of the multi body system shown in figure 2 using more than two cooridinates.
which can be put in the form Kk 2 q 2 2 6 4 3 7 5: ð5:31Þ Note that these equations for the unconstrained system are almost trivial to obtain.To model the system shown in figure 2 using these two separate subsystems, we connect the two sub-systems by using the constraint Differentiating this constraint twice we get By equations (5.31) and (5.32), we thus have A Z ½1 K1 0 ð 5:35Þ and b Z 0: ð5:36Þ Note that in choosing more than the minimum number of coordinates to describe the configuration of the system, and treating them as being independent, we obtain a mass matrix, M, that is singular.
We assume that the constraint is ideal so that CZ0.Since We next obtain the equations of motion for the constrained systems for three different cases.
Case 5.1.m 1 , m 2 O0 Note that although the matrix M is singular, the matrix M has full rank.Therefore, the equation of motion of the composite system is unique and is explicitly given by equation (2.15).Our modelling of the two degree of freedom system shown in figure 2 in terms of the three chosen coordinates is appropriate.
We now obtain the equation of motion for the constrained system as Case 5.2.m 1 Z0, m 2 O0 Since in this case m 1 Z0, M given by equation (5.38) does not have full rank.However, as noted from the physical system, the equation of motion must be unique.In order to get a unique equation of motion, our attention is drawn to the need for an additional constraint, so that the resulting M matrix has full rank.Considering the location where m 1 Z0, we obtain the following constraint relation ð5:40Þ Differentiating equation (5.40) twice with respect to time t and then putting it together with equation (5.32), we get Since A has full rank,  Hence, the equation of motion for the system can be expressed as Since in this case m 1 Z0, one can imagine that the system is now composed of only the mass m 2 , and the springs k 1 and k 2 connected to it in series.To see that equation (5.47) is correct, it may be easier to see the equation of motion in terms of x 2 .By equation (5.47), the acceleration € x 2 of the mass m 2 is given by which is, as anticipated, the correct equation, where x 1 C q 2 is the total extension of both the springs k 1 and k 2 .The first equation in the set (5.47) simply states that € x 1 Z ðk 2 =k 1 Þ € q 2 , as is obvious from equation (5.40).
Case 5.3.m 1 O0, m 2 Z0 Since m 2 Z0, M given by equation (5.38) does not have full rank.But for the equation of motion to be unique, as shown before, the matrices M and M must have full rank, and so one constraint is additionally required.It is given by k 2 q 2 Z 0; ð5:49Þ because when the mass m 2 is zero, there can be no force in the spring with stiffness k 2 .Differentiating equation (5.49) twice with respect to time t and writing it together with equation (5.32) in matrix form, we obtain We note that M now has full rank, as does M .
Therefore, the equation of motion of the composite system with m 2 Z0 becomes which, as expected, is the correct equation of motion.

Conclusions
The main contributions of this paper are as follows.
(i) The explicit equations of motion that are available (Udwadia & Kalaba 2001, 2002) to date are not applicable to mechanical systems whose unconstrained equations of motion have singular mass matrices.In this paper, we present a new, general and explicit form of the equations of motion for constrained mechanical systems applicable to such systems with singular mass matrices.These equations explicitly provide the acceleration of the constrained system and apply to systems with holonomic and/or non-holonomic constraints, as well as constraints that may or may not be ideal.It may be noted that our entire development does not require, nor use, the notion of Lagrange multipliers.(ii) We show that, in general, when the mass matrix is singular, the acceleration vector of the constrained system may not be unique.As shown in one of the examples, this can occur with the addition of a zero mass particle to a mechanical system.(iii) When the mass matrix is singular, a unique equation of motion is obtained if and only if M has full rank.Since M involves the determination of generalized inverses, we show that this condition is equivalent to the matrix having full rank.Hence, for a system whose unconstrained equation of motion has a singular mass matrix to have a unique equation of motion we require that it be suitably constrained so that the columns of M are linearly independent.
(iv) As shown in the examples, mass matrices could be singular when (i) the number of generalized coordinates used is more than the minimum number necessary to describe the configuration of the system, and when (ii) we include particles with zero mass.Both conditions can be handled by using the new, explicit equations of motion.As illustrated in an example, the former situation may arise when a complex multi-body mechanical system is sub-structured into its component sub-systems.It is because of this that the new equation of motion developed here has immediate applicability to engineering systems.(v) The equation obtained herein leads to the following principle that is of considerable practical value in the modelling of complex multi-body systems met with in civil, mechanical and aerospace engineering where zero-mass particles generally do not occur, and in which we expect the equation of motion to be unique.
It may be appropriate to model a given complex mechanical system by using more than the minimum number of coordinates required to describe its configuration, pretending that these coordinates are independent, and therefore that the system is unconstrained.Doing this for complex systems will generally lead to a much simpler procedure for writing out Lagrange's equations for the unconstrained system.To model the given physical system properly, one must then add to these equations the requisite constraints between the redundant coordinates and treat the problem as one of constrained motion.The explicit equation describing the motion of the system is then given by (2.15).Whether the modelling procedure is appropriate or not that is, whether the choice of coordinates and the corresponding constraints that are used in the mathematical formulation are suitable for correctly describing the dynamics , is shown to be dictated by whether the matrix M has full rank (see Remarks on Result 3).If M has full rank, the correct, unique equations of motion will result, and the mathematical model will correctly describe the motion of the given multi-body system.
(vi) The methodology proposed herein allows complicated mechanical systems to be decomposed into smaller sub-systems that are suitably connected together through appropriate constraints to yield the composite system.The new equation can deal with all these constituent sub-systems, whether or not their mass matrices are singular, to obtain the explicit equation for the acceleration, € qðtÞ, of the composite system in a simple, systematic and straightforward manner that is readily suitable for computer implementation.Hence, new ways of modelling complex mechanical systems, hereto not available, are opened up.(vii) The general, explicit equation of motion obtained in this paper that is applicable to systems with singular mass matrices with general, holonomic and non-holonomic constraints that may or may not be ideal, appears to be first of a kind in classical mechanics.It is shown to reduce to the known and familiar explicit equations of motion (Udwadia & Kalaba 2001, 2002) when the mass matrices are restricted to being symmetric and positive definite.
. That is, M has full rank5 M has full rank.Proof.(a) We first prove that if M does not have full rank, then M does not have full rank.If full rank, there exists an n-vector ws0 such that M if M does not have full rank, then M does not have full rank.(b) We next show that if M does not have full rank, then M does not have full rank.If the matrix M does not have full rank, there exists an n-vector ws0 such that ðI KA C AÞM A " # w Z 0; ð3:3Þ which can also be written as the two relations ðI KA C AÞMw Z 0 ð3:4Þ and Aw Z 0: ð3:5Þ Solving equation (3.5), we have w Z ðI KA C AÞd; ð3:6Þ where d is any arbitrary n-vector.Substituting equation (3.6) in equation (3.4), we obtain ðI KA C AÞM ðI KA C AÞd Z 0; ð3:7Þ from which it follows that d T ðI KA C AÞM 1=2 M 1=2 ðI KA C AÞd Z 0: ð3:8Þ Since relation (3.8) can be rewritten as ðM 1=2 ðI KA C AÞdÞ T ðM 1=2 ðI KA C AÞdÞ Z 0; ð3:9Þ we must then have, by equation (3.6), M 1=2 ðI KA C AÞd Z M 1=2 w Z 0; ð3:10Þ so that Mw Z 0: ð3:11Þ Hence, by equations (3.11) and (3.5) there exists an n-vector ws0 such that MwZ0 and AwZ0, so that M wZ 0. We have thus shown what we set out to: M does not have full rank0 M does not have full rank.Combining the results proved in parts (a) and (b) above, we have thus deduced that M does not have full rank if and only if M does not have full rank.And hence, the matrix M has full rank if and only if M has full rank.& Lemma 3.1 gives us the following important result.Result 3. The equation of motion of the constrained mechanical system described by relations (2.1) (2.3) is unique if and only if the matrix :3Þ and 0€ y Kmg Z 0; ð5:4Þ which can be written in matrix form as Figure 1.A wheel rolling down an inclined plane under gravity.
b and C, which are given, as before, by equations (5.11) and (5.12), respectively.Using equation (2.14) we then obtain the constrained equation of motion of the system well defined, we have ðAM K1=2 ÞðM 1=2 ðI KA C AÞÞ Z 0; When M is non-singular, the ranks of B and V are k and nKk, respectively.