Flexible-Link Multibody System Eigenvalue Analysis Parameterized with Respect to Rigid-Body Motion

: The dynamics of ﬂexible multibody systems (FMBSs) is governed by ordinary differential equations or differential-algebraic equations, depending on the modeling approach chosen. In both the cases, the resulting models are highly nonlinear. Thus, they are not directly suitable for the application of the modal analysis and the development of modal models, which are very useful for several advanced engineering techniques (e.g., motion planning, control, and stability analysis of ﬂexible multibody systems). To deﬁne and solve an eigenvalue problem for FMBSs, the system dynamics has to be linearized about a selected conﬁguration. However, as modal parameters vary nonlinearly with the system conﬁguration, they should be recomputed for each change of the operating point. This procedure is computationally demanding. Additionally, it does not provide any numerical or analytical correlation between the eigenpairs computed in the different operating points. This paper discusses a parametric modal analysis approach for FMBSs, which allows to derive an analytical polynomial expression for the eigenpairs as function of the system conﬁguration, by solving a single eigenvalue problem and using only matrix operations. The availability of a similar modal model, which explicitly depends on the system conﬁguration, can be very helpful for, e.g., model-based motion planning and control strategies towards to zero residual vibration employing the system modal characteristics. Moreover, it allows for an easy sensitivity analysis of modal characteristics to parameter uncertainties. After the theoretical development, the method is applied and validated on a ﬂexible multibody system, speciﬁcally using the Equivalent Rigid Link System dynamic formulation. Finally, numerical results are presented and discussed.

