A projection continuation approach for minimal coordinate set constrained dynamics

The formulation of constrained system dynamics using coordinate projection onto a sub-space locally tangent to the constraint manifold is revisited using the QR factorization of the constraint Jacobian matrix to extract a suitable subspace and integrating the evolution of the QR factorization along with that of the constraint Jacobian matrix, as the solution evolves. A true continuation algorithm is thus proposed to track the evolution of the subspace of independent coordinates. It does not visibly affect the quality of the solution but avoids the artiﬁcial algorithmic irregularities or discontinuities in the generalized velocities that could otherwise result from arbitrary reparameterizations of the coordinate set. The characteristics of the proposed subspace evolution approach are exempliﬁed by solving simple single-and multi-degree-of-freedom problems.


Introduction
In mechanical system dynamics, unconstrained dynamics problems are usually formulated as a set of second-order ordinary differential equations (ODE) that depend on a corresponding set of coordinates. A convenient approximation to describe the interaction between the parts of the system is often their idealization at a purely kinematic level, as algebraic relationships between the coordinates of those parts. The addition of these algebraic equations turns the problem into a system of differential-algebraic equations (DAE). The original coordinates are no longer independent; the actual number of independent coordinates reduces to that of the truly independent degrees of freedom. The enforcement of the constraints results in constraint reaction generalized forces, usually formulated as Lagrange multipliers. They represent the (unknown) internal forces whose value is whatever is required to guarantee the constraints' enforcement. A review of the possible approaches can be found in [1].
The constrained dynamics problem can be solved either directly, as a system of DAEs, where the original coordinates, augmented by the Lagrange multipliers, represent the unknowns in the so-called redundant coordinate set (RCS) formulation or, through manipulations that will be detailed in a later section, it can be transformed into the corresponding underlying ODE problem, reducing the set of coordinates to the truly independent (Lagrangian) ones, following the so-called minimal coordinate set (MCS) formulation. For a review of the possible approaches, see for example [2]. This paper focuses on this latter approach.
It is worth noticing that a third approach is possible, i.e., to somehow embed the constraints in the unconstrained problem, formally preserving its original structure and unknowns. See for example the so-called augmented Lagrangian approach [3] or the force projection method [4], as discussed in [5]. In this case, the problem formally reduces as well to ODE, with all the related implications, opportunities, and limitations in terms of approaches for its numerical integration: explicit methods can be used, subjected to conditional stability limitations.
It is recognized that the reduction of the original coordinates into the MCS may be a challenging task and that their nature is local, i.e., there may not exist a generally valid choice, which works for all configurations of the system [6]. Such a choice being local, when the coordinates need to be redefined with respect to a new configuration, discontinuities are expected in the generalized coordinates (specifically, in their derivatives), although such discontinuities are not related to any physical discontinuity in the kinematics or dynamics of the system. Indeed, the motion with respect to the original coordinates is not expected to show any discontinuity, the latter being mere artifacts of a redefinition of the local coordinates across time steps.
As originally discussed in [7], this work presents a well-known and effective method for selecting a subspace of independent coordinates that is intrinsically tangent to the constraint manifold at a specific configuration, based on the QR factorization of the constraint Jacobian matrix, and discusses how to operate the redefinition of the coordinates' subspace in a continuous manner, whenever possible, to eliminate those unnecessary, formal discontinuities.
What this work does not is advocate any superiority of the MCS over other approaches, nor any superiority of the proposed coordinate selection algorithm over others, especially in terms of performances. Being at its core a mere recombination of the coordinates, it does not substantially impact the effectiveness nor the efficiency of the solution, the required or allowed time step, nor the computational time, with the possible exception of the additional operations needed by the proposed coordinate selection, which likely add some minimal computational burden.

