## 2 Method Description  
### Modal Approximation Model  
The motion of an undamped FLMS can be decomposed into the sum of the large rigid-body motion of an equivalent system having rigid links, called ERLS, and the small elastic deformation of the flexible links with respect to the ERLS. The ERLS generalized coordinates are denoted by \(\mathbf{q}\in\mathbb{R}^{n_{q}}\), while the elastic deformation, modeled through finite element methods, is denoted by \(\mathbf{u}\in\mathbb{R}^{n_{u}}\). Since the rigid and elastic motion are intertwined, the dynamics of the FLMS is represented by the coupled ordinary differential equations  
\[\underbrace{\left[ \begin{array}{cc}\mathbf{M}_{e}&\mathbf{M}_{ e}\mathbf{S}\\ \mathbf{S}^{T}\mathbf{M}_{e}&\mathbf{S}^{T}\mathbf{M}_{e}\mathbf{S}\end{array}  \right]}_{\mathbf{M}}\left\{ \begin{array}{cc}\ddot{\mathbf{u}}\\ \ddot{\mathbf{q}}\end{array} \right\}+\underbrace{\left[ \begin{array}{cc}2 \mathbf{M}_{c_{G}}&\mathbf{M}_{e}\dot{\mathbf{S}}\\ 2\mathbf{S}^{T}\mathbf{M}_{c_{G}}&\mathbf{S}^{T}\mathbf{M}_{e}\mathbf{S}\end{array}  \right]}_{\mathbf{C}}\left\{ \begin{array}{cc}\dot{\mathbf{u}}\\ \dot{\mathbf{q}}\end{array} \right\}+\\ \underbrace{\left[ \begin{array}{cc}\mathbf{K}_{e}&0\\ 0&0\end{array} \right]}_{\mathbf{K}}\left\{ \begin{array}{cc}\mathbf{u}\\ \mathbf{q}\end{array} \right\}=\left[ \begin{array}{cc}\mathbf{M}_{e}&\mathbf{ I}\\ \mathbf{S}^{T}\mathbf{M}_{e}&\mathbf{S}^{T}\end{array} \right]\left\{ \begin{array}{ cc}\mathbf{g}\\ \mathbf{f}\end{array} \right\}. \tag{1}\]  
Let \(\mathbf{y}\in\mathbb{R}^{n_{y}}\) be the vector of design parameters of the system to assign the desired eigenfrequencies, which can include, e.g., geometrical and material properties that have consequence on inertial and stiffness terms. Therefore, the matrices and vectors in (1) depend nonlinearly not only on the ERLS coordinates, but also on the design variables. The system matrices here are defined as follows:  
* \(\mathbf{M}_{e}=\mathbf{M}_{e}(\mathbf{q},\mathbf{y})\) and \(\mathbf{K}_{e}=\mathbf{K}_{e}(\mathbf{q},\mathbf{y})\) are the mass and stiffness matrices of the assembled finite element matrices of all the flexible links, respectively;
* \(\mathbf{M}_{c_{G}}=\mathbf{M}_{c_{G}}(\mathbf{q},\mathbf{q},\mathbf{y})\) is the matrix of centrifugal and Coriolis terms;
* \(\mathbf{S}=\mathbf{S}(\mathbf{q},\mathbf{y})\) is the matrix which represents the relation between velocities at the nodes of the ERLS and the generalized velocities \(\dot{\mathbf{q}}\), and \(\dot{\mathbf{S}}=\dot{\mathbf{S}}(\mathbf{q},\dot{\mathbf{q}},\mathbf{y})\) is its time derivative;* \(\mathbf{g}=\mathbf{g}(\mathbf{q},\mathbf{y})\) is the gravity vector and \(\mathbf{f}=\mathbf{f}(\mathbf{q},\dot{\mathbf{q}},\mathbf{y})\) is the generalized external force vector. The square matrices \(\mathbf{M}=\mathbf{M}(\mathbf{q},\mathbf{y})\), \(\mathbf{C}=\mathbf{C}(\mathbf{q},\dot{\mathbf{q}},\mathbf{y},)\) and \(\mathbf{K}=\mathbf{K}(\mathbf{q},\mathbf{y})\) defined in (1) have order \(n_{\text{dof}}=n_{u}+n_{q}\). Such matrices can be used for numerical modal analysis for any given configuration of the system. Assuming that the motion of the ERLS is sufficiently slow, i.e., \(\dot{\mathbf{q}}\approx 0\), then all the velocity-dependent terms in Equation (1) can be neglected, hence matrix \(\mathbf{C}\). As a consequence, the eigenvalues and eigenvectors are solutions of the linear eigenvalue problem \[\mathbf{K}(\mathbf{q},\mathbf{y})\boldsymbol{\phi}_{i}(\mathbf{q},\mathbf{y}) =\lambda_{i}(\mathbf{q},\mathbf{y})\mathbf{M}(\mathbf{q},\mathbf{y}) \boldsymbol{\phi}_{i}(\mathbf{q},\mathbf{y}),\qquad i=1,\dots,n_{\text{dof}},\] (2) where \(\lambda_{i}\) is the \(i\)-th eigenvalue and \(\boldsymbol{\phi}_{i}\) is the associated eigenvector. The functions in the latter equation can be approximated in a neighborhood of \((\mathbf{q},\mathbf{y})=(\mathbf{q}_{0},\mathbf{y}_{0})\) by means of the Taylor polynomial of order \(p\), as \[