A Fumagalli

P Masarati
email: masarati@aero.polimi.it

Efficient application of Gauss' principle to generic mechanical systems

Keywords: analytical dynamics, Gauss' principle, least constraint principle, Cholesky decomposition, non-symmetric inertia matrix

This article discusses one of the solutions proposed in the literature to Gauss' principle known as principle of 'least constraint', proposing a clarifying interpretation that allows significant computational improvements when considering its application to problems with non-diagonal mass matrix. The case of non-symmetric mass matrix is briefly discussed as well.

INTRODUCTION

In computational mechanics, the importance of constrained systems modelling is highlighted by its use in many fields of application. Although the problem of formulating the laws that systems of bodies must obey during their motion has been solved many years ago by great scientists and mathematicians like Newton, d'Alembert, Lagrange, and Euler, the very same principles can be viewed from different perspectives and formulated in different, yet equivalent, manners, as pointed out by Gauss himself in reference [START_REF] Gauss | Ueber ein neues algemeines Grundgesetz der Mechanik[END_REF]. This is the case, e.g. of Hamilton's, Gauss', Jourdain's, and Gibbs-Appell's principles, and of Maggi-Kane's equations, as illustrated in many textbooks and articles on analytical dynamics [START_REF] Finzi | Meccanica razionale[END_REF][START_REF] Schiehlen | Multibody system dynamics: roots and perspectives[END_REF].

One of the reasons those principles have been formulated for is the capability to describe the motion of constrained systems without the need to explicitly take into account the details of the kinematics of the constraints.

This article discusses the principle of 'least constraint', as formulated by Gauss [START_REF] Gauss | Ueber ein neues algemeines Grundgesetz der Mechanik[END_REF], and its solution proposed in recent years by Udwadia and Kalaba [START_REF] Udwadia | A new perspective on constrained motion[END_REF][START_REF] Udwadia | Analytical dynamics[END_REF]. The latter formulation is presented in the literature as an elegant means of simulating mechanical systems with redundant or varying constraints [START_REF] De Falco | The Udwadia-Kalaba formulation: a report on its numerical efficiency in multibody dynamics simulations and on its teaching effectiveness[END_REF]. In reference [START_REF] Arabyan | An improved formulation for constrained mechanical systems[END_REF], Arabyan and Wu state that the greatest merit of the Udwadia-Kalaba formulation is that it 'is more advantageous for systematic generation of equations of motion because the user of a program based on this formulation does not have to worry about the possible redundancy of a pair of constraints'.

Although redundant constraints do not need special care, the solution proposed in references [START_REF] Udwadia | A new perspective on constrained motion[END_REF] and [START_REF] Udwadia | Analytical dynamics[END_REF] needs some care in the choice of the coordinates and in formulating the inertia properties of the model. In fact, the Udwadia-Kalaba formulation makes intense use of the square root of the mass matrix. Its computation is straightforward when the matrix is diagonal. However, when dealing with non-diagonal matrices, the computation of the square root of the matrix can be very expensive, possibly reducing the appeal of the formulation. This may happen, e.g. when modelling rigid and/or flexible multibody systems, using a consistent inertia formulation, or accounting for the rigid-body motion of bodies whose inertia tensor is not isotropic, or whose motion is not referred to the centre of mass.

This work illustrates a modification to the original solution that reformulates it in a computationally efficient manner in those cases where a non-diagonal, possibly non-symmetric mass matrix must be dealt with, while maintaining all the qualities of the original one.

In section 2, Gauss' principle of least constraint is presented, while section 3 illustrates the expression of the corresponding accelerations as given by Udwadia and Kalaba in references [START_REF] Udwadia | A new perspective on constrained motion[END_REF] and [START_REF] Udwadia | Analytical dynamics[END_REF]. A computationally more efficient formulation when dealing with non-diagonal mass matrices is presented in section 4, which helps giving a deeper insight into the meaning of Udwadia-Kalaba's solution. Section 5 illustrates how the solution of section 3 can be obtained by a constrained minimization of the Gaussian functional, augmented by the constraint equations in a Lagrangian multipliers' sense. Appendix 2 shows how the principle and the proposed solution can be applied to those cases characterized by a non-symmetric inertia matrix, illustrated by practical examples.

GAUSS' PRINCIPLE

Gauss (Carl Friedrich Gauss (1777-1855), the great German mathematician), in his 1829 paper [START_REF] Gauss | Ueber ein neues algemeines Grundgesetz der Mechanik[END_REF], formulated the principle of least constraint in mechanics, a generalization of an intuition of Maupertuis (Pierre-Louis Moreau de Maupertuis (1698-1759) was a French philosopher and mathematician who, in a dissertation presented in 1741 to the Academy of Sciences of Paris, 'Loi du repos des corps' (Law of bodies at rest), formulated an essentially philosophical principle about a quantity, generically called 'action', that tends to be minimized in all natural phenomena, thus the name of 'least action').