Constrained dynamics problem formulation
A generic constrained system dynamics problem is formulated by adding m (holonomic in the present case, without excessive loss of generality, and ideal) independent kinematic constraints as the set of algebraic equations with c ∈ R m , to a set of n ordinary differential equations (n ≥ m, but usually n > m) that express the dynamics of an unconstrained system of n coordinates x ∈ R n , subjected to a set of generalized forces f ∈ R n , energetically conjugated to a virtual perturbation of the coordinates δx, where M ∈ R n×n is the symmetric, positive-definite mass matrix. These equations are modified by the addition of the constraint reactions f c = −c T /x λ as follows: where c /x = A ∈ R m×n is the partial derivative of the constraint equations c with respect to the coordinates x, namely the constraint Jacobian matrix, which is expected to be full-rank thanks to the previously assumed independence of the constraints, and λ ∈ R m is the vector of the corresponding Lagrange multipliers. The rank of the constraint Jacobian matrix could reduce in case singular configurations are reached. This condition would be critical irrespective of the formulation in use, and thus it is not specific to the present discussion. For this reason, it is not explicitly treated.

Minimal coordinate set approach
The minimal coordinate set approach consists in defining a suitable subspace T ∈ R n×(n−m) of the space spanned by the coordinates x that is tangent to the constraint manifold, namely where q ∈ R n−m are local, truly independent coordinates, with β and β defined accordingly, the former being nonzero only in case of rheonomous constraints such that as c /x T = AT ≡ 0 holds, and analogously Suitable expressions of β and β are determined later. The constrained dynamics problem, projected in such a subspace, yields The solution is sought by first integrating Eq. (7) to obtain the generalized velocitiesq; then, Eq. (4a) is integrated to obtain an estimate of x, which needs to be subsequently refined by enforcing the constraint at the position level, Eq. (1). The cancellation of the term T T A T in Eq. (7) is guaranteed by the constraint ideality assumption. When this is not the case, e.g., in the presence of friction, the Lagrange multipliers need to be computed from local equilibria. This problem is common to all MCS formulations and not specific to the present discussion.

QR factorization for optimal minimal coordinate set selection
Among the several approaches proposed in the literature [8], a suitable choice for matrix T is obtained through the QR factorization [9] of the transpose of the constraint Jacobian matrix where matrix Q ∈ R n×n is orthogonal and submatrix R 1 ∈ R m×m is upper triangular. Submatrix Q 2 ∈ R n×(n−m) , although not uniquely determined, represents an optimal choice for T, as discussed in the following. The velocities can then be expressed aṡ with Q 1 p = β , such that which yields whereas the accelerations can be expressed as with Q 1 p = β , such that which yields Notice that the inversion of matrix R 1 is straightforward since it is triangular and nonsingular when the constraints are independent; thus, only a back substitution is needed. According to Eqs. (7) and (9), the problem becomes a form that resembles the one originally devised by Maggi [10,11] and subsequently reproposed by Kane [12], nowadays known as Maggi-Kane equations.
Its integration from time t k to t k+1 yields where the superscript (·) (0) indicates an estimate of the final value, pending verification that it complies with the constraint of Eq. (1). The final value of the unknown p results from the iterative solution of namely i.e., is the constraint Jacobian matrix at time t k+1 during the ith constraint enforcement iteration, evaluated as a function of x (i) k+1 . Operator += signifies the incremental update of the left-hand side by way of the right-hand side.

Other choices for subspace selection
Other methods have been conceived to determine a suitable subspace of the unconstrained coordinates space that is intrinsically tangent to the constraint manifold. Two of them are reported here for completeness. The continuation algorithm proposed in Sect. 2.5 could be easily adapted to them.

Singular value decomposition (SVD)
The singular value decomposition [13] decomposes a matrix A T into three matrices: where matrices U and V are orthonormal (namely U T U ≡ I, V T V ≡ I), whereas submatrix 1 is square and diagonal and contains the singular values of A on the main diagonal. Submatrix U 2 plays the role of the projection matrix T, as originally presented in [14].

