# Assessment of Linearization Approaches for Multibody Dynamics Formulations 线性化方法的多体动力学公式评估 Francisco Gonz´alez∗ Juan de la Cierva Post-Doctoral Fellow Laboratorio de Ingenierı´a Meca´nica Universidade da Coru˜na 15403 Ferrol, Spain Email: f.gonzalez@udc.es Pierangelo Masarati Professor Dipartimento di Scienze e Tecnologie Aerospaziali Politecnico di Milano 20156 Milano, Italy Email: pierangelo.masarati@polimi.it Javier Cuadrado Professor Laboratorio de Ingenierı´a Meca´nica Universidade da Coru˜na 15403 Ferrol, Spain Email: javicuad@cdf.udc.es Miguel A. Naya Associate Professor Laboratorio de Ingenierı´a Meca´nica Universidade da Coru˜na 15403 Ferrol, Spain Email: minaya@udc.es 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. 采用多体动力学方法推导机械系统动力学方程通常会导致一组高度非线性的微分代数方程 (DAEs)。虽然这种运动方程的形式适用于广泛的实际应用，但在某些情况下，需要访问线性化的系统动力学。当需要进行稳定性分析和模态分析时，以及为某些控制算法和状态估计器定义植物和系统模型时，都需要线性化的动力学表达式。 文献中可以找到许多多体动力学线性化方法。这些方法在处理运动方程的方式和呈现结果的方式上存在差异，而这些差异又取决于用于描述机械系统的广义坐标的选择。这种选择与处理系统运动学约束的方式密切相关。可以区分出三种主要方法，可用于对迄今为止发表的大多数线性化方法进行分类。在这项工作中，我们将展示每种方法在静态平衡系统的线性化中的特性，并通过研究两个具有代表性的例子加以说明。 # 1 Introduction 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. 多体系统 (MBS) 动力学在过去几十年里经历了显著增长，这得益于计算机体系结构和软件的进步。MBS 动力学最具吸引力的特征之一是研究人员和分析师在处理特定问题时可以选择的各种建模和公式化方法[1]。一般来说，系统可以用一组独立的广义坐标进行建模，这会导致一组常微分方程 (ODEs)，或者使用相关的依赖坐标，其运动由一个或多个运动学约束相关联。在后一种情况下，动力学方程用一组 DAEs 表示，必须使用适当的算法来处理，例如 [2,3]。MBS 研究社区不断提出和开发新的公式和方法，这使得对它们的性能进行基准测试和表征变得必要。这项工作的最终目标是提供指导方针，帮助分析师在处理特定问题时选择最合适的公式。 对于那些需要线性近似运动方程的应用，也需要选择线性化方法。这适用于某些经典的控制算法，它们需要用线性模型来表示要控制的系统[4]，以及一些状态和输入估计器，如卡尔曼滤波器[5]。即使不需要运动学方程的线性表达式，它的可用性可以使人们更容易了解系统的行为。使用从原始微分代数方程的线性化中提取的降阶模型，可以简化结构的和气弹的问题分析[6]；线性化模型在稳定性分析中也很有用[7–9]。模态分析和确定系统的固有频率和振动模态是线性化动力学的另一个重要应用[10]。 在多体文献中，已经提出了几种方法来获得动力学方程的线性形式。它们具有不同的特征，主要取决于描述机械系统的坐标的选择。其中一些是用于机器人系统的递归算法[11]，另一些是基于 Maggi-Kane 坐标[12]。拉格朗日乘子是表示约束系统的常用方法，并且存在线性化所得方程的方法[6, 10, 13]。本文将多体动力学的线性化方法分为三大类，这些分类由广义坐标的选择定义。通过对静态平衡机械系统的动力学进行线性化，展示了它们的特性，从而保证获得线性时不变 (LTI) 问题。将每类方法的代表性应用于测试问题的线性化，并从效率、易用性和准确性方面进行比较。此外，在 [14] 中，确定了两种类型的线性化问题：高度约束的系统，其自由度数量远低于运动学变量的数量，以及（通常是柔性机构）自由度和运动学变量数量相似的系统。研究了线性化方法的行为，并对每种类型的系统给出了一个示例。文章最后向读者呈现了关于为每个情况选择最合适技术的实用指南。 # 2 Coordinates Selection 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 一个多体系统可以用一组 $n$ 个广义坐标 $\mathbf{q}=\{q_{1},\dots,q_{n}\}^{\mathrm{~T~}}$ 来描述，它们之间由 $m$ 个运动约束 $\Phi(\mathbf{q},t)=\mathbf{0}$ 相关联。为了本论文的目的，假设这些约束是全纯的且线性无关；因此，$m\times n$ 约束函数 $\Phi$ 关于广义坐标 $\mathbf{q}$ 的雅可比矩阵 $\Phi_{\mathbf{q}}$ 可以假设具有满行秩 $m$ 。运动方程可以表示为一组非线性DAE系统，如下所示： $$ \begin{array}{c}{\mathbf{M}\ddot{\mathbf{q}}=\mathbf{f}+\mathbf{f}_{c}}\\ {\Phi\left(\mathbf{q},t\right)=\mathbf{0}}\end{array} $$ 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. 其中 $\mathbf{M}\left(\mathbf{q}\right)$ 是 $n\times n$ 质量矩阵，$\mathbf{f}$ 是广义施加力和速度相关力的数组，$\mathbf{f}_{c}$ 是由运动学约束引入的反应力。 方程 (1) 通常被转化为文献中已有的多体形式中的一种，这些形式可分为三大类 [14]： 1. 最小坐标系 (MCS)，该形式的坐标数量等于系统的自由度数量，即 $n-m$； 2. 冗余坐标系 (RCS)，由 $n$ 个坐标和一个包含 $m$ 个拉格朗日乘子的集合组成； 3. 无约束坐标系 (UCS)，在该系统中，使用原始的 $n$ 个无约束问题的坐标，并将约束嵌入到形式化中。 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. 速度投影技术，其中系统速度被投影到允许运动的子空间 [15]，属于 MCS 类别。基于图论 [16,17] 和传递矩阵 [18] 的多体算法也应包含在这个类别中。 如果采用拉格朗日方法，约束反作用力可以表示为约束雅可比矩阵的函数，即 $\mathbf{f}_{c}=-\Phi_{\mathbf{q}}^{\mathrm{T}}\pmb{\lambda}$ ，其中 $\lambda$ 包含 $m$ 个拉格朗日乘子。使用 Baumgarte 稳定化来消除约束漂移 [19] 将直接得到一种 RCS 方法。或者，消除约束会得到一种 UCS 方法。UCS 方法也可以通过惩罚公式 [20] 或 [21] 中的力投影方法获得。 如前所述，选择一种方法或另一种方法决定了线性化方法的输出、效率以及它们传递关于机械系统信息的方式。MCS 格式化结果是一个具有约束系统精确谱的线性问题；对该问题的特征分析可以直接得到原始系统的 $2\left(n-m\right)$ 个特征值。RCS 和 UCS 方法可能会产生虚假或近似特征值；另一方面，它们的使用可以为了效率或易用性而方便。在后续部分，将从有效性和效率方面比较每个类别的代表，并讨论其主要特征。两个主要因素，即区分虚假特征值的难易程度以及线性化方程的简单性，也将从计算的角度来考虑。 # 3 Linearization Methods 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]. 本节介绍三种线性化 Eqs. (1) 的方法。第一种是基于速度投影得到的 MCS 格式化 [22]。第二种直接线性化了 RCS 对动态的描述 [6]。最后一种是基于 UCS 格式化，通过用惩罚系统替代运动学约束得到 [20]。 # 3.1 MCS Formulation: Velocity Projection The system dynamics can be expressed in terms of a set of $n-m$ independent velocities $\dot{\mathbf{z}}$ using the transformation 系统动力学可以表示为一组 $n-m$ 个独立的速率 $\dot{\mathbf{z}}$，使用变换： $$ {\dot{\mathbf{q}}}=\mathbf{R}{\dot{\mathbf{z}}} $$ 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 其中 $\mathbf{R}\left(\mathbf{q}\right)$ 是一个 $n\times(n-m)$ 速度变换矩阵，它满足 $\mathbf{R}^{\mathrm{T}}\boldsymbol{\Phi}_{\mathbf{q}}^{\mathrm{T}}=\mathbf{0}$ 。矩阵 $\mathbf{R}$ 可以通过多种方法确定：例如，坐标选择（将 $\mathbf{q}$ 分解为独立变量和相关变量），零特征值定理 [23]，奇异值分解 (SVD) [24,25]，QR 分解 [26]，以及 Gram-Schmidt 正交化 [27, 28]。所有这些方法都通过对雅可比矩阵 $\Phi_{\mathbf{q}}$ 进行操作，来确定满足约束条件的的速度子空间 $\mathbf{R}$ 。应用公式 (2) 中的变换，系统动力学变为 $(n-m)$ 阶常微分方程组。 $$ \mathbf{H}_{1}=\mathbf{R}^{\mathrm{T}}\mathbf{M}\mathbf{R}\ddot{\mathbf{z}}-\mathbf{R}^{\mathrm{T}}\left(\mathbf{f}-\mathbf{M}\dot{\mathbf{R}}\dot{\mathbf{z}}\right)=\mathbf{0} $$ 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 如果输入表示为 $\mathbf{u}$，则系统可以围绕平衡构型 $\mathbf{z}_{0},\,\dot{\mathbf{z}}_{0},\,\ddot{\mathbf{z}}_{0},\,\mathbf{u}_{0}$ 进行线性化，如下所示： $$ \begin{array}{c}{{{\bf{H}}_{1}\left({{{\bf{z}}_{0}}+\delta{{\bf{z}},{{\dot{\bf{z}}}}_{0}}+\delta{{\dot{\bf{z}}},{{\dot{\bf{z}}}}_{0}}+\delta{{\ddot{\bf{z}}},{{\bf{u}}_{0}}}+\delta{\bf{u}}\right)\cong}}\\ {{{\frac{{\partial{\bf{H}}_{1}}}{{\partial{\bf{z}}}}}|_{0}\delta{\bf z}+\frac{{\partial{\bf{H}}_{1}}}{{\partial{\dot{\bf{z}}}}}\bigg\vert_{0}\delta{\dot{\bf z}}+\frac{{\partial{\bf{H}}_{1}}}{{\partial{\dot{\bf z}}}}\bigg\vert_{0}\,\delta{\dot{\bf z}}+\frac{{\partial{\bf{H}}_{1}}}{{\partial{\bf{u}}}}\bigg\vert_{0}\,\delta{\bf u}}\end{array} $$ The terms in Eq. (4) take the form $$ \begin{array}{r l}&{\frac{\partial\mathbf{H}_{1}}{\partial\mathbf{z}}=\mathbf{K}_{r}=\left(\frac{\partial\mathbf{R}^{\mathsf{T}}}{\partial\mathbf{z}}\mathbf{M}\mathbf{R}+\mathbf{R}^{\mathsf{T}}\frac{\partial\mathbf{M}}{\partial\mathbf{z}}\mathbf{R}+\mathbf{R}^{\mathsf{T}}\mathbf{M}\frac{\partial\mathbf{R}}{\partial\mathbf{z}}\right)\hat{\boldsymbol{z}}}\\ &{\quad\quad-\frac{\partial\mathbf{R}^{\mathsf{T}}}{\partial\mathbf{z}}\left(\mathbf{f}-\mathbf{M}\mathbf{R}\hat{\boldsymbol{z}}\right)-\mathbf{R}^{\mathsf{T}}\left(\frac{\partial\mathbf{f}}{\partial\mathbf{z}}-\frac{\partial\mathbf{M}}{\partial\mathbf{z}}\mathbf{R}\hat{\boldsymbol{z}}-\mathbf{M}\frac{\partial\mathbf{R}\hat{\boldsymbol{z}}}{\partial\mathbf{z}}\right)}\\ &{\frac{\partial\mathbf{H}_{1}}{\partial\hat{\boldsymbol{z}}}=\mathbf{C}_{r}=-\mathbf{R}^{\mathsf{T}}\left(\frac{\partial\mathbf{f}}{\partial\hat{\boldsymbol{z}}}-\mathbf{M}\frac{\partial\mathbf{R}\hat{\boldsymbol{z}}}{\partial\hat{\boldsymbol{z}}}\right)}\\ &{\frac{\partial\mathbf{H}_{1}}{\partial\hat{\boldsymbol{z}}}=\mathbf{M}_{r}=\mathbf{R}^{\mathsf{T}}\mathbf{M}\mathbf{R}}\\ &{\frac{\partial\mathbf{H}_{1}}{\partial\mathbf{u}}=-\mathbf{F}_{r}=-\mathbf{R}^{\mathsf{T}}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}}\end{array} $$ 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 数值评估偏导数 $\mathbf{R}$ 的表达式详见 [29]。线性化动力学方程可以写为： $$ \mathbf{M}_{r}\delta{\ddot{\mathbf{z}}}+\mathbf{C}_{r}\delta{\dot{\mathbf{z}}}+\mathbf{K}_{r}\delta\mathbf{z}=\mathbf{F}_{r}\delta\mathbf{u} $$ 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. 方程 (6) 具有 $2\left(n-m\right)$ 个特征值，它们代表问题的精确谱。 此外，其首导矩阵 $\mathbf{M}_{r}$ 是对称正定矩阵。 因此，方程 (6) 可用于带有显式积分方案的逆动力学问题的求解。 # 3.2 RCS Formulation: Generalized Eigenanalysis广义特征分析 Following a Lagrangian approach, the dynamics equations (1) can be expressed in the form 采用拉格朗日方法，动力学方程 (1) 可以表达为： $$ \mathbf{H}_{2}=\left\{\begin{array}{c}{\mathbf{M}\ddot{\mathbf{q}}+\boldsymbol{\Phi}_{\mathbf{q}}^{\mathrm{T}}\boldsymbol{\lambda}-\mathbf{f}}\\ {\boldsymbol{\Phi}\left(\mathbf{q},t\right)}\end{array}\right\}=\mathbf{0} $$ which can be linearized directly and cast in descriptor form 可以被直接线性化并转化为描述符形式。 $$ \mathbf{E}_{q}\delta\dot{\mathbf{y}}=\mathbf{A}_{q}\delta\mathbf{y}+\mathbf{B}_{q}\delta\mathbf{u} $$ $$ \mathbf{E}_{q}=\left[\begin{array}{c c}{\mathbf{I}}&{\mathbf{0}}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{M}_{q}}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}\end{array}\right];\qquad\mathbf{A}_{q}=\left[\begin{array}{c c}{\mathbf{0}}&{\mathbf{I}}&{\mathbf{0}}\\ {-\mathbf{K}_{q}-\mathbf{C}_{q}\mathbf{\Lambda}-\mathbf{\Phi}\mathbf{\Phi}\mathbf{\Phi}\mathbf{0}}\\ {-\mathbf{\Phi}\mathbf{\Phi}\mathbf{q}}&{\mathbf{0}}&{\mathbf{0}}\end{array}\right] $$ with $$ {\bf M}_{q}={\bf M}\,;\quad{\bf C}_{q}=-\frac{\partial{\bf f}}{\partial\dot{\bf q}}\,;\quad{\bf K}_{q}=\frac{\partial{\bf M}}{\partial{\bf q}}\ddot{\bf q}+\frac{\partial\Phi_{\bf q}^{\mathrm{T}}}{\partial{\bf q}}\lambda-\frac{\partial{\bf f}}{\partial{\bf q}} $$ 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]. 需要注意的是，评估 $\mathbf{K}_{q}$ 需要计算在平衡状态下的拉格朗日乘子 $\lambda$。 除了受约束问题谱中的 $2\left(n-m\right)$ 个特征值外，对 $(2n+m)$ 个线性化方程 (8) 的广义特征分析会引入 $3m$ 个虚假特征值，它们与受约束的运动学变量和拉格朗日乘子相关。 这些虚假特征值很容易识别，因为它们的值要么是无穷大 [6]，要么是零，如果 $\Phi_{\mathbf{q}}$ 和 $\Phi_{\mathbf{q}}^{\mathrm{T}}$ 中的一个或两个从 $\mathbf{A}_{q}$ 传递到 $\mathbf{E}_{q}$。 由于方程 (7) 代表一个微分代数问题，矩阵 $\mathbf{E}_{q}$ 具有结构奇异性。 因此，方程 (8) 不能用于带有显式数值积分器的正向动力学问题。 然而，可以通过对矩阵 $\mathbf{E}_{q}$ 和 $\mathbf{A}_{q}$ 进行 QZ 分解来克服这个问题 [30]。 该方法的细节见 [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 $$ \pmb{\lambda}=\Xi\left(\ddot{\Phi}+\Theta\dot{\Phi}+\Omega\Phi\right) $$ 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 $$ \mathbf{H}_{3}=\mathbf{M}\ddot{\mathbf{q}}+\boldsymbol{\Phi}_{\mathbf{q}}^{\mathrm{T}}\boldsymbol{\Xi}\left(\ddot{\boldsymbol{\Phi}}+\boldsymbol{\Theta}\dot{\boldsymbol{\Phi}}+\boldsymbol{\Omega}\boldsymbol{\Phi}\right)-\mathbf{f}=\mathbf{0} $$ 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 $$ \begin{array}{l}{{\displaystyle{\bf H}_{3}\left({\bf q}_{0}+\delta{\bf q},{\dot{\bf q}}_{0}+\delta{\dot{\bf q}},{\ddot{\bf q}}_{0}+\delta{\dot{\bf q}},{\bf u}_{0}+\delta{\bf u}\right)\cong}}\\ {{\displaystyle\left.\frac{\partial{\bf H}_{3}}{\partial{\bf q}}\right\rvert_{0}\delta{\bf q}+\left.\frac{\partial{\bf H}_{3}}{\partial{\dot{\bf q}}}\right\rvert_{0}\delta{\dot{\bf q}}+\left.\frac{\partial{\bf H}_{3}}{\partial{\dot{\bf q}}}\right\rvert_{0}\delta{\ddot{\bf q}}+\left.\frac{\partial{\bf H}_{3}}{\partial{\bf u}}\right\rvert_{0}\delta{\bf u}}}\end{array} $$ or in the more compact form $$ \mathbf{M}_{p}\delta\ddot{\mathbf{q}}+\mathbf{C}_{p}\delta\dot{\mathbf{q}}+\mathbf{K}_{p}\delta\mathbf{q}=\mathbf{F}_{p}\delta\mathbf{u} $$ using the terms $$ \begin{array}{l}{\displaystyle\frac{{\partial{\bf{H}}_{3}}}{{\partial{\bf{q}}}}={\bf{K}}_{p}=\frac{{\partial{\bf{M}}}}{{\partial{\bf{q}}}}{\dot{\bf{q}}}+\frac{{\partial\Phi_{\bf{q}}^{\mathrm{T}}}}{{\partial{\bf{q}}}}\Xi\left({\ddot{\Phi}}+\Theta{\dot{\Phi}}+\Omega{\Phi}\right)\quad}\\ {\displaystyle\quad\quad+\ \Phi_{\bf{q}}^{\mathrm{T}}\Xi\left({\Omega{\Phi_{\bf{q}}}+\frac{{\partial{\Phi_{\bf{q}}}}}{{\partial{\bf{q}}}}\left({\ddot{\bf{q}}}+{\bf{\Theta}}{\dot{\bf{q}}}\right)+\frac{{\partial{\dot{\Phi}}_{\bf{q}}}}{{\partial{\bf{q}}}}{\dot{\bf{q}}}\right)\quad}\\ {\displaystyle\quad\quad+\ \Phi_{\bf{q}}^{\mathrm{T}}\Xi\left({\Theta\frac{{\partial{\Phi_{t}}}}{{\partial{\bf{q}}}}+\frac{{\partial{\dot{\Phi}}_{t}}}{{\partial{\bf{q}}}}}\right)-\frac{{\partial{\bf{f}}}}{{\partial{\bf{q}}}}}\end{array} $$ $$ \begin{array}{r l}&{\frac{\partial{\bf H}_{3}}{\partial{\bf\dot{q}}}={\bf C}_{p}}\\ &{\quad\quad=\Phi_{{\bf q}}^{\mathrm{T}}\Xi\left(\frac{\partial{\bf\dot{\Phi}}_{{\bf q}}}{\partial{\bf\dot{q}}}{\bf\dot{q}}+\frac{\partial{\bf\dot{\Phi}}_{t}}{\partial{\bf\dot{q}}}+\dot{\bf\Phi}_{{\bf q}}+\Theta{\bf\Phi}_{{\bf q}}\right)-\frac{\partial{\bf f}}{\partial{\bf q}}}\\ &{\quad\quad\frac{\partial{\bf H}_{3}}{\partial{\bf\dot{q}}}={\bf M}_{p}={\bf M}+\Phi_{{\bf q}}^{\mathrm{T}}\Xi\Phi_{{\bf q}}}\\ &{\quad\quad\frac{\partial{\bf H}_{3}}{\partial{\bf u}}=-{\bf F}_{p}=-\frac{\partial{\bf f}}{\partial{\bf u}}}\end{array} $$ Equation (14) is a system of $n$ ODEs. The method delivers an approximation of the $2\left(n-m\right)$ true system eigenvalues, ![](images/4dbaa92b67aa376a9c626b28d98d1f4d1a3a9ee5c377f60bc1a0d6a3d7c22fac.jpg) 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$ . 第3节中描述的方法被应用于机械系统在静态平衡构型附近的线性化。选择了两个例子：第一个是一个带有弹簧元件的 $N$ 环四杆机构；第二个是一个柔性双摆。四杆机构受到大量约束，只有一个自由度；它代表了在运动学约束的数量与广义坐标数 $\left(n-m\ll n\right)$ 相似的机械系统。另一方面，双摆包括与广义坐标 $\mathbf{q}$ 之间的柔性相关的自由度。在这种情况下，$n-m\approx n$。 # 4.1 Multiple Loop Four-Bar Linkage 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)$. 图1所示为由长度为$l_{f}=1\;\mathrm{m}$且质量均匀分布为$m_{f}=1~\mathrm{kg}$的等长杆组成的$N$环四杆机构。其在重力作用下运动，重力加速度为$g=9.81$ $\mathrm{m}/\mathrm{s}^{2}$。机构中每个环 $i$ 都有一个刚度为$k_{f}\,{=}\,25\,\mathrm{N/m}$且自然长度为$l_{f0}=\sqrt{2}\,\mathrm{m}$的弹簧连接点 $B_{i}$ 和 $A_{i-1}$。当 $N_{\mathrm{}}=1$ 时，该系统等效于[14]中讨论的测试案例。 该机构以点 $B_{0}–B_{N}$ 的 $x$ 和 $y$ 坐标以及连接点 $A_{N}$ 和 $B_{N}$ 的杆与全局 $x_{\mathrm{{}}}$ 轴之间的角度 $\boldsymbol{\upvarphi}$ 作为变量进行建模（共${n=2}N+3$个变量）。运动学约束项 $\Phi$ 由方程组成，这些方程强制每个杆的末端之间的距离在运动过程中保持恒定，外加一个方程，将 $\boldsymbol{\upvarphi}$ 与 $x_{N}$ 和 $y_{N}$ 相关联（$m=2N+2$）。 # 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 the velocity projection become $$ \begin{array}{l l}{{\displaystyle{\bf K}_{r}=-\frac{\partial{\bf R}^{\mathrm{T}}}{\partial{\bf z}}{\bf f}-{\bf R}^{\mathrm{T}}\frac{\partial{\bf f}}{\partial{\bf z}}\;;\;\;\;\;\;\;}}&{{\displaystyle{\bf C}_{r}=-{\bf R}^{\mathrm{T}}\frac{\partial{\bf f}}{\partial\dot{\bf z}}}}\\ {~~}&{{\displaystyle{\bf M}_{r}={\bf R}^{\mathrm{T}}{\bf M}{\bf R}\;;\;\;\;}}&{{\displaystyle{\bf F}_{r}={\bf R}^{\mathrm{T}}\frac{\partial{\bf f}}{\partial{\bf u}}}}\end{array} $$ ![](images/ce1c48625d6a75b393e72dafab788aea716830d364e810bc9e34d64ce2979baa.jpg) Fig. 2: A flexible double pendulum The terms for the linearization of the Lagrangian formulation in Section 3.2 are $$ \mathbf{M}_{q}=\mathbf{M}\,;\quad\mathbf{C}_{q}=-{\frac{\partial\mathbf{f}}{\partial{\dot{\mathbf{q}}}}}\,;\quad\mathbf{K}_{q}={\frac{\partial\Phi_{\mathbf{q}}^{\mathrm{T}}}{\partial\mathbf{q}}}{\boldsymbol{\lambda}}-{\frac{\partial\mathbf{f}}{\partial\mathbf{q}}} $$ ![](images/7399fa4b679180befc20e9c026e0f78ad0a0b1b0185c1509a99665b835fc1eec.jpg) 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 $$ \begin{array}{l}{{{\bf{K}}_{p}=\displaystyle\frac{\partial\hat{\bf{\Phi}}_{\bf{q}}^{\mathrm{T}}}{\partial{\bf q}}\Xi\Omega\Phi+\Phi_{\bf{q}}^{\mathrm{T}}\Xi\Omega\Phi_{\bf{q}}-\frac{\partial\mathrm{f}}{\partial{\bf q}}}}\\ {~~}\\ {{{\bf{C}}_{p}=\displaystyle\Phi_{\bf{q}}^{\mathrm{T}}\Xi\Theta\Phi_{\bf{q}}-\frac{\partial\mathrm{f}}{\partial\hat{\bf{q}}}}}\\ {~~}\\ {{{\bf{M}}_{p}={\bf{M}}+\Phi_{\bf{q}}^{\mathrm{T}}\Xi\Phi_{\bf{q}}~;~~~~~~~~~~~~~~~{\bf{F}}_{p}=\frac{\partial\mathrm{f}}{\partial{\bf u}}}}\end{array} $$ 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$
elapsed (ms)
NndS1MCSRCSUCS
152.2343±2.1477i0.1150.1080.185
5131.8922±1.5455i0.1880.3560.548
10231.8454±1.4352i0.3130.9541.684
15331.8230±1.3955i0.5561.8723.612
20431.8217±1.3750i0.7524.5176.152
Table 3: First eigenvalues of the flexible double pendulum discretized with $n_{p}=50$
S1S2S3S4
±1.3357i±3.5244i±9.4175i±22.4054i
Table 2: Eigenvalues of the damped $N_{}$ -loop four-bar linkage with $c_{f}=1\:\mathrm{Ns/m}$
N S1
1 -0.2424±2.1340i
5 -0.2350±1.5276i
10 -0.2325±1.4162i
15 -0.2316±1.3761i
20 -0.2312±1.3554i
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. ![](images/e31eba505a8a9f8ed7b341ad196ce6cbf1dde26927ccc4a5a81b942b376560bc.jpg) 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)
np10255075100
MCS0.0150.1851.3386.41117.940
RCS0.0200.1761.2295.80016.072
UCS0.0190.1701.2195.38814.228
Table 5: Flexible double pendulum: maximum differences with respect to the eigenvalues evaluated with the MCS method
npmax |Re (si)-Re(s*)Il RCS UCSmax Im (si) - Im (s*)Il RCS UCS
107.6.10-11 1.4·10-46.7.10-11 3.6·10-5
254.0·10-9 1.4·10-43.9·10-9 3.6·10-5
505.6·10-8 1.4·10-48.2·10-9 3.6·10-5
759.6·10-8 1.4·10-49.9·10-7 3.7·10-5
1001.5.10-7 1.4·10-42.5·10-6 3.8·10-5
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. Table 6: Summary of method properties
FeatureMCSRCSUCS
Efficiency (n ≈ m)
Efficiency (n 》m)
Forward dynamics with explicit integrationYesNoYes
Exact eigenvaluesYesYesNo
Spurious eigenvaluesNoYesYes
Handles singular configurationsNoYesYes
Degrades sparsityYesNoModerately
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. ![](images/b7c62f9b940c0204f1d552509ca0ed24d79ea86f8bba56142f4ebe17201cb819.jpg) 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.