The principle states that the functional

G( ẍ) = 1 2 ( ẍ -a) T M( ẍ -a) (1)
is minimal, where

a = M -1 f (2)
are the accelerations caused by the external forces f applied to the corresponding unconstrained system of inertia matrix M. The functional G( ẍ) is called the Gaussian function.

A strict minimum of G( ẍ) must exist, since matrix M is positive-definite (matrix M is positive definite iff v T Mv > 0 ∀v = 0; if the matrix is semi-definite, namely v T Mv 0 ∀v = 0, then the order of the system can be reduced by eliminating the degrees of freedom with no associated inertia, resulting in a positive-definite reduced matrix). In the case of an unconstrained system, the accelerations ẍ are coincident with those of equation [START_REF] Finzi | Meccanica razionale[END_REF], and the functional is exactly zero.

If the system is constrained by ideal bilateral constraints in the form

A(x, ẋ, t) ẍ = b(x, ẋ, t) (3)
the minimization of G( ẍ) can be used to find the right value of ẍ among those that comply with equation ( 3). This constraint equation can either be the second derivative (in this case, matrix A would possibly depend on x and t, but not on ẋ) of an equation

φ(x, t) = 0 (4a) (φ /x ) ẋ = -φ /t (4b) φ /x A ẍ = -(φ /x ẋ + φ /t ) /x ẋ -(φ /x ẋ + φ /t ) /t b (4c)
describing a set of holonomic constraints, or the first derivative (In this case, matrix A could also depend on ẋ if the constraint expressed by ψ is non-linear in ẋ; this extension is attributed to Udwadia and Kalaba [START_REF] Udwadia | An alternate proof for the equation of motion for constrained mechanical systems[END_REF]) of an equation

ψ(x, ẋ, t) = 0 (5a) ψ / ẋ A ẍ = -ψ /x ẋ -ψ /t b (5b)
describing a set of non-holonomic constraints.

The principle is equivalent to d'Alembert-Lagrange equations, but presents them in a different, interesting variational light. A formal verification is given in section 5.

UDWADIA-KALABA'S FORMULATION

According to the formulation proposed by Udwadia and Kalaba [START_REF] Udwadia | A new perspective on constrained motion[END_REF][START_REF] Udwadia | Analytical dynamics[END_REF], the only acceleration ẍ that simultaneously complies with the constraints of equation (3), and minimizes Gauss' functional expressed by equation [START_REF] Gauss | Ueber ein neues algemeines Grundgesetz der Mechanik[END_REF], is

ẍ = a + M -1/2 (AM -1/2 ) + (b -Aa) (6)
where (•) + indicates the Moore-Penrose pseudoinverse of a matrix [START_REF] Moore | On the reciprocal of the general algebraic matrix[END_REF][START_REF] Penrose | A generalized inverse for matrices[END_REF], whose properties are well described in reference [START_REF] Ben-Israel | Generalized inverses: theory and applications[END_REF]. This formulation requires the computation of the square root of the inertia matrix M, and of its inverse, an expensive operation when the matrix is not diagonal, although symmetric and positive definite.

The forces applied by the constraints to the otherwise unconstrained system are

f c = M 1/2 (AM -1/2 ) + (b -Aa) (7)
The solution proposed by Udwadia and Kalaba appears interesting because its use of a pseudoinversion allows to transparently handle overconstrained systems, which are characterized by a rank-deficient matrix A. This formulation is appealing because it allows to compute the accelerations directly from the external loads as if the problem were not constrained but, unlike other approaches, it does not need to preliminarily reduce the coordinate set. Since the resulting equations are formally independent and differential, numerical integration techniques for ordinary differential equations (ODE) can be used. In this case, however, the integration may suffer from drift, since only the second derivative of the constraints is enforced, so only the accelerations are guaranteed to be compliant; constraint stabilization is thus required. This can be provided by a variety of techniques, ranging from the early work of Baumgarte [START_REF] Baumgarte | Stabilization of constraints and integrals of motion in dynamical systems[END_REF], extended by other authors (e.g. reference [START_REF] Amirouche | Regularization and stability of the constraints in the dynamics of multibody systems[END_REF]), to augmented Lagrangian [START_REF] De Jalon | Kinematic and dynamic simulation of multibody systems[END_REF] and projection methods [START_REF] Schiehlen | Multibody system dynamics: roots and perspectives[END_REF]. An interesting review is presented in reference [START_REF] Neto | Stabilization methods for the integration of DAE in the presence of redundant constraints[END_REF].

The square root of a matrix