Zero-eigenvalue theorem Consider the spectral representation of matrix A T A,
where the eigenvalue and eigenvector matrices U and have been partitioned to isolate the blocks related to zero-valued eigenvalues from those of the nonzero ones. Submatrix U 2 plays the role of the projection matrix T [15].

Tangent subspace selection and continuation
Submatrices Q 1 and R 1 are uniquely determined 1 once matrix A is known. Submatrix Q 2 , instead, is only subjected to matrix Q's general constraint of being orthogonal, namely ×m , but otherwise undefined; the dashed lines originating from the solution at point O in Fig. 1(a) represent a possible, arbitrary choice of the directions represented by the columns of submatrix Q 2 in a 3D point mass pendulum problem. Specifically, it is defined in excess of post-multiplication by an arbitrary orthogonal matrix, P ∈ R (n−m)×(n−m) :Q 2 = Q 2 P also complies with the orthogonality requirement sinceQ T In fact, the QR factorization produces a "local" representation of the constraint Jacobian matrix; as such, the generalized coordinates associated with the subspace T = Q 2 , which do not have any specific physical meaning, represent a local reparameterization of the subspace of the coordinates that is tangent to the constraint manifold; for example, local axes x and y originating from point O in Fig. 1(b). When the QR factorization is computed at different time steps t k , if n − m > 1, then the columns of the resulting Q 2 k are completely unrelated, their resulting value being solely dictated by the internal intricacies of the QR factorization algorithm.
This work proposes a simple and intuitive algorithm that tracks the evolution of the subspace spanned by Q 2 using a form of differential "continuation" to preserve some sort of spatial continuity of the generalized coordinates q by minimizing the amount of deviation of the subspace that is intrinsically required to maintain Q 2 tangent to the constraint manifold across time steps; for example, the axes originating from point O after the solution moved there from point O in Fig. 1(c), without altering the quality of the solution.
Consider the time derivative of the transpose of the constraint Jacobian matrix in its QR factorized form:Ȧ The derivative of matrix Q may be expressed asQ = Q , where the skew-symmetric nature of matrix ∈ R n×n , namely T = − , descends from the orthogonality of matrix Q: When the problem is integrated numerically, the solution from time step t k to time step t k+1 is computed. The QR factorization at time t k yields submatrices Q 1 k and R 1 k . The generalized velocities at time t k are computed with reference to the subspace spanned by Q 2 k . After computing the solution at the new time step, the Jacobian matrix at time t k+1 , A k+1 , is known. As such, through the economy QR factorization of its transpose, submatrices Q 1 k+1 and R 1 k+1 are determined. Instead of computing also submatrix Q 2 k+1 through the full QR factorization, the proposed continuation algorithm is used as illustrated in the following. Consider MatrixṘ 1 R −1 1 is the product of two upper triangular matrices, thus it is itself an upper triangular matrix. Matrix Q T 1Q 1 = 1 ∈ R m×m is skew-symmetric for the orthogonality of matrix Q 1 ; it can be seen as 1 ( 1 ) is the strictly lower triangular part of matrix 1 , which can be obtained as 1 being upper triangular. From Eq. (23), one can show that the derivative of matrix Q, is entirely known (the top right block of the rightmost matrix contains matrix − T 21 instead of matrix 12 since 12 ≡ − T 21 owing to the skew-symmetry of ). In fact, the bottom right block of matrix should contain an unknown skew-symmetric contribution Q T 2Q 2 = 2 ∈ R (n−m)×(n−m) . Such a matrix is arbitrarily set to zero, which corresponds to requiring that the subspace Q 2 is modified as little as possible; specifically, 2 may be interpreted as the angular velocity of subspace Q 2 , the rate of reorientation with respect to itself, whereas 21 expresses the rate of reorientation with respect to Q 1 that is required to remain orthogonal to it during its evolution.
Thus, the subspace Q 2 can be integrated starting from an arbitrary value, provided it is orthogonal to the initial value of Q 1 (the value resulting from Matlab's implementation of the QR factorization was used in the numerical examples of Sect. 3), taking appropriate measures (e.g., using Munthe-Kaas' method [16]) to guarantee that the resulting matrix Q preserves orthogonality and submatrix Q 1 matches that resulting from the factorization of the transpose of the constraint Jacobian matrix. For example, for constant across a time step of duration t k+1 − t k = h, or the latter being only a first-order approximation of the former since the intrinsic skewsymmetric structure of the exponent matrix is lost. Submatrix Q 2 k+1 resulting from the proposed integration, e.g., from Eq. (28), may need to be corrected to guarantee its orthogonality with respect to submatrix Q 1 k+1 obtained from the economy QR factorization of A T k+1 ; a Gram-Schmidt reorthogonalization [9] may be used.

