Formulating the dynamics equations of a mechanical system following a multibody dynamics approach often leads to a set of highly nonlinear Differential Algebraic Equations (DAEs). While this form of the equations of motion is suitable for a wide range of practical applications, in some cases it is necessary to have access to the linearized system dynamics. This is the case when stability and modal analysis are to be carried out; the definition of plant and system models for certain control algorithms and state estimators also requires a linear expression of the dynamics.
A number of methods for the linearization of multibody dynamics can be found in the literature. They differ in both the approach that they follow to handle the equations of motion and the way in which they deliver their results, which in turn are determined by the selection of the generalized coordinates used to describe the mechanical system. This selection is closely related to the way in which the kinematic constraints of the system are treated. Three major approaches can be distinguished and used to categorize most of the linearization methods published so far. In this work, we demonstrate the properties of each approach in the linearization of systems in static equilibrium, illustrating them with the study of two representative examples.
Multibody system (MBS) dynamics has experienced a significant growth over the last decades, boosted by advances in computer architecture and software. One of the most attractive features of MBS dynamics is the wide variety of modeling and formulation approaches from which researchers and analysts can choose when dealing with a particular problem [1]. Generally speaking, the system can be modeled with a set of independent generalized coordinates, which leads to a system of Ordinary Differential Equations (ODEs), or using dependent coordinates whose motion is related by one or more kinematic constraints. In the latter case the dynamics equations are expressed with a set of DAEs, which must be handled with appropriate algorithms, e.g., [2,3]. New formulations and methods are proposed and developed by the MBS research community, making it necessary to benchmark their performance and characterize their features. The ultimate goal of this effort is to provide guidelines that help the analyst select the most appropriate formalisms when dealing with a particular problem.
Such guidelines are also necessary to select linearization methods for those applications that require a linear approximation of the equations of motion. This is the case of certain classical control algorithms, which need to represent the plant to be controlled with linear models [4], and some state and input estimators, like Kalman filters [5]. Even when a linear expression of the dynamics is not strictly required, its availability can make it simpler to gain insight into the system behavior. Using reduced order models, extracted from the linearization of the original differential-algebraic equations of motion, simplifies the analysis of structural and aeroelastic problems [6]; linearized models are also useful in stability analysis [7–9]. Modal analysis and the determination of natural frequencies and vibration modes of a system are another important application of linearized dynamics [10].
Several methods have been proposed in the multibody literature to arrive at linearized forms of the dynamics equations. They present different features, depending mainly on the selection of the coordinates with which they describe the mechanical system. Some of them are recursive algorithms for robotic systems [11], others are built on Maggi-Kane’s coordinates [12]. Lagrange multipliers are a popular way to represent constrained systems and methods to linearize the resulting equations exist as well [6, 10, 13]. This paper categorizes linearization methods for multibody dynamics into three main groups, defined by the selection of generalized coordinates. Their properties are demonstrated in the linearization of the dynamics of mechanical systems in static equilibrium, for which Linear Time Invariant (LTI) problems are guaranteed to be obtained. Representative methods of each category were applied to the linearization of test problems and compared in terms of efficiency, ease of use, and accuracy. Moreover, in [14] two types of linearization problems were identified: heavily constrained systems in which the number of degrees of freedom is much lower than the number of kinematic variables, and systems (typically flexible mechanisms) in which both numbers are of similar magnitude. The behavior of the linearization methods was studied with an example of each type. Practical guidelines for the selection of the most appropriate technique for each case are presented to the reader as conclusions.
A multibody system can be described with a set of $n$ generalized coordinates $\mathbf{q}=\{q_{1},\dots,q_{n}\}^{\mathrm{~T~}}$ related by $m$ kinematic constraints $\Phi(\mathbf{q},t)=\mathbf{0}$ . For the purposes of this paper, these are assumed to be holonomic and linearly independent; accordingly, the $m\times n$ Jacobian matrix of the constraints $\Phi_{\mathbf{q}}$ can be assumed to have full row rank $m$ . The equations of motion can be expressed as a nonlinear system of DAEs as
where $\mathbf{M}\left(\mathbf{q}\right)$ is the $n\times n$ mass matrix, $\mathbf{f}$ is the array of generalized applied and velocity-dependent forces, and $\mathbf{f}_{c}$ are the reaction forces introduced by the kinematic constraints. Equations (1) are often cast into one of the multibody formulations that exist in the literature, which can be classified into three main groups [14]:
1. Minimal Coordinate Set (MCS), in which the formulation has as many coordinates as degrees of freedom the system has, $n-m$ ;
2. Redundant Coordinate Set (RCS), consisting of the $n$ coordinates and a set of $m$ Lagrange multipliers; and
3. Unconstrained Coordinate Set (UCS), in which the system is modeled with the original $n$ coordinates of the unconstrained problem and the constraints are embedded in the formulation.
Velocity projection techniques, in which the system velocities are projected onto the subspace of admissible motion [15], fall into the MCS category. Graph-theory-based [16,17] and Transfer Matrix [18] multibody algorithms could also be included in this category.
If a Lagrangian approach is followed, the constraint reactions are expressed in terms of the constraints Jacobian matrix as $\mathbf{f}_{c}=-\Phi_{\mathbf{q}}^{\mathrm{T}}\pmb{\lambda}$ , where $\lambda$ contains $m$ Lagrange multipliers. Using Baumgarte stabilization to remove constraints drift [19] one would directly arrive at a RCS method. Alternatively, eliminating the constraints delivers a UCS approach. UCS methods can also be obtained through penalty formulations [20] or with the force projection approach in [21].
As mentioned, selecting one approach or another determines the output of the linearization methods, their efficiency and the way in which they convey information about the mechanical system. MCS formulations result in a linear problem with the exact spectrum of the constrained system; the eigenanalysis of this problem directly yields the $2\left(n-m\right)$ eigenvalues of the original system. RCS and UCS approaches may give rise to spurious or approximate eigenvalues; on the other hand, their use can be convenient for the sake of efficiency or ease of use. In the subsequent sections, representatives of each category will be compared in terms of effectiveness and efficiency, and their main features will be discussed. Two main factors, namely how easy it is to discriminate spurious eigenvalues and the simplicity of the linearized equations, also from the computational point of view, will be considered in the discussion.
In this section three methods to linearize Eqs. (1) are introduced. The first one is a MCS formulation obtained via velocity projections [22]. The second directly linearizes a RCS description of the dynamics [6]. The last one is a UCS formulation obtained replacing the kinematic constraints with penalty systems [20].
where $\mathbf{R}\left(\mathbf{q}\right)$ is an $n\times(n-m)$ velocity transformation matrix that verifies $\mathbf{R}^{\mathrm{T}}\boldsymbol{\Phi}_{\mathbf{q}}^{\mathrm{T}}=\mathbf{0}$ . Matrix $\mathbf{R}$ can be determined via several methods: e.g., coordinate selection (splitting $\mathbf{q}$ into independent and dependent variables), the zero eigenvalue theorem [23], Singular Value Decomposition (SVD) [24,25], QR decomposition [26], and Gram-Schmidt orthonormalization [27, 28]. All these methods determine the subspace $\mathbf{R}$ of the velocities that complies with the constraints through operations on the Jacobian $\Phi_{\mathbf{q}}$ . Applying the transformation in Eq. (2), the system dynamics become a system of $(n-m)$ ODEs
If the inputs are represented by $\mathbf{u}$ , the system can be linearized about an equilibrium configuration $\mathbf{z}_{0},\,\dot{\mathbf{z}}_{0},\,\ddot{\mathbf{z}}_{0},\,\mathbf{u}_{0}$ as follows
The expressions for the numerical evaluation of the partial derivatives of $\mathbf{R}$ are detailed in [29]. The linearized dynamics can then be written as
Equation (6) has $2\left(n-m\right)$ eigenvalues that represent the exact spectrum of the problem. Moreover, its leading matrix $\mathbf{M}_{r}$ is symmetric and positive-definite. As a consequence, Eq. (6) can be used in the solution of forward-dynamics problems with explicit integration schemes.
It must be noted that the evaluation of $\mathbf{K}_{q}$ requires the computation of the Lagrange multipliers $\lambda$ at equilibrium. Besides the $2\left(n-m\right)$ eigenvalues in the spectrum of the constrained problem, the generalized eigenanalysis of the $(2n+m)$ linearized equations (8) introduces $3m$ spurious eigenvalues, related to the constrained kinematic variables and the Lagrange multipliers. These are easily identifiable because their value is either infinity [6] or zero, if either or both $\Phi_{\mathbf{q}}$ and $\Phi_{\mathbf{q}}^{\mathrm{T}}$ are transferred from $\mathbf{A}_{q}$ to $\mathbf{E}_{q}$ . Because Eq. (7) represents a differential-algebraic problem, matrix $\mathbf{E}_{q}$ is structurally singular. Therefore, Eq. (8) cannot be used in forward-dynamics problems with explicit numerical integrators. It is possible though to overcome this issue applying a QZ decomposition to matrices $\mathbf{E}_{q}$ and $\mathbf{A}_{q}$ [30]. The details of this method are discussed in [6].
# 3.3 UCS Formulation: Penalty Method
Penalty-based relaxation of the constraints can be used to transform the system of DAEs in Eq. (1) into a set of $n$
ODEs [20]. This approach makes the Lagrange multipliers $\lambda$ proportional to the constraints violation at the configuration, velocity, and acceleration levels
where $\Xi$ is an $n\times n$ matrix that contains the penalty factors; $\Theta$ and $\pmb{\Omega}$ are also $n\times n$ terms with stabilization parameters that have a similar function to those used in Baumgarte stabilization [31]. The derivatives of the constraints with respect to time are $\dot{\Phi}\,{=}\,\Phi_{{\bf q}}\dot{\bf q}+\Phi_{t}$ and $\ddot{\Phi}=\Phi_{\mathbf{q}}\ddot{\mathbf{q}}+\dot{\Phi}_{\mathbf{q}}\dot{\mathbf{q}}+\dot{\Phi}_{t}$ , where $\Phi_{t}=\partial\Phi/\partial t$ . This approach is equivalent to replacing the constraints with mass-spring-damper systems. The resulting dynamics equations are
It must be noted that the equilibrium configuration of Eq. (12) is not completely equivalent to that of Eq. (1). In most cases a certain violation of the kinematic constraints $\Phi={\bf0}$ will be required to balance the applied forces f. The linearized dynamics can be expressed as
Fig. 1: An $N$ -loop four-bar linkage with spring elements
together with $2m$ spurious eigenvalues related to the constrained coordinates. True and spurious eigenvalues can be told apart by checking if their associated eigenvectors v verify the velocity-level kinematic constraints Φ˙ΦΦ = 0. The violation of such constraints remains close to zero for true eigenvalues, while it is not negligible for the spurious ones.
As it happened in Section 3.1, if the penalty matrix $\Xi$ is correctly chosen, $\mathbf{M}_{p}$ is symmetric and positive-definite, so Eq. (14) can be directly used in forward-dynamics applications.
# 4 Numerical Examples
The methods described in Section 3 were applied to the linearization of mechanical systems about static equilibrium configurations. Two examples were selected: the first was an $N$ -loop four-bar linkage with spring elements along the diagonals. The second was a flexible double pendulum. The four-bar linkage is heavily constrained and only has one degree of freedom; it is representative of mechanical systems in which the number of kinematic constraints is similar to the number of generalized coordinates $\left(n-m\ll n\right)$ . The double pendulum, on the other hand, includes degrees of freedom associated with flexibility among the generalized coordinates $\mathbf{q}$ . In this case, $n-m\approx n$ .
Figure 1 shows an $N$ -loop four-bar linkage made up of equal rods of length $l_{f}=1\;\mathrm{m}$ and uniformly distributed mass $m_{f}=1~\mathrm{kg}$ . It moves under gravity effects with $g=9.81$ $\mathrm{m}/\mathrm{s}^{2}$ . Each loop $i$ in the linkage has a spring connecting points $B_{i}$ and $A_{i-1}$ of stiffness $k_{f}\,{=}\,25\,\mathrm{N/m}$ and natural length $l_{f0}=\sqrt{2}\,\mathrm{m}$ . If $N_{\mathrm{}}=1$ , the system is equivalent to the test case discussed in [14].
The linkage is modeled with the $x$ and $y$ coordinates of points $B_{0}–B_{N}$ as variables plus the $\boldsymbol{\upvarphi}$ angle from the global $x_{\mathrm{{}}}$ - axis to the rod that connects points $A_{N}$ and $B_{N}$ $\begin{array}{r}{{n=2}N+3}\end{array}$ ). The kinematic constraints term $\Phi$ is composed by equations that enforce that the distances between the tips of each rod remain constant during motion, plus one extra equation that relates the value of $\boldsymbol{\upvarphi}$ to $x_{N}$ and $y_{N}$ $\!\!\!/m=2N+2\!\!\!$ ).
Equation (14) is a system of $n$ ODEs. The method delivers an approximation of the $2\left(n-m\right)$ true system eigenvalues,
# 4.2 Flexible Double Pendulum
Figure 2 shows a planar flexible double pendulum. It is composed of two links with length $l_{p}=4~\mathrm{m}$ and uniformly
Fig. 3: Coordinates used to describe the motion of each link
distributed mass $m_{p}$ , density $\uprho=1500\ \mathrm{kg}/\mathrm{m}^{3}$ , and Young modulus $E=10^{8}\,\,\mathrm{\dot{N}}/\mathrm{m}^{2}$ . The rods have a circular section with radius $r_{p}=0.15\;\mathrm{m}$ . Rod flexibility is modeled following the finite element approach described in [32]. Each link has a tangent floating frame of reference attached to its first articulation, point $o$ for the first link and point $P_{1}$ for the second one, and is discretized with $n_{p}$ finite elements.
The motion of each link is represented through coordinates $x,y$ , and θ, which describe the translation and rotation of the body-fixed reference frame with respect to the global one, and nodal coordinates $q_{f x}^{i},q_{f y}^{i}$ , and $q_{f\Theta}^{i}$ , which represent the elastic deformations at node $i$ , as shown in Fig. 3. The total number of generalized coordinates is $n=6n_{p}+12$ .
Regarding kinematic constraints, four equations are used to fix point $o$ to the ground and to ensure that both links remain connected at point $P_{1}$ during motion. Three more constraints per link enforce that the body-fixed reference frame remains attached to the first node of the rod and tangent to its center line, by making q1fx = q1fy $q_{f x}^{1}=q_{f y}^{1}=q_{f\theta}^{1}=0$ Thus, the total number of constraint equations is $m=10$ .
# 5 Results and Discussion
The linearization of the examples in Section 4 takes place about static equilibrium configurations, in which ${\dot{\mathbf{q}}}=\mathbf{0}$ and $\ddot{\bf q}={\bf0}$ . Under these conditions, terms in Eqs. (5), (10), and (15) adopt simpler expressions. The terms required by
The terms required by the penalty formulation in Section 3.3 become
The linearization methods were implemented in MATLAB; the computation times reported in the following were obtained running the code in an Intel i7-4790K processor at $4.0\ \mathrm{GHz}$ with 8 GB RAM under Windows $10~\mathrm{Pro}$ . These times include the evaluation of the linearization terms in Eqs. (16)–(18) and the solution of the corresponding eigenvalue problem, but not the computation of the equilibrium configuration. In the case of the penalty formulation, they also include the time elapsed in the discrimination of the spurious eigenvalues.
In all the numerical experiments reported in this section, the penalty and stabilization terms required by the penalty formulation were diagonal matrices, proportional to a single scalar value: $\Xi=\alpha\mathbf{I}$ , $\Theta=2\xi\omega\mathbf{I}$ , and $\Omega=\omega^{2}\mathbf{I}$ , where I is the $n\times n$ identity matrix.
# 5.1 Multiple Loop Four-Bar Linkage
Regardless of $N$ , the exact spectrum of the problem is composed of a double complex eigenvalue $s_{1}$ . Angle $\boldsymbol{\upvarphi}$ was selected as independent coordinate with the velocity transformation method. The values selected for the parameters of the penalty method were $\upalpha=10^{6}$ , $\xi=1$ , and ${\mathfrak{o}}=10$ . The MCS and RCS methods delivered the exact value of $s_{1}$ , while the penalty formulation introduced some error in its value as a consequence of the modification of the original system. Table 1 shows the system properties and the time elapsed in linearization with each method.
In order to evaluate the suitability of the methods when dealing with systems with dissipative elements, dampers with $c_{f}=1~\mathrm{Ns/m}$ were introduced alongside the four-bar
Table 1: Properties and spectrum of the $N$ -loop four-bar linkage, and elapsed times in linearization, for different values of $N$
The penalty formulation did not experience any problems in the singular configuration and correctly identified the two double eigenvalues $s_{1}=\pm3.4310i$ and $s_{2}=\pm4.4294\,i$ of the system without any special provisions. The generalized eigenanalysis in Section 3.2 also found the correct eigenvalues. The determination of the static solution required pseudo-inversion, since $\Phi_{\mathbf{q}}$ is no longer full-rank. Indeed, the multipliers $\pmb{\lambda}$ are now underdetermined; nonetheless, the solution remains determined. The MCS formulation, however, failed unless a new selection of independent coordinates was carried out at the singular configuration. The increase in the degree of freedom means that two independent velocities instead of one are necessary to fully define the system motion at this point.
# 5.2 Flexible Double Pendulum
In the flexible double pendulum example, the number of eigenvalues that composes the spectrum of the problem increases with the number of finite elements $n_{p}$ used in the discretization. The dynamics of this system was linearized at the equilibrium configuration in which both rods were aligned along the negative $y$ -axis, i.e., $\theta_{1}=\theta_{2}=-\pi/2$ . Table 3 shows the first four system eigenvalues computed using the MCS method.
Fig. 4: A singular static equilibrium configuration for the one-loop four-bar linkage with $k_{f}=0$
The RCS and UCS methods introduced some deviations from the eigenvalues reported in Table 3; both yielded eigenvalues with a small real component. The values of the parameters of the penalty method were $\upalpha=10^{7}$ , $\xi=1$ , and ${\mathfrak{o}}=10$ . Tables 4 and 5 show the elapsed times in the linearization of the dynamics and the maximum differences between the first four eigenvalues computed with the MCS approach $\left({{s}_{i}}\right)$ and the other methods $(s_{i}^{*})$ . It should be mentioned that, while the differences yielded by the UCS method are a consequence of the modification of the original system, the RCS should deliver exact eigenvalues. The errors reported in Table 5 are quite small and likely the result of numerical tolerances in the QZ algorithm used by MATLAB to carry out the eigenanalysis, and in the determination of $\pmb{\lambda}$ required by Eq. (17).
springs in a second set of numerical simulations. All the methods correctly evaluated the system eigenvalues for the damped case, shown in Table 2; the elapsed times were similar to those reported in Table 1.
In all cases the error in the eigenvalues introduced by the penalty formulation remained below $3\cdot10^{-6}$ rad/s. This error depends on the value of the penalty factor $\upalpha$ and the stabilization parameters $\mathbf{\omega}_{\infty}$ and $\xi$ [14].
A particular case of the static equilibrium of the four-bar linkage takes place when $k_{f}=0$ and gravity acts along the global $x_{\mathrm{{}}}$ -axis. In this scenario, the equilibrium configuration occurs at $\upvarphi=0$ , as illustrated in Fig. 4 for $N=1$ . This is a singular configuration at which the Jacobian $\Phi_{\mathbf{q}}$ instantaneously loses rank, leading to the sudden introduction of one extra degree of freedom [33].
# 5.3 Methods Assessment
The study of the two examples in Sections 5.1 and 5.2 highlighted the features of the three linearization approaches considered in this research.
The MCS velocity projection approach in Section 3.1 delivered the exact spectrum of the problems. This method requires selecting a set of independent velocities $\dot{\bf z}$ and the evaluation of the transformation matrix R, which in the general case is accomplished numerically. In problems with $n\approx m$ like the $N$ -loop four-bar linkage, composed mostly of rigid links, the method significantly reduces the number of coordinates and is computationally advantageous. However, if $n\gg m$ , as in the case of the flexible pendulum, the reduction in problem size is not so significant and may be outweighed by the evaluation of R. The method will fail if the selection of independent velocities does not remain valid during the whole motion, as in the case of systems with singular configurations.
Table 4: Linearization of the flexible double pendulum: elapsed times (s)
The RCS and UCS methods added spurious eigenvalues to the problems spectra. The direct linearization described in Section 3.2 is much simpler than evaluating the terms required by the penalty and velocity transformation methods in Eqs. (5) and (15). Moreover, the spurious terms are easily identifiable. However, the linearized dynamics obtained this way cannot be directly used in explicit timeintegration schemes. The decomposition process required to circumvent this drawback can be computationally expensive, although subspace methods, e.g., the Implicitly Restarted Arnoldi Method in ARPACK [34], can be used to alleviate this issue. Indeed, in practical applications only a subspace of the entire eigenspace is usually needed. Furthermore, without even considering the spurious eigenvalues, computing the entire spectrum in practical problems with thousands to millions degrees of freedom would be too computationally expensive and of little use, if at all feasible.
The penalty method in Section 3.3 delivers approximated eigenvalues, as a consequence of relaxing the original kinematic constraints and replacing them with penalty systems, thus modifying the equilibrium point of the system. It also introduces spurious eigenvalues in the spectrum. The parameters of the method must be adjusted [33,35] to ensure that the linearized dynamics in Eq. (14) can be considered an adequate approximation of the constrained system’s exact one.
Each approach also affects the sparsity pattern of the system matrices in a different manner. The sparsity patterns of matrices $\mathbf{M}_{r}$ , $\mathbf{M}_{q}$ , and $\mathbf{M}_{p}$ in the double pendulum example with $n_{p}=5$ are compared in Fig. 5, both in their original form and after a Cholesky factorization. Direct eigenanalysis preserves matrix sparsity. The velocity transformation of the MCS method reduces the system size, but at the cost of increasing the fill-in of the system matrix after factorization. The penalty method introduces some additional structural non-zeros, but its effect on sparsity is less severe. In practice, both the RCS and the UCS approaches preserve the structure of the mass matrix in terms of separation between the inertia coefficients of each structural component. As a consequence, the Cholesky factors of Figs. 5d and 5f produce a characteristic block-diagonal structure. It must be noted that in systems with a more complicated, non openloop topology, the constraints responsible for loop closure would introduce some coupling terms off the main block diagonal in matrix $\Phi_{\mathbf{q}}^{\mathrm{T}}\Xi\Phi_{\mathbf{q}}$ and, consequently, in matrix $\mathbf{M}_{p}$ . On the contrary, the MCS structurally couples the degrees of freedom of the components, thus completely filling the Cholesky factor of Fig. 5b.
The features of each method are summarized in Table 6.
# 6 Conclusions
When dependent coordinates are used to model a multibody system, kinematic constraints can be handled with a wide variety of methods, which can be grouped into three major categories: Minimal Coordinate Set (MCS) methods, in which the dynamics are expressed in terms of a set of independent velocities; Redundant Coordinate Set (RCS) methods, in which the set of coordinates is enlarged with as many Lagrange multipliers as constraints the system has; and Unconstrained Coordinate Set (UCS) methods, which embed the constraints in the dynamics equations. This work assessed the properties of each family of methods in the linearization of two examples in static equilibrium configurations: an $N$ -loop four-bar linkage with spring elements and a flexible double pendulum modeled with finite elements.
Fig. 5: Sparsity pattern of the mass matrix of the flexible double pendulum with $n_{p}=5$
Results confirmed that the method choice affects the linearization process in terms of the way in which information about the linearized system behavior is delivered, accuracy, and computational efficiency. While MCS methods obtained the exact problem spectrum, RCS and UCS introduced spurious eigenvalues in the spectra that needed to be properly identified and dealt with. The RCS approach featured a simple expression of the linearization terms, while the resultant descriptor form equations required additional processing before they could be used in forward-dynamics settings with explicit integration. MCS and UCS methods, on the other hand, provided compact expressions that could be readily used to this end. Regarding efficiency, the MCS method proved to be advantageous in the linearization of heavily constrained problems. Conversely, if the number of constraints was small compared to the total number of variables used in the modeling, RCS and UCS methods showed a superior performance. As a consequence, the linearization method selection must be carried out taking into consideration the characteristics of the mechanical system and the intended use of the linearized equations of motion.
# Acknowledgements
The first author would like to acknowledge the support of the Spanish Ministry of Economy through its postdoctoral research program Juan de la Cierva, contract No. JCI-2012-12376.
# References
[1] Bauchau, O. A., 2011. Flexible Multibody Dynamics. Springer.
[2] Mariti, L., Belfiore, N. P., Pennestr\`ı, E., and Valentini, P. P., 2011. “Comparison of solution strategies for multibody dynamics equations”. International Journal for Numerical Methods in Engineering, 88(7), pp. 637– 656.
[3] Marques, F., Souto, A., and Flores, P., 2016. “On the constraints violation in forward dynamics of multibody systems”. Multibody System Dynamics, Online First, pp. 1–35.
[4] Ogata, K., 2001. Modern Control Engineering, 4th ed. Prentice Hall.
[5] Grewal, M. S., and Andrews, A. P., 2015. Kalman Filtering: Theory and Practice with MATLAB, 4th ed. Wiley-IEEE Press.
[6] Ripepi, M., and Masarati, P., 2011. “Reduced order models using generalized eigenanalysis”. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics, 225(1), pp. 52–65.
[7] Escalona, J. L., and Chamorro, R., 2008. “Stability analysis of vehicles on circular motions using multibody dynamics”. Nonlinear Dynamics, 53(3), pp. 237– 250.
[8] Masarati, P., 2013. “Estimation of Lyapunov exponents from multibody dynamics in differentialalgebraic form”. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics, 227(4), pp. 23–33.
[9] Masarati, P., and Tamer, A., 2015. “Sensitivity of trajectory stability estimated by Lyapunov characteristic exponents”. Aerospace Science and Technology, 47, pp. 501–510.
[10] Masarati, P., 2009. “Direct eigenanalysis of constrained system dynamics”. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics, 223(4), pp. 335–342.
[11] Gontier, C., and Li, Y., 1995. “Lagrangian formulation and linearization of multibody system equations”. Computers & Structures, 57(2), pp. 317–331.
[12] Peterson, D. L., Gede, G., and Hubbard, M., 2015. “Symbolic linearization of equations of motion of constrained multibody systems”. Multibody System Dynamics, 33(2), pp. 143–161.
[13] Negrut, D., and Ortiz, J. L., 2006. “A practical approach for the linearization of the constrained multibody dynamics equations”. Journal of Computational and Nonlinear Dynamics, 1(3), pp. 230–239.
[14] Gonz´alez, F., Masarati, P., and Cuadrado, J., 2016. “On the linearization of multibody dynamics formulations”. In Proceedings of the ASME 2016 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, no. DETC2016-59227.
[15] Serna, M. A., Avil´es, R., and Garc´ıa de Jal´on, J., 1982. “Dynamic analysis of plane mechanisms with lower pairs in basic coordinates”. Mechanism and Machine Theory, 17(6), pp. 397–403.
[16] Jain, A., 2011. “Graph theoretic foundations of multibody dynamics. Part II: analysis and algorithms”. Multibody System Dynamics, 26(3), pp. 307–333.
[17] Jain, A., 2012. “Multibody graph transformations and analysis – Part II: Closed-chain constraint embedding”. Nonlinear Dynamics, 67(3), pp. 2153–2170.
[18] Rong, B., Rui, X., and Wang, G., 2010. “Modified finite element transfer matrix method for eigenvalue problem of flexible structures”. Journal of Applied Mechanics, 78(2), p. 021016.
[19] Baumgarte, J., 1972. “Stabilization of constraints and integrals of motion in dynamical systems”. Computer Methods in Applied Mechanics and Engineering, 1(1), pp. 1–16.
[20] Bayo, E., Garc´ıa de Jal´on, J., and Serna, M. A., 1988. “A modified Lagrangian formulation for the dynamic analysis of constrained mechanical systems”. Computer Methods in Applied Mechanics and Engineering, 71(2), pp. 183–195.
[21] Masarati, P., 2011. “Adding kinematic constraints to purely differential dynamics”. Computational Mechanics, 47(2), pp. 187–203.
[22] Garc´ıa de Jalo´n, J., and Bayo, E., 1994. Kinematic and Dynamic Simulation of Multibody Systems. The Real– Time Challenge. Springer–Verlag.
[23] Kamman, J. W., and Huston, R. L., 1984. “Dynamics of constrained multibody systems”. Journal of Applied Mechanics, 51(4), December, pp. 899–903.
[24] Singh, R. P., and Likins, P. W., 1985. “Singular value decomposition for constrained dynamical systems”. Journal of Applied Mechanics, 52(4), pp. 943– 948.
[25] Mani, N. K., Haug, E. J., and Atkinson, K. E., 1985. “Application of singular value decomposition for analysis of mechanical system dynamics”. Journal of Mechanisms, Transmissions, and Automation in Design, 107(1), pp. 82–87.
[26] Kim, S. S., and Vanderploeg, M. J., 1986. “QR decomposition for state space representation of constrained mechanical dynamic systems”. Journal of Mechanisms, Transmissions, and Automation in Design, 108(2), pp. 183–188.
[27] Liang, C. G., and Lance, G. M., 1987. “A differentiable null space method for constrained dynamic analysis”. Journal of Mechanisms, Transmissions, and Automation in Design, 109(3), pp. 405–411.
[28] Agrawal, O. P., and Saigal, S., 1989. “Dynamic analysis of multi-body systems using tangent coordinates”. Computers & Structures, 31(3), pp. 349–355.
[29] Dopico, D., Zhu, Y., Sandu, A., and Sandu, C., 2014. “Direct and adjoint sensitivity analysis of ordinary differential equation multibody formulations”. Journal of Computational and Nonlinear Dynamics, 10(1), p. paper 011012.
[30] Moler, C. B., and Stewart, G. W., 1973. “An algorithm for generalized matrix eigenvalue problems”. SIAM Journal on Numerical Analysis, 10(2), pp. 241–256.
[31] Gonza´lez, F., and K¨ovecses, J., 2013. “Use of penalty formulations in dynamic simulation and analysis of redundantly constrained multibody systems”. Multibody System Dynamics, 29(1), pp. 57–76.
[32] Shabana, A. A., 1998. Dynamics of Multibody Systems, 2nd ed. Cambridge University Press.
[33] Gonza´lez, F., Dopico, D., Pastorino, R., and Cuadrado, J., 2016. “Behaviour of augmented Lagrangian and Hamiltonian methods for multibody dynamics in the proximity of singular configurations”. Nonlinear Dynamics, 85(3), pp. 1491–1508.
[34] Lehoucq, B., Sorensen, D. C., and Vu, P., 1995. ARPACK: An implementation of the Implicitly Restarted Arnoldi Iteration that computes some of the eigenvalues and eigenvectors of a large sparse matrix. Tech. rep. Available from netlib $@$ ornl.gov under the directory scalapack.
[35] Flores, P., Machado, M., Seabra, E., and Tavares da Silva, M., 2011. “A parametric study on the Baumgarte stabilization method for forward dynamics of constrained multibody systems”. Journal of Computational and Nonlinear Dynamics, 6(1), p. paper 011019.