Modal analysis relies on the solution of an eigenvalue problem, which seeks the eigenvalues and eigenvectors associated to a linear system of equations. However, the adoption of such an analysis may provide a useful insight also for the study of mechanical systems whose dynamics is not governed by linear time-invariant equations. A significant class of mechanical systems that falls in this folder is the one of the flexible-link multibody systems (FMBSs). Such systems are robots or mechanisms that can deflect due to external loads or internal body forces, whose motion is described by means of kineto-elastodynamic models, hereafter referred to as dynamic models [17]. Several contributions can be found in the literature on the modeling of such systems as well as survey papers [18,19] and books [20,21]. The system elastic behavior is represented by continuous ordinary and partial differential equations. Such equations, so as to be simplified and solved, are discretized by means of lumped parameters, assumed modes, or finite element methods [19]. The common approach for modeling FMBSs consists in the use of the finite element method to discretize the flexible links and to represent their elastic deformations and in superposing such deformations to a known rigid body motion. Based on the set of coordinates chosen to model the rigid body motion, i.e., a minimum set of independent coordinates (representing the system degrees of freedom) or a redundant set of coordinates (including both the independent and dependent coordinates), the equations of motion are formulated as a set of ordinary differential equations (ODEs) or a coupled set of differential and algebraic equations (DAEs) to be solved simultaneously [20], respectively. In both the aforementioned cases, due to the large displacements to which FMBSs are subjected, their dynamic models are highly nonlinear and depend on the system configuration. Therefore, the related mass, stiffness, and damping matrices are in general nonconstant, as well as the resulting modal parameters.
A common practice used in the literature to apply the modal analysis to FLMBSs consists of linearizing the nonlinear dynamic model about a selected configuration so that the modal parameters can be computed. In particular, if the system model is formulated by means of ODEs, the eigenvalue analysis can be straightforwardly applied to the linearized system equations [12,22]. Conversely, if the system model is formulated by means of DAEs, two different strategies are possible for computing the eingensolutions: (a) transform the motion equations from DAEs to ODEs, linearize the resulting model, and then compute the eigenpairs [5,23,24]; (b) perform a direct eigenanalysis, i.e., the eigenanalysis for the system of equations resulting from the direct linearization of the DAEs [25,26]. The main difference between the eigensolutions obtained from linearized ODEs or DAEs is that the linearized ODEs allow to obtain the exact problem spectrum, while the linearized DAEs are affected by the linearization method and introduce spurious eigenvalues in the spectrum [27]. Due to such an evidence, this paper will focus on FMBSs modeled by means of a minimum set of ODEs.
Although the adoption of linearized models, on the one hand, allows the computation of the eigenpairs of FMBSs, on the other hand, it does not take directly into account the variability of the modal parameters due to the system configuration change. This last point is typically addressed by discretizing a given system motion/task in a certain number of operative points about which the nonlinear model is linearized and a new eigenvalue problem solved [5,6]. Following such an approach, several eigenvalue problems are to be solved, which is typically a computational expensive operation; additionally, unless interpolation techniques are used, the eigenpairs computed at different operating points are not related among them.
An approach that could help in overcoming such an issue has been proposed by Wittmuess et al. [28]. In such a paper, a method to get a parametric representation of the eigenvalues and eigenvectors of an undamped second-order mechanical system, whose model is analytically known, has been presented. In [29], Wittmuess et al. extended the method to proportionally damped systems. In this approach, the parametric representation of the eigenpairs is inferred from an iterative Taylor series expansion of the eigenvalue problem associated to the linearized system matrices about a parametric operating point. A first extension and application of the approach proposed by Wittmuess et al. to FMBSs characterized by small deformation and negligible damping and velocity-dependent terms has been proposed by the authors in [30,31]. In particular, in [31], the authors performed a preliminary investigation on the method capability to approximate the modal content over a wide range of the FMBS parameters, including not only the rigid motion coordinates, but also the payload handled by a two-degree-of-freedom (dof) planar robot carrying out a pick-and-place trajectory.
As promising results have been obtained in these preliminary studies, this paper aims at providing a comprehensive dissertation on the extension of such a method to FMBSs characterized by small deformations and non negligible damping and velocity-dependent terms. Indeed, the knowledge of an analytical relationship between the system motion condition and the natural frequencies, damping factors and modal shapes of, at least, the main vibrational modes, can be fruitfully exploited for, e.g., the development and implementation of more efficient model-based motion planning and control strategies towards to vibration minimization or the set-up of optimization problems based on the system modal characteristics. Additionally, a similar analytical expression allows for an easy sensitivity analysis of the system model to parameter uncertainties.
The paper is set out as follows. Starting from a nonlinear dynamic model formulated by means of a minimum set of ODEs, Section 2 derives the linearized dynamic model of a FMBS about a dynamic equilibrium configuration and discusses the eigenvalue problem for systems having nonsymmetric matrices, as it is the case of the linearized models of FMBSs. Section 3 outlines the method to derive the polynomial expressions for the eigenpairs of a FMBS in its generalized coordinates. In Section 4, the effectiveness of the method is proved by applying it to a flexible planar robot following a predefined trajectory. Finally, Section 5 gives concluding remarks.

Motion Equations
The total motion of a FMBS undergoing large rigid body motion and small elastic deformations can be modeled as an elastic motion, due to the link flexibility superposed onto a rigid body motion. Different strategies exist to model both the elastic and the rigid motions. The widespread approach to model link flexibility is by means of finite element model. Such a model is then embedded into a floating frame [20], in a corotational frame [32], or in a moving reference configuration [33], which represents the rigid motion. Depending on the selected model strategy and the type of coordinates used (i.e., relative, absolute, natural, minimal, or nodal coordinates) the motion equations result in a nonlinear set of DAEs or ODEs. Despite the large use of the DAE formulation, in this paper dynamic models formulated by means of ODEs will be adopted, as we are interested in the system modal parameters. Indeed, models employing ODEs lead to the exact system spectrum, unlike DAE formulation that may lead to approximate eigenvalues or introduce spurious ones [27].
The motion equations of a FMBS represented by means of a minimum set of second-order ODEs expressed in terms of the system degrees of freedom (i.e., independent coordinates) take the following matrix form, M(q r )q + C(q r ,q r )q + K(q r )q = f(q r , u) where q = q T f q T r T is the vector of the independent coordinates, including the vector q f ∈ R n f of the elastic coordinates of the flexible links and the vector q r ∈ R n r of the rigid body variables. The number of rigid coordinates, n r , is equal to number of rigid body motions of the system, whereas the number of the elastic coordinates, n f , depends on the number of finite elements employed to discretize the flexible links; their sum gives the total number of the system dofs, n do f = n f + n r . In Equation (1), M is the symmetric, positive definite mass matrix; C is the matrix containing the damping, centrifugal, and Coriolis terms; K is the system stiffness matrix, which is semi-positive definite as rigid-body motion is allowed. According to the model formulation employed, the K matrix can depend or not on the rigid body coordinates q r . Indeed, e.g., if the floating frame of reference formulation [20] is employed, K is constant, whereas in the case of the Equivalent Rigid Link System (ERLS) formulation [34], it depends nonlinearly on q r . The term on the right-hand side of Equation (1) represents the vector of the gravity, friction, and generalized external forces, it depends on the system coordinates, q r , and on the external inputs, u.