Results
Simple examples are analyzed to illustrate how the proposed method produces a more regular and intuitive choice of the projection subspace during the integration of the solution.
Single-degree-of-freedom problems represent a trivial case in the context of the present discussion; in fact, the corresponding submatrix Q 2 consists of a single column, which is thus exactly determined, except for its sign.

Single-degree-of-freedom problems: planar pendulum
A simple, single-degree-of-freedom example is considered to exemplify how the continuation algorithm avoids the occasional reversing of the sign of Q 2 that occurs when the QR factorization of matrix A T is blindly performed.
Consider the equations of motion of a simple point mass pendulum of mass M and length , subjected to a uniform gravity field g along the negative z direction: Mz + 2zλ + Mg = 0 (30b) where x and z are the horizontal and vertical components of the point mass position, the unconstrained coordinates in the present context, and λ is the Lagrange multiplier associated with the constraint of Eq. (30c). The corresponding constraint Jacobian matrix is The QR factorization of its transpose can be written as where matrix Q is structurally orthonormal and defined by the single parameter θ , with as a possible determination of the parameters θ and R 1 . Submatrices Q 1 and Q 2 can then be written as The velocity and acceleration vectors of the point mass are thus The projected equation of motion is One may observe that redefiningq = θ the velocity vector can be written as which is integrable, yielding It is worth recalling that when solving for θ and R 1 from Eq. (32), To extract R 1 , consider the norm of both sides 4 x 2 + y 2 = R 2 1 cos 2 θ + sin 2 θ → or where we arbitrarily chose R 1 = +2 , but the choice with the negative sign is also legitimate. To extract θ , consider the ratio of the second and first elements: The culprit lies in the fact that in Eq. (32) matrix Q could have been alternatively defined as i.e., changing the sign of its second column, submatrix Q 2 , at the only cost of no longer representing a rotation matrix (since its determinant would now be −1 instead of +1), but without impacting its ability to represent the transpose of the constraint Jacobian matrix as Q 1 R 1 , nor that of Q 2 to represent a suitable subspace tangent to the constraint manifold.
If an intermittent change of sign occurs between consecutive QR factorizations, the sign of the derivative of the generalized coordinateq also changes, resulting in unnecessary discontinuities, although harmless for what concerns the result in terms of physical variables. Consider s 1 = ±1 as the (arbitrary) sign of R 1 and s 2 = ±1 as the (arbitrary, independent from s 1 ) sign of submatrix Q 2 such that s 2 1 = s 2 2 = 1. The QR factorization of A T then becomes 2x 2z = −s 1 cos θ −s 2 sin θ s 1 sin θ −s 2 cos θ In this case, the proposed continuation algorithm yields the "angular velocity" submatrix Considering, from Eq. (37), one obtains As a consequence, matrix from Eq. (23) by construction becomes