Equations ( 6) and ( 7) require the square root of the mass matrix. This section briefly discusses the implications of computing the square root of a matrix, in view of eliminating this requirement.

The square root of a matrix H is another matrix H 1/2 that, multiplied by itself, yields the original one, namely

H 1/2 H 1/2 = H (8)
The requirements for the square root of a matrix to exist and be real are usually met by inertia matrices, typically symmetric positive-definite.

Several techniques have been developed to compute the square root of a matrix. For example, the one based on matrix diagonalization, and other, more efficient algorithms, like those based on the matrix sign function method [START_REF] Denman | The matrix sign function and computations in systems[END_REF][START_REF] Benner | Solving stable generalized Lyapunov equations with the matrix sign function[END_REF][START_REF] Stykel | A modified matrix sign function method for projected Lyapunov equations[END_REF], or on Schur decomposition [START_REF] Björck | A Schur method for the square root of a matrix[END_REF][START_REF] Higham | Computing real square roots of a real matrix[END_REF].

The spectral decomposition of the original matrix, when it is diagonalizable (a mass matrix, usually symmetric positive-definite, is) yields

H = VDV -1 (9)
where V is a complete basis for H and D is a diagonal matrix containing the eigenvalues of matrix H. This formulation can be extended to generic non-singular (it is assumed in this work that problems with singular inertia matrices can be reduced in order to eliminate the singularity) matrices; in that case matrix D will be in Jordan canonical form instead of diagonal [START_REF] Higham | Computing real square roots of a real matrix[END_REF].

The square root of a diagonal matrix consists of a diagonal matrix whose diagonal coefficients are the square roots of the coefficients of the original matrix. If H is positive definite, or at least positive semidefinite, its eigenvalues are real non-negative; in this case, the coefficients of the square root matrix are real as well. So matrix H 1/2 is

H 1/2 = VD 1/2 V -1 (10)
where the positive roots of matrix D are arbitrarily chosen for D 1/2 . The computational complexity of this algorithm goes from 12n 3 for a symmetric matrix to 21n 3 in the general case.

Similarly, since the inverse of matrix H is

H -1 = VD -1 V -1 (11)
where D -1 is a diagonal matrix whose diagonal contains the inverses of the diagonal coefficients of matrix D, the inverse of matrix H 1/2 is

H -1/2 = VD -1/2 V -1 (12)
The Denman and Beavers iterative algorithm [START_REF] Denman | The matrix sign function and computations in systems[END_REF] consists of recursively computing matrices

Y k+1 = 1 2 (Y k + Z -1 k ) (13a) Z k+1 = 1 2 (Z k + Y -1 k ) (13b)
initialized as Y 0 = M and Z 0 = I. Matrices Y and Z converge quadratically to M 1/2 and M -1/2 , respectively. The cost of this algorithm is that of several matrix factorizations, repeated as many times as needed for quadratic convergence.

Another method [START_REF] Björck | A Schur method for the square root of a matrix[END_REF] consists of computing the Schur decomposition of matrix M as

M = QSQ T ( 14
)
where Q is a unitary matrix and S is upper triangular (quasi-upper triangular when matrix M is not diagonalizable). Then, a matrix T, such that

T 2 = S (15)
can be easily and efficiently computed based only on the square root of the diagonal coefficients (the matrix square root of the diagonal blocks when matrix M is not diagonalizable) of S. Thus the square root of matrix M is given by

M 1/2 = QTQ T (16)
This algorithm has been extended in reference [START_REF] Higham | Computing real square roots of a real matrix[END_REF] to make use only of real arithmetic when applied to real matrices that have real square roots, also considering numerical stability issues [START_REF] Higham | Stable iterations for the matrix square root[END_REF]. The computational complexity of the real arithmetic version of this algorithm, according to reference [START_REF] Higham | Computing real square roots of a real matrix[END_REF], is 15n 3 for the Schur decomposition, n 3 /6 for the computation of T, and 3n 3 /2 for equation [START_REF] Denman | The matrix sign function and computations in systems[END_REF]. The overall complexity is thus 50n 3 /3. In conclusion, the formulation proposed by Udwadia and Kalaba could require the solution of a full eigenproblem, or in any case a similarly expensive matrix decomposition, with little chances to be able to exploit any mass matrix sparsity, unless matrix M is already diagonal. As a consequence, the direct application of equation ( 6) could be quite expensive.

A BETTER COMPUTATIONAL FORM

A slightly different formulation is proposed, which is shown to possess the same properties of that presented by Udwadia and Kalaba, but requires much less computational effort when matrix M is not diagonal.

What property of the matrix square root is actually required?

In Udwadia-Kalaba's formulation, the square root of matrix M is essentially used because of its property expressed by equation [START_REF] Udwadia | An alternate proof for the equation of motion for constrained mechanical systems[END_REF]. Since matrix M is symmetric positive-definite, if the rows of matrixV are normalized with unit norm, then equations ( 9) and ( 10) can be, respectively, written as