Linearization of the Equations of Motion
Due to the strong coupling between the gross rigid body motion and the fine motion (vibration) of the flexible links, the system motion equations are highly nonlinear differential equations in the coordinates and velocities.
To compute the modal characteristics of the system, the motion equations are to be linearized. Let us rewrite Equation (1) as follows.
Equation (2) can be linearized about an operating point P 0 = (q 0 ,q 0 ,q 0 , u 0 ) as follows. where In Equations (4) and (5), the symbol ⊗ indicates the inner product between the partial derivative of a matrix A ∈ R n do f ×n do f with respect to vector s ∈ R n do f by vector b ∈ R n do f : Equations (4)-(7) represent the most general formulation of the linearized system matrices, i.e., when the operating point is chosen as a dynamic equilibrium state of the system. As it is typically the case, whenever the system model is linearized about a static equilibrium position, i.e., P 0 = (q 0 , 0, 0, u 0 ), some of the terms in the linearized stiffness K and damping C matrices vanish. Then, it holds that Finally, the system motion equations linearized about an operating point can be written as

Modal Analysis for Systems with Nonsymmetric Matrices
The linearized damping matrix C in Equation (11) is, in general, neither proportional to the mass and stiffness matrices nor symmetric (due to the velocity-dependent terms); therefore, the eigenvalue problem associated to the model (11) takes the form of a quadratic eigenvalue problem: where λ is a system eigenvalue, and φ ∈ R n do f and ψ ∈ R n do f are the corresponding right and left eigenvectors, respectively. In Equation (13), symbol * denotes the conjugate transpose, as λ, φ, and ψ may be complex valued. The quadratic eigenvalue problem in Equations (12) and (13) has 2n do f eigenvalues, which are symmetric with respect to the real axis of the complex plane, being the system matrices real [35]. This means that the system eigenvalues can be either real or complex, but in the last case, they occur in complex conjugate pairs as well as the corresponding eigenvectors: if λ 1 = λ r + iλ i is a system eigenvalue, then λ 2 = λ r − iλ i is a system eigenvalue too, and the corresponding eigenvectors are φ 1 = φ r + iφ i and φ 2 = φ r − iφ i , respectively. Although the quadratic eigenvalue problem in Equations (12) and (13) could be directly solved, it is usually transformed in a generalized eigenvalue problem by reducing the n do f -dimensional set of second-order homogeneous differential equations in Equation (11) to a 2n do f -dimensional first-order set of equations: where N ∈ R n do f ×n do f can be any nonsingular matrix; here it is considered N = M. The quadratic (Equations (12) and (13)) and the generalized (Equations (14) and (15)) eigenvalue problems have the same 2n do f eigenvalues, whereas the corresponding eigenvectors are correlated by the following relations, If all the system eigenvalues are distinct and the left and the right eigenvectors are normalized so that w * i Bz i = 1, the following biorthonormality relations hold,