Single-degree-of-freedom problems: three-dimensional slider crank
Here the proposed continuation scheme is applied to the simulation of the spatial slidercrank mechanism shown in Fig. 2. This is a rigid multibody benchmark proposed by IFToMM's Technical Committee for Multibody Dynamics 2 [17] (https://www.iftommmultibody.org/) and analyzed for example in [18]. The mechanism consists of a rigid crank AB of length 0.08 m, a connecting rod BC of length 0.3 m, and a rigid sliding block. The crank, connected to the ground by revolute joint A, can rotate freely from the initial position, corresponding to an angle θ = 0 rad with an initial angular velocity of 6 rad/s. The block is constrained to the ground by a translational joint D that allows it to slide along the x axis. The problem is described by 21 coordinates, r n ∈ R 3 and e n ∈ R 4 (n = 1, 2, 3), and 20 constraint equations, thus possesses one degree of freedom q 1 . During the simulation, a time step size h = 0.00001 s is considered. The motion of the slider along the x direction is compared to the result proposed by Ramin Masoudi in the mentioned website and to those obtained with the traditional implementation of the QR projection method, also with h = 0.00001 s, which is shown in Fig. 3. The corresponding crank angle is shown in Fig. 4. The resulting q 1 andq 1 are compared in Fig. 5 and Fig. 6 to those obtained with the traditional implementation of the QR projection method. One can notice that the result of the traditional   implementation differs from the proposed one only by occasional sign changes. Since the two algorithms were initialized in the same manner, using Matlab's QR factorization and the same constraint Jacobian matrix, the initial value of q 1 is negative in both solutions. With the proposed implementation it remains such, whereas with the traditional one it occasionally turns positive and then again negative, as dictated by the intricacies of the blind execution of the QR factorization algorithm.

Multi-degree-of-freedom problems: spatial pendulum
Consider a simple point mass spherical pendulum of mass M and length , subjected to a uniform gravity field g = 9.81 m/s 2 , directed along the negative z axis. Its equations of motion are ⎡ The unconstrained problem has three degrees of freedom x, y, and z, and one constraint, Eq. (49a); thus, the constrained problem has two degrees of freedom. Consequently, Q 1 ∈ R 3×1 and Q 2 ∈ R 3×2 . The constraint Jacobian matrix and its time derivative are A = 2x 2y 2z (50) The QR factorization of A T in a given initial configuration (x, y, z) = (x 0 , y 0 , z 0 ) yields where the six coefficients q ij , i = 1, 2, 3, j = 2, 3 are related by five orthogonality conditions, leaving only one undetermined parameter. Without loss of generality, let us assume that (x 0 , y 0 , z 0 ) = ( , 0, 0), which complies with the constraint equation; this yields where α 0 is an arbitrary parameter. Clearly, when α 0 = 0, the two vectors that span the subspace of Q 2 are the coordinate axes y and z.
The projected equations of motion in the initial configuration are Without loss of generality, it is assumed that (ẋ 0 ,ẏ 0 ,ż 0 ) = (0, v 0 , 0), which complies with the derivative of the constraint equation in the initial configuration in this case, one obtains as one would expect for a diagonal element of a skew-symmetric matrix, and The trajectory of the mass simulated using the explicit Runge-Kutta scheme proposed by Dormand and Prince [19] and implemented in Matlab's ode45 function with h = 0.001 s is compared to that resulting from the integration of the original DAE governing equations with the free general-purpose multibody solver MBDyn 3 [20] using a second-order accurate Fig. 7 Trajectory of the spatial pendulum's center of mass: traditional and proposed QR approach results are compared to those obtained from DAE integration using MBDyn implicit linear multistep integration method with algorithmic dissipation (asymptotic spectral radius ρ ∞ = 0.6) [21,22] and a projection method based on the QR factorization of the transpose of the Jacobian matrix performed at each time step, without any knowledge of its evolution, using ode45 with a time step h = 0.000005 s to act as a reference solution, as shown in Fig. 7. The resulting minimal set generalized coordinates q 1 and q 2 and the projected generalized velocitiesq 1 andq 2 (q i = Q 2 (:, i) Tẋ ) are compared to the results obtained from what is here termed "traditional QR method" (i.e., without subspace continuation), as shown in Fig. 8 and Fig. 9, respectively. The same time step is used for integration with both the traditional and proposed coordinate selections not only for fairness of comparison, but also because the latter merely consists of a recombination of the coordinates' subspace, without any impact on the characteristics of the equations that are integrated. One may observe that the coordinates q i resulting from the proposed method are much more regular than those resulting from the traditional QR factorization. Specifically, those resulting from the proposed method appear to be continuous and differentiable, whereas those resulting from the traditional QR factorization show discontinuities in their first derivativesq i . This is well explained by the continuity and regularity of the evolution of each column of matrix Q 2 for the proposed method, compared to the discontinuity of those resulting from the traditional QR factorization, as depicted in Fig. 10.
Furthermore, from Fig. 10 one can observe that for t = 0 the first column of matrix Q 2 corresponds to [0, 1, 0] T , i.e., the unit vector along the y axis, whereas the second column of matrix Q 2 corresponds to [0, 0, 1] T , i.e., the unit vector along the z axis, i.e., the QR algorithm chose α 0 = 0 in Eq. (53) when initializing the subspace Q 2 . Indeed, considering the initial velocity of the mass, one can observe thatq 1 (0) ≡ v 0 andq 2 (0) ≡ 0, which is consistent with the given initial conditions.