H = VDV T (17)
and

H 1/2 = VD 1/2 V T ( 18
)
which shows that both are symmetric by construction. As a consequence, in this case equation ( 8) can also be written as

(H 1/2 ) T H 1/2 = H (19)
In fact, what is basically required by the weighted leastsquares interpretation of Gauss' principle is to be able to describe equation [START_REF] Gauss | Ueber ein neues algemeines Grundgesetz der Mechanik[END_REF] in the form of the dot product of a vector by itself, namely

G( ẍ) = 1 2 ( ẍ -a) T M( ẍ -a) = 1 2 [M 1/2 ( ẍ -a)] T [M 1/2 ( ẍ -a)] (20)
This suggests that any other convenient decomposition that allows to write matrix M as the product of a matrix by its transpose could be used as well.

A very convenient decomposition is the one proposed by Cholesky [START_REF] Benoit | Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés à un systéme d'équations linéaires en nombre inférieur à celui des inconnues (procédé du commandant Cholesky)[END_REF], namely

M = U T U ( 21
)
where U is an upper-triangular matrix (this decomposition is occasionally found in the literature in the alternative form M = LL T , where L = U T is a lowertriangular matrix). This transformation was already suggested in reference [START_REF] Redon | Gauss' least constraints principle and rigid body simulations[END_REF] as a means to turn a constrained dynamics problem into a minimization problem for the Euclidean norm of a variable w = U ẍ, corresponding to a non-Euclidean norm for the variable ẍ.

The decomposition of equation ( 21) is perfectly suitable to matrix M because the matrix is positive-definite by definition, and usually symmetric, yielding

G( ẍ) = 1 2 [U( ẍ -a)] T [U( ẍ -a)] (22)
Moreover, this decomposition can also be used in equation ( 2) to efficiently compute the accelerations of the unconstrained system required in equations ( 1) and ( 6), and thus only needs to be computed once for two purposes, allowing one to 'kill two birds with a stone', a highly appreciated art when developing computationally intensive scientific software. The computational complexity of this decomposition is n 3 /3, thus more than 30 times less than any of the algorithms illustrated in section 3.1 to extract the matrix square root. Moreover, it can be dramatically reduced when matrix sparsity is exploited.

Is the proposed modification to Udwadia-Kalaba's formula fit?

In order to demonstrate that the proposed modification to Udwadia-Kalaba's formulation is correct, one needs to show that the accelerations resulting from the formula