Taylor Series for Multivariate Function
Let E be a matrix depending nonlinearly on a set of p parameters collected in the vector x = if the parameter dependency is analytically known, it is possible to approximate the multivariate function E(x) with a polynomial by means of a Taylor series truncated at order t: where The last term on the right-hand side of Equation (19) has been obtained by using a multi-index formulation. In particular, σ is the multi-index, i.e., a tuple of non-negative integers in a number equal to the one of system parameters: σ = (σ 1 , σ 2 . . . , σ p ) with σ i ∈ N + , i = 1, . . . , p. In Equation (19) and in the remaining of the paper, the following definitions for a multi-index are adopted [36]: (5) number of ordered combinations of p positive integers whose sum is equal to |σ|: if and only if σ 1 = ς 1 , σ 2 = ς 2 , . . . , σ p = ς p .

Eigensolution Expansion
The dynamic model of a FMBS depends nonlinearly on the values assumed by the rigid motion coordinates and their derivatives. Therefore, the solution of the eigenvalue problem obtained by means of a linearized model has only a local validity, i.e., it holds for the expansion point P = P 0 . The validity of such a solution can be extended by deriving a polynomial representation of the eigenpairs in the system rigid motion coordinates and velocities. To this end, the Taylor series for multivariate functions (see Section 3.1) and the eigenvalue problem formulation of Equation (12) are here exploited. Let x ∈ R p≥2n r be the vector of model parameters including the rigid motion coordinates q r ; their velocitiesq r ; and, possibly, the parameters which the system model can depend on, such as payload mass, stiffness varying with the operative conditions, design parameters, and so on. Let us approximate, by means of Taylor series (see Equation (19)), the linearized system matrices, which are functions of x, about the operative point x 0 : Now, let us assume that also the eigenpairs can be approximated by means of a Taylor expansion: Note that having assumed to analytically know the system matrices, the Taylor series in Equation (20) are known expressions; conversely, the derivatives of the eigensolutions (λ with respect to x in Equation (21) are unknown quantities and must be determined. To this end, let us rewrite the eigenvalue problem of Equation (12) as a Taylor series expansion truncated to the t-th order through Equations (20) and (21): where α = (α 1 , α 2 , α 3 , α 4 ), β = (β 1 , β 2 , β 3 ) and γ = (γ 1 , γ 2 ) are three multi-indices. The components of such multi-indices set the order of each factor of the three addends in Equation (22), so that their whole order is equal to |α| = |β| = |γ|. Note that the entries of α, β and γ are multi-indices too. They play the same role of the multi-index σ in Equation (19), i.e., they are the multi-indices related to the system parameters x, and include p non-negative integer entries: If the terms in Equation (22) are now ordered by powers of δx, by defining Equation (22) can be rearranged as It results in a sum of n p,τ terms for each value of τ (see Section 3.1) made by two coefficients each, i.e., the term between brackets and δx τ . Equation (25) has to hold for all δx ∈ R p , meaning that for a given |τ|, all the n p,τ terms between brackets resulting from Equation (25) have to be equal to zero: Equation (26) contains the eigenvalue λ (order) i and eigenvector φ (order) i derivatives from order 0 to t. Let us rearrange such an equation so that the highest eigenpair derivatives are outside the sums: By solving Equation (27) for the highest derivatives of the eigenpairs it holds that where S i is a rectangular matrix with n do f rows and n do f + 1 columns. Symbol " † " indicates the pseudo-inverse matrix. Note that on the right-hand side of Equation (28) only the eigenpair derivatives up to order t − 1 appear. Therefore, by proceeding iteratively from the first order (that requires the only knowledge of the eigenpairs computed at the expansion point λ i ) up to the desired Taylor order, it is possible to compute all the coefficients of the Taylor series expansion in Equation (22). In such a way, the eigenvectors and eigenvalues of a FMBS are expressed as polynomial functions of the system rigid-body motion and of other possible variables of interest.
Although the method has been derived for the right eigenvectors (φ), its extension to the left (ψ) ones is trivial: the same procedure has to be repeated by substituting the linearized system matrices with their transposes, i.e., the starting eigenvalue problem is the following one, Note that the derivatives of each eigenpair (λ i , φ i ) are independent from the others 2n do f − 1 pairs; therefore, the problem can be solved for just the eigenpairs of interest, as well as if the eigenpairs are all complex conjugate, it is enough to solve the problem for n do f eigenpairs. The order t at which to truncate the Taylor series must be chosen as a trade-off between accuracy and complexity of the approximated eigensolutions. Indeed, a higher truncation order better approximates the real solutions, but the number of addends of the polynomial functions representing the eigensolutions rapidly increases with it, compromising the computational benefits of the method. An evidence of this is provided in Table 1. It shows the number of addends of the approximated functions expanded up to a certain order t and parameterized on the rigid-body motion of systems having n r rigid dofs. The parameterization about dynamic equilibrium configurations is computationally more demanding with respect to the one made about static equilibrium configurations, as the first uses twice the number of parameters. Therefore, a computationally efficient approximation of the eigenpairs about dynamic equilibrium configurations can be obtained only for systems with few rigid dofs. Conversely, if, as it is common in practice, the eigenpairs are computed about static equilibrium configurations, they can be efficiently approximated with the proposed method also if they depend on several rigid dofs. Regardless the computational efficiency and hence the possibility to employ the approximated eigenpairs in real-time applications, the method provides a function that puts in correlation the modal parameters with the system configurations; these could be used in offline methods such as optimization techniques.

Results
The effectiveness of the method in approximating the variability of the eigenpairs of FMBS due to system configuration change has been validated by means of the test case shown in Figure 1. It is an open-chain, planar mechanism with two flexible links and two revolute joints, and therefore two rigid dofs. The system has been modeled by means of the ERLS formulation [34], which leads to a dynamic model consisting of a minimum set of nonlinear ODEs (the interested reader can find the analytical nonlinear model of the studied system in the Appendix A.1). The link flexibility has been modeled through finite elements, in particular, uniform two-node and six-dof beam elements (see Figure 1) have been employed with the properties shown in Table 2. The system dynamic model has 23 dofs, of which two are the rigid motion coordinates and 21 are elastic dofs. The two rigid motion coordinates have been chosen as the absolute angular position of the shoulder and the elbow joints, respectively. They have been denoted q 1 and q 2 in Figure 1.
To test the correctness of the method in approximating the eigenpairs under motion condition, the system is forced to follow a given trajectory. In particular, the end-effector (point E in Figure 1) moves along a horizontal linear path of 0.3 m in 0.5 s. The corresponding trajectory in the joint coordinates is shown in Figure 2.  To apply the method under investigation, the system dynamic model has been linearized about a dynamic equilibrium configuration: P 0 = (q 0 = q(T/2),q 0 =q(T/2),q 0 =q(T/2), u 0 = u(T/2)), where T is the time period of the trajectory, the numerical values of the dynamic equilibrium configuration are stated in Appendix A.2. Then, the linearized system has been parameterized on the positions and velocities of the rigid motion coordinates, i.e., x = q 1 q 2q1q2 T .

Rigid dofs Elastic dofs
Starting from the solution of the eigenvalue problem computed at the expansion point x = x 0 = q 1 (T/2) q 2 (T/2)q 1 (T/2)q 2 (T/2) T , the eigenpair derivatives up to the fourth Taylor series expansion have been computed by means of Equation (28). Finally, fourth-grade polynomials, representing the eigenvalues and eigenvectors as function of the angular positions and velocities of the shoulder and the elbow joints, have been obtained.       The approximated eigenpairs have been compared with the exact ones, which are computed linearizing the model at each sample time (∆t = 0.01 s) and solving the corresponding eigenvalue problem. In particular, for comparison, the following indices have been employed, • relative percentage error on the undamped natural frequency, f : Let λ = λ r + iλ i be a system eigenvalue, it holds: • relative percentage error on the damping factor, ξ : • modal assurance criterion (MAC): In Equations (30), (32) and (33), the over-set tilde denotes the approximated quantities. The target value for the comparison of the adopted indices is 1 for the MAC and 0 for f and ξ .
The comparison indices for the first three elastic pairs of complex conjugate eigenvalues and eigenvectors are shown in Figures 3-5. The same figures also show the trend of the undamped eigenfrequencies and of the damping ratio along the trajectory, to provide an evidence of the variability of the modal parameters due to system configuration change. All the three figures show a good agreement between the exact and the approximated modal parameters along almost the entire trajectory. The biggest discrepancies (but still bounded) occur on the tails of the trajectory (at the beginning and the end), where the expansion parameters differ significantly from the values they assumed at the expansion point. Overall, the results appear very satisfying for all the three comparison metrics adopted.
Finally, note that, although the system has been parameterized around a dynamic equilibrium configuration, the computational time has been drastically reduced. Indeed, the solution of the eigenvalue problem in the 51 operating points in which the trajectory has been discretized requires approximately 0.135 s with a Matlab implementation on a machine with one processor of the type Intel Core i5-6200U CPU at 2.4 GHz with 8 gigabyte of RAM. Conversely, the computation of the approximations of the analyzed eigenpairs in the same points takes on average 0.008 s, leading to a reduction of the computation time of 94%.

Conclusions
This paper has outlined the entire procedure to get an analytical polynomial expression of the eigenpairs of a flexible-link multibody system (FMBS) as function of its configuration, in terms of rigid-body motion coordinates and velocities. The method is suitable for FMBSs characterized by small deformations, whose dynamic model is analytically known and formulated by means of a minimum set of ordinary differential equations. Starting from a dynamic model linearized about a static or dynamic equilibrium configuration, the eigenvalue problem in such a configuration is computed. By iteratively expanding the eigenvalue problem in a Taylor series with respect to the rigid motion coordinates, their velocities, and, possibly, the parameters which the system model can depend on (e.g., payload mass), the eigenpair derivatives at the expansion points are inferred. These terms represent the coefficients of the polynomial functions describing the approximated eigenpairs.
The correctness of the method in approximating the system eigenpairs under motion condition (while reducing the computational time), has been numerically proved with satisfactory results by employing a two-dof flexible planar manipulator forced to follow a given trajectory. An experimental validation of the method both in static and dynamic conditions will be addressed in future works by means of an impact analysis and an operational modal analysis, respectively.
The availability of a polynomial function for the eigenpairs explicitly depending on the system configuration can result very helpful for model-based techniques that need an accurate knowledge of the system modal characteristics, such as motion planning and control strategies towards to zero residual vibration.
The method allows to reach a desired accuracy in the eigenpair approximation by selecting the proper expansion order of the Taylor series. However, the selection of a high expansion order could compromise the computational benefits, as the number of addends of the polynomial functions that approximate the eigenpairs rapidly increases with the expansion order. Further research will focus on the development of an index that allows a priori estimation of the proper order of the truncated Taylor series to reach a desired accuracy on the approximated eigensolutions while preserving the computational efficiency.

Funding:
The research leading to these results can be framed within the COVI project (TN200Y) Free University of Bozen-Bolzano.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A. Test Case Implementation Details
Appendix A. 1

. Nonlinear Dynamic Model
The nonlinear dynamic model of the studied planar flexible-link mechanism, obtained through the equivalent rigid-link mechanism formulation, takes the form of Equation (1). The coordinate vector contains the following 21 elastic coordinates and 2 rigid coordinates, whose physical meaning can be inferred from Figure 1.
The force vector (see Equation (A2)) includes the effects of the gravity forces and the two torques, u 1 and u 2 , acting on the shoulder and elbow joints, respectively.

Appendix A.2. Dynamic Equilibrium Configuration
For numerically validating the method, the system dynamic model has been linearized about the following operating point P 0 = (q 0 ,q 0 ,q 0 , u 0 ).