Multi-degree-of-freedom problems: spin top
Consider an axisymmetric spin top whose tip is constrained to be at unit distance from the origin of the global coordinate system, i.e., lying on a sphere of unit radius centered in the origin. This problem was recently proposed by Haug [6]. The tip of the spin top is 1 m far away from its center of mass. The problem is sketched in Fig. 11. The inertia properties of the spin top are m = 30 kg and J = diag(90, 90, 30) kg · m 2 . The initial position of the spin top is r 0 = [0, −1, 0] T m. The local coordinate system x -y -z is initially coincident with the global coordinate system x-y-z. A uniform gravity field of magnitude 9.81 m/s 2 is assumed in the negative z direction. In this model, the Euler parameters e are used to describe the rotation of the spin top; their initial value is e 0 = [1, 0, 0, 0] T . The initial velocity isṙ 0 = where f and t are the applied force and torque at the mass center, respectively. Since only gravity is applied, f = (0, 0, −9.81m) and t = 0 with   since u T 0 u 0 ≡ 1. The problem is described by seven coordinates, r ∈ R 3 and e ∈ R 4 , and two constraint equations, Eqs. (62) and (61d), thus possesses five degrees of freedom, q i (i = 1, . . . , 5).
The trajectory of the centroid resulting from the simulation using the previously mentioned method with h = 0.0001 s is compared in Fig. 12 to those obtained by integrating the original DAE system using MBDyn and by using the traditional QR method, also with h = 0.0001 s. The projection motion of q i (i = 1, . . . , 5) is compared to the results of the traditional QR method without any projection in Fig. 13, whereas that ofq i (i = 1, . . . , 5) is compared to the results of the traditional QR method without any projection in Fig. 14. Again, one can notice the much greater regularity of the coordinates and their derivatives as they result from the proposed method.

Conclusions
This paper presented a continuation algorithm for the redefinition of the subspace of minimal coordinates that is tangent to the constraint manifold. It is based on the full QR factorization of the constraint Jacobian matrix to initialize the subspace through the portion of the space defined by the orthogonal matrix Q that is orthogonal to the constraint Jacobian matrix. The economy QR factorization is then used to exactly factor the subspace in which the constraint Jacobian matrix lies, while the evolution of the tangent subspace is tracked by integrating the time derivative of matrix Q, eventually re-orthogonalizing the result to eliminate possible drift from the integrated tangent subspace. Numerical examples show that the analysis results are unchanged, but the generalized velocities no longer show the discontinuities that occasionally characterize them when the tangent subspace is recomputed without considering its value at the previous time step.