ẍ = a + U -1 (AU -1 ) + (b -Aa) (23) simultaneously:
(a) comply with the constraints of equation ( 3

= Aa + AU -1 (AU -1 ) + (b -Aa) = Aa + (AU -1 )(AU -1 ) + b -(AU -1 )(AU -1 ) + (AU -1 ) AU -1 (Ua) = (AU -1 )(AU -1 ) + b ( 24
)
The requirement is strictly satisfied only when AU -1 is full row-rank, implying that A is full row-rank, and thus (AU -1 )(AU -1 ) + = I. However, if equation ( 24) is not satisfied not only are the constraints redundant, but also the term b is inconsistent. The use of Gauss' principle of least constraint to address inconsistently constrained systems is illustrated in reference [START_REF] Franklin | Least-squares solution of equations of motion under inconsistent constraints[END_REF]. Proof: if matrix A is not full rank, equation (3) can only be satisfied when its right-hand side, vector b, belongs to the rangespace of matrix A. Consider a singular value decomposition (SVD) of matrix

A A = [U m×r |U m×(m-r) ] r×r 0 r×(n-r) 0 (m-r)×r 0 (m-r)×(n-r) × [V n×r |V n×(n-r) ] T ( 25
)
where submatrices

U (•)×(•) ,V (
with

A * = U T m×r A = r×r V T n×r (28)
which is full rank. The problem of equation ( 27) is formally equivalent to the original one when matrix A is full rank .

Compliance with requirement (b) is achieved by showing that any other acceleration that differs from the one of equation ( 23), while complying with equation ( 3), implies a larger value of Gauss' functional, equation [START_REF] Gauss | Ueber ein neues algemeines Grundgesetz der Mechanik[END_REF]. Proof: consider an acceleration

ü = ẍ + ÿ (29)
that complies with the constraints of equation ( 3), namely

b = A ü = A( ẍ + ÿ) (30)
Since ẍ already complies with equation (3), this implies that ÿ belongs to the nullspace of A; as a consequence

A ÿ = 0 (31)
This in turn implies that

AU -1 U ÿ = 0 (32)
The pseudoinverse of a null vector is the transpose of the null vector itself

0 + = 0 T (33)
while the pseudoinverse of a non-null vector is

v + = 1 v T v v T ( 34
)
As a consequence, the pseudoinversion of equation (31) yields

(A ÿ) + = 0 T (35)
while that of equation (32) yields

0 T = U ÿ T AU -1 + (36)
under the assumption that U ÿ = 0 for ÿ = 0, guaranteed by the positive definiteness of matrix M. The evaluation of Gauss' functional of equation ( 1) with the accelerations of equation (29) yields

G( ü) = 1 2 ( ü -a) T M( ü -a) = 1 2 [(AU -1 ) + (b -Aa) + U ÿ] T × [(AU -1 ) + (b -Aa) + U ÿ] = 1 2 [(AU -1 ) + (b -Aa)] T [(AU -1 ) + (b -Aa)] + 1 2 (U ÿ) T (AU -1 ) + = 0 (according to equation (36)) (b -Aa) + 1 2 (b -Aa) T [(AU -1 ) + ] T (U ÿ) = 0 (according to equation (36)) + 1 2 (U ÿ) T (U ÿ) = G( ẍ) + 1 2 ÿT M ÿ (37)
The value of G( ü) in equation (37) can only be larger than that of G( ẍ), since the value of ÿT M ÿ can only be positive when ÿ = 0, as matrix M is positive-definite .

Is it really the same formula?

The question whether the proposed formula, equation [START_REF] Redon | Gauss' least constraints principle and rigid body simulations[END_REF], and the one by Udwadia and Kalaba, equation [START_REF] De Falco | The Udwadia-Kalaba formulation: a report on its numerical efficiency in multibody dynamics simulations and on its teaching effectiveness[END_REF], are the same arises. In fact, their proof, which inspired the one presented in this work, clearly indicates that the formula of equation ( 6) is unique.

The next paragraph shows that the two formulae are exactly the same. The Moore-Penrose pseudoinverse of a matrix H ∈ R m×n , with n > m and rank r = m, is

H + = H T (HH T ) -1 (38)
The requirement on r = m can be relaxed to r m; this means that the constraints expressed by matrix H are not independent; in this case, the pseudoinverse is defined as the limit

H + = lim ε→0 H T (HH T + εI) -1 (39)
and this limit exists for every matrix H. So, in the case of equation ( 6), the expression

M -1/2 (AM -1/2 ) + yields M -1/2 (AM -1/2 ) + = M -1/2 (AM -1/2 ) T [(AM -1/2 )(AM -1/2 ) T ] + = M -1/2 (M -1/2 ) T A T [AM -1/2 (M -1/2 ) T A T ] + = M -1 A T (AM -1 A T ) + (40)
In the case of equation ( 23), the expression U -1 (AU -1 ) + yields

U -1 (AU -1 ) + = U -1 (AU -1 ) T [(AU -1 )(AU -1 ) T ] + = U -1 U -T A T (AU -1 U -T A T ) + = M -1 A T (AM -1 A T ) + (41)
The two expressions in equations ( 40) and (41) are clearly the same. So the use of Cholesky's decomposition simply represents a much more computationally efficient way to write the same correction to the unconstrained accelerations proposed by Udwadia and Kalaba.

WHERE DOES IT COME FROM?

The principle of least constraint requires the minimization of the functional of equation ( 1), subjected to the linear constraint of equation ( 3). Consider the functional of equation ( 1) augmented by the linear constraint of equation ( 3) using the Lagrange multipliers approach

G( ẍ, λ) = 1 2 ( ẍ -a) T M( ẍ -a) + λ T (A ẍ -b) (42)
and minimize it by imposing the vanishing of its gradient with respect to both the unknown accelerations ẍ and the Lagrange multipliers λ

∂ G ∂ ẍ = M( ẍ -a) + A T λ = 0 (43a) ∂ G ∂λ = A ẍ -b = 0 (43b) Solve equation (43a) with respect to ẍ ẍ = a -M -1 A T λ (44)
and use it to solve equation (43b) with respect to

λ λ = -(AM -1 A T ) + (b -Aa) (45)
where the pseudoinverse of matrix AM -1 A T is used instead of its inverse, to indicate the fact that the inverse either exists as is, or in the limit sense of equation (39). This allows to account for the case of redundant constrains, where the rank of matrix A ∈ R m×n is r m. When the multipliers are finally eliminated by exploiting equation (45) in equation (44), one obtains

ẍ = a + M -1 A T (AM -1 A T ) + (b -Aa) (46)
which is Udwadia-Kalaba's formula after removing the need to compute the non-strictly necessary matrix M 1/2 , and the pseudoinversion is confined to the place where it could be required, if ever, as indicated in equations ( 40) and (41). If r < m, the multipliers are not determined; this is not an issue, because the dynamics of the system are not affected. In fact, the constraints of equation ( 3) are ideal, and thus the constraint reaction forces do not do work for any admissible motion. The pseudoinversion of matrix AM -1 A T computes the minimal norm value of the multipliers, weighted by matrix M -1 , among the values that yield the accelerations of equation ( 46).

Note that equation ( 46) can be conveniently reworked as

ẍ = (I -P)a + M -1 A T (AM -1 A T ) + b ( 47
)
where

P = M -1 A T (AM -1 A T ) + A (48)
is the non-orthogonal projection (P 2 = P, P T = P) that maps the accelerations in the space of the physical coordinates onto the rangespace of the constraints A, weighted by the inertia matrix M. An analogous interpretation was given in reference [START_REF] Tseng | Efficient numerical solution of constrained multibody dynamics systems[END_REF], in a different context.

The physics of the projection provided by matrix P, interpreting equation (48) right to left, is as follows:

(a) matrix A computes the accelerations a in the rangespace of the constraints; (b) matrix (AM -1 A T ) + is the mass matrix projected in the space of the constraints, so it computes the inertia forces in the space of the constraints resulting from the previously described accelerations; (c) matrix A T projects the inertia forces in the space of the constraints back onto the space of the physical coordinates; (d) finally, matrix M -1 extracts the accelerations from the previously described inertia forces.

Since P is a projection matrix, the norm of the constrained accelerations, ẍ , is always less than the norm of the unconstrained ones, a , when b = 0, as can be shown by considering a spectral decomposition of matrix P.

It is worth noticing that the symmetry of matrix M is not required by equation ( 46), nor anywhere in this section. Although unusual, there might be cases where the equations of motion of a system can conveniently be written in a manner that yields a non-symmetric mass matrix, as pointed out in Appendix 2. What is described in this section applies in those cases as well.

COMPUTATIONAL COMPLEXITY

The claim for more efficiency comes from a series of considerations, reported in the following. They are based on the assumption that full matrices are considered. This is not usually the case for typical mass matrices related to the physical coordinates of systems of bodies. Apart from the trivial case of diagonal mass matrix, typically the mass matrix of a system of rigid bodies is block diagonal. Coupling might be more complex when inertia is formulated in a consistent manner, although typically resulting in a banded matrix. In any case, a discussion addressing issues related to sparsity handling would probably exceed the scope of this work.

The key step required to directly compute the constrained system accelerations ẍ is represented by the inversion of the matrix that yields the correction to the unconstrained system accelerations. The proposed formulation of equation ( 23) exploits an efficient factorization of matrix M that is required anyway by equation ( 2). The original formulation of equation ( 6) needs matrix M -1/2 ; its efficient computation would provide a relatively efficient factorization of matrix M as well (e.g. the spectral decomposition of M, or its Cholesky's factorization for the first iteration of Denman and Beavers' algorithm, or its Schur's decomposition).

The matrices on the left-hand sides of equations ( 40) and ( 41) do not actually need to be computed; in fact, what needs to be computed is their product by vector (Aab). As a consequence, most of the computational complexity is represented by computing and pseudoinverting the respective matrices (AM -1/2 ) and (AU -1 ), while the remaining operations would essentially be O(nm) at most.

The pseudoinverse of the previously mentioned matrices can be computed using the SVD. Its computational complexity amounts to 4n 2 m + 8nm 2 + 9m 3 . With respect to pseudoinversion, the two approaches would have the same complexity; the one proposed in this work would still be slightly more efficient since the cost of AU -1 is n 2 m, while that of AM -1/2 is 2n 2 m if M -1/2 is available as a n × n matrix, or higher if it is available as any of the factorizations indicated previously.

It is worth noticing that the 'regularized' Lagrange's formulation of equation ( 46) outperforms both the above mentioned formulations. This is especially true when n m, and thus a problem with many bodies has a limited number of constraints. In fact, equation (46) requires the pseudoinversion of matrix (AM -1 A T ). The computation of this matrix can exploit the Cholesky's factorization of matrix M, and the symmetry of the resulting matrix, for a complexity of n 2 m + nm 2 . The complexity of the SVD for the pseudoinversion reduces to 12m 3 , thanks to the symmetry of the matrix.

CONCLUDING REMARKS

This work illustrates how the solution proposed by Udwadia and Kalaba to the problem of constrained system dynamics, based on Gauss' least constraint principle, does not actually need to deal with the square root of the mass matrix, nor it requires the mass matrix to be symmetric.

A computational approach that preserves the original formalism of Udwadia and Kalaba is proposed, based on Cholesky's decomposition. Its advantages in terms of computational complexity are highlighted.

An alternative, generic projection that directly computes the constrained system accelerations from the active forces is also presented. The latter formulation does not require the symmetry of the mass matrix.

Points of strength are:

(a) no coordinate partitioning is required; (b) the problem is naturally written in the form of ODEs, thus it can be solved by explicit integration schemes; (c) redundant and inconsistent constraints can be handled seamlessly.

Points of weakness are:

(a) the solution resulting from the integration of ẍ may suffer from drift, and needs constraint stabilization; (b) the formulation needs the second derivative (for holonomic constraints) or the first derivative (for non-holonomic ones) of the constraint equation, which might not be trivial to formulate and evaluate; (c) the matrix pseudo inversion required by the proposed formula can be computationally expensive, although less than the original one especially for large bodies to constraint ratios; (d) the basic formulation does not handle non-ideal constraints; extensions have been proposed, limited to the case of independent non-ideal constraints, but no definitive solution exists, to the authors' knowledge.

Examples of non-symmetric mass matrices resulting from problem formulations of general usefulness, including rigid and deformable body dynamics, are presented and discussed in Appendix 2.

w ∈ R n normalized accelerations (w = U ẍ) x ∈ R n coordinates

Non-Symmetric Inertia Matrix

To further generalize the discussion of Gauss' principle of least constraint, consider an inertia matrix M , * that, although non-singular, is not symmetric. Nonsymmetric inertia matrices can result from Petrov-Galerkin approaches where different trial and test functions are used; examples are given at the end of this section. Since in Lagrangian mechanics an inertia matrix M results from a consistent definition of the kinetic coenergy of a mechanical system, it must be symmetric. However, there might be reasons to express the equations of motion of the system using a set of trial variables, x * , that differ from the test variables x. They are related by the invertible relationship

x * = x * (x) ( 49
)
whose second-order derivative is linear in the secondorder derivative of x

ẍ * = ∂x * ∂x ẍ + ∂ ∂x ∂x * ∂x ẋ ẋ = X -1 ẍ + (x * /x ẋ) /x ẋ (50)
Since the transformation is invertible, matrix X is nonsingular. The non-symmetric mass matrix results from this transformation, namely

M , * = MX (51)
If equations (50) were known, the symmetry of the matrix would be restored by substituting it into the definition of the fraction of the inertia forces that depend on the second derivative of x * , namely

M , * ẍ * -MX(x * /x ẋ) /x ẋ = MXX -1 ẍ (52)
and thus

M = M , * X -1 = X -T (M , * ) T (53)
Similarly, the same fraction of the inertia forces could be rewritten as generalized forces conjugated to a virtual variation of x * δx * = x * /x δx (54) by premultiplying them by the transpose of the virtual displacement δx, namely

δx T M , * ẍ * = (δx * ) T X T M , * ẍ * (55)
the corresponding inertia matrix would be

M * , * = X T M , * = (M , * ) T X (56)
To be consistent with the relationship between the trial and the test functions of equation ( 49), the augmented Gaussian of equation ( 42) needs to be redefined as

G( ẍ * (x), λ) = 1 2 ( ẍ -a) T M , * ( ẍ * -a * ) + λ T AX A , * ẍ * -b ( 57
)
with

a * = (M , * ) -1 f + (x * /x ẋ) /x ẋ = X -1 a + (x * /x ẋ) /x ẋ (58)
Recalling equations ( 56) and (50), equations (43) become

M , * ( ẍ * -a * ) + A T λ = 0 (59a) A , * ẍ * = b (59b)
Equations ( 59) are identical to equations (43) when equations ( 50) and (58) are taken into account. After substituting ẍ * from equation (59a) into equation (59b) to compute the multipliers λ, one obtains

ẍ * = a * -(M , * ) -1 A T (A , * (M , * ) -1 A T ) + (A , * a * -b) (60)
which, premultiplied by X, would yield equation (46). Note that matrix M does not need to be actually built from equation (53); in general, this is not possible at all, since matrix X might be unavailable. In this case equation (60) needs to deal with a non-symmetric mass matrix, thus requiring a linear solver for generic matrices. For typical multibody problems, an efficient sparse LU factorization would fit. Few examples follow, which illustrate how nondiagonal, non-symmetric inertia matrices can arise from the formulation of non-trivial constrained dynamics problems.

Example: Newton-Euler equations

Consider the Newton-Euler formulation of the equations of motion of a rigid body m ẍCM = f (61a)

J CM ω + ω × J CM ω = m CM ( 61b
)
where m is the mass, x CM is the location of the centre of mass, f are the external forces, J CM is the inertia tensor with respect to the centre of mass, projected in an inertial reference frame and thus depending on the parameters used to describe the orientation of the body, ω is the angular velocity of the body, and m CM are the external moments, including those of the external forces with respect to the centre of mass. Consider now an arbitrary parametrization g of the rotations, such that

R = R(g ) (62a) ω = (g ) ġ (62b)
where R is the orthonormal rotation tensor and ω is the angular velocity, while matrix (g ) expresses the relationship between the angular velocity and the derivatives of the rotation parameters, with (The operator (•)× is used to indicate the matrix that produces the vector cross-product when multiplies a vector: a× produces a × b when multiplied by vector b)

ω× = ṘR T (63)
The equations of motion, in the variables x CM and g , become

mI 0 0 J CM (g ) ẍCM g = f m CM -J CM ˙ (g ) ġ -ω × J CM ω ( 64
)
where the inertia matrix is clearly non-diagonal and non-symmetric because of the J CM (g ) block-diagonal term. In this case the transformation matrix X is known; for each body

x * = x CM g (65) M = mI 0 0 J CM (66) X = I 0 0 (g ) (67)
A symmetric inertia matrix would result when applying the Virtual Work Principle, Lagrange's equations of the first kind or similar formulations. In that case, the moment equations would be premultiplied by ω T / ġ = (g ) T , resulting in

mI 0 0 (g ) T J CM (g ) ẍCM g = f (g ) T (m CM -J CM ˙ (g ) ġ -ω × J CM ω) (68)
In this latter case, the inertia matrix is symmetric, but the moment equations no longer truly express the equilibrium of the moments applied to the body; they rather express the equilibrium of the generalized forces conjugated to the orientation parameters g . Their linearization is further complicated by the need to take into account the (g ) T term that premultiplies the moment equations. Also for this reason, implicit formulations of the equations of motion that would in any case yield a non-symmetric Jacobian matrix may find the use of Newton-Euler equations appealing.

Example: finite volume beam

Consider the finite volume beam formulation proposed in [START_REF] Ghiringhelli | A multi-body implementation of finite volume beams[END_REF]. It basically consists of writing the force and moment equilibrium of finite portions of beam with respect to the node they are related to, including the internal forces and moments at the boundaries of the finite portion. The equilibrium of an infinitesimal beam slice subjected to distributed loads yields

f m /ξ + 0 0 p /ξ × 0 f m = τ µ ( 69
)
where f and m are the internal force and moment vectors, respectively, while τ and µ are the distributed forces and moments, and p(ξ ) is the position of the generic point at abscissa ξ along the reference line p of the beam. An interesting weak form is obtained by integrating equation (69) weighted by piecewise constant functions. Integration by parts shows that this corresponds to writing the equilibrium of finite portions of beam, subjected to the external loads and to the internal forces and moments evaluated at the span-wise boundaries of each finite portion. The latter, in turn, can be expressed as functions of the beam kinematics by means of appropriate interpolation functions.

The distributed inertia forces and moments τ and µ acting on the node n, located at ξ n and comprised between abscissae a and b of each beam portion, yield They result from multiplying the reference linear and angular accelerations of the beam section, respectively, p and ω, by the sectional inertia properties (mass m and first and second inertia moments s and J per unit length, the latter evaluated with respect to the reference line, which need not be coincident with the centre of mass). In this case, equation (70) yields the nodal inertia forces and moments. For example, when the linear and angular accelerations at the generic point of abscissa ξ are linearly this expression naturally leads to a non-symmetric form for the inertia matrix of the corresponding finite element of beam. In fact, even when s = 0 equation ( 70) clearly shows that the nodal force f only depends on the distributed inertia force, while the nodal moment m n depends on both the distributed inertia force and moment because of the cross product by the arm p(ξ )p(ξ n ).

In this latter example, the transformation matrix X that would allow to symmetrize the inertia matrix cannot be explicitly formulated in a trivial manner. For this reason, the capability to deal with non-symmetric, yet consistently formulated inertia matrices can be important.

  ); (b) minimize the functional of equation (1). Compliance with requirement (a) is ensured by the fact that, after substituting equation (23) into equation (3), one obtains b

  •)×(•) are unitary, while submatrix r×r is diagonal and contains the non-null singular values of A. Vector b must belong to the space represented by submatrix U m×r b = U m×r b * (26) b * ∈ R r . The problem can be reformulated by leftmultiplying equation (3) by the transpose of matrix U m×r , yielding A * ẍ = b *

  interpolated from nodal values, indicated by the subscript n, by means of regular shape functions N n (ξ ), namelyp = N n (ξ )p n (73a) g = N n (ξ )g n (73b)where p and g are generic position and orientation parameters, and thusp = N n (ξ ) pn (74a) ω = [g (ξ )]N n (ξ ) g n + γ[g (ξ ), ġ (ξ )] (74b)

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the fruitful discussions they had with Professors Ettore Pennestrì and Paolo Mantegazza on the subject of this work.

APPENDIX 1

Notation