The Mdof methods are based on the assumption that the system can be represented as a multiple degree of freedom system. They are generally based on the construction of a high-order polynomial matrix whose coefficients are estimated from the measured FRFs. As for the SDOF methods, the estimation of the modal parameters is divided into two steps:
Natural frequencies and damping ratios extraction
The poles of the system are then obtained as the roots of the characteristic polynomial associated with this polynomial matrix. Three methods are implemented in this package:
- the Least-Squares Complex Exponential (LSCE) method,
- the Least-Squares Complex Frequency-domain (LSCF) method.
- the Polyreference Least-Squares Complex Frequency-domain (pLSCF) method,
LSCE method
The LSCE method is based on the assumption that the system can be modeled as a sum of complex exponentials in the time domain. The impulse response function (IRF) between an output \(o\) (\(o = 1, \ldots, n_o\)) and an input \(i\) (\(i = 1, \ldots, n_i\)) can be modeled as: \[
h_{oi}(t) = \sum_{k = 1}^n \left(c_k e^{\lambda_k t} + c_k^\ast e^{-\lambda_k^\ast t}\right) = 2 Re\left(\sum_{k = 1}^n c_k e^{\lambda_k t}\right).
\]
When sampling the IRF at a sampling period \(T_s\), one obtains the following discrete-time model: \[
h_{oi}[m] = 2 Re\left(\sum_{k = 1}^n c_k z_k^m\right) \text{ with } z_k = e^{\lambda_k T_s},
\] for \(m = 0, 1, \ldots, 2n\), (\(n\) is the model order). According to the Prony’s method, the poles \(\lambda_k\) can be calculated from the roots of the following characteristic polynomial: \[
\sum_{m = 0}^{2n} \beta_m z_k^m = 0 \text{ with } \beta_{2n} = 1.
\]
The coefficients \(\beta_l\) of the characteristic polynomial can be estimated from the sampled IRF \(h_{oi}[m]\). Todo so, one first multiply \(h_{oi}[m]\) by the corresponding coefficients \(\beta_m\) and then sums the resulting equations for \(m = 0, 1, \ldots, 2n\). This leads to the following homogeneous set of equations: \[
\sum_{m = 0}^{2n} \beta_m h_{oi}[m] = 2 Re\left(\sum_{k = 1}^{n} c_k \sum_{m = 0}^{2n} \beta_m z_k^m\right).
\]
When \(z_k\) is a root of the characteristic polynomial, the right-hand side of the previous equation is equal to zero. Therefore, one can estimate the coefficients \(\beta_m\) by solving the following set of linear equations: \[
\sum_{m = 0}^{2n} \beta_m h_{oi}[m] = 0.
\]
The previous equation can be derived using different set of data points, \(h_{oi}[m]\), \(h_{oi}[m+1]\), \(\ldots\), \(h_{oi}[m + 2n]\) to obtain an overdetermined set of linear equations that can be solved in the least-squares sense. In this package, the coefficients \(\beta_m\) are estimated by solving the following set of equations: \[
\begin{bmatrix}
h_{oi}[0] & h_{oi}[1] & \ldots & h_{oi}[2n-1] \\
h_{oi}[1] & h_{oi}[2] & \ldots & h_{oi}[2n] \\
\vdots & \vdots & \ddots & \vdots \\
h_{oi}[2n-1] & h_{oi}[2n] & \ldots & h_{oi}[4n-2]
\end{bmatrix}
\begin{bmatrix}
\beta_0 \\
\beta_1 \\
\vdots \\
\beta_{2n-1}
\end{bmatrix}
= -\begin{bmatrix}
h_{oi}[2n] \\
h_{oi}[2n+1] \\
\vdots \\
h_{oi}[4n-1]
\end{bmatrix}
\]
Once the coefficients \(\beta_m\) are estimated, the poles \(\lambda_k\) can be obtained from the roots \(z_k\) of the characteristic polynomial.
LSCF method
According to Guillaume et al. , the transfer function between an output \(o\) (\(o = 1, \ldots, n_o\)) and an input \(i\) (\(i = 1, \ldots, n_i\)) can be modeled in the frequency domain as: \[
\widehat{H}_k(\omega) = \frac{N_k(\omega)}{d(\omega)},
\] where \(k = 1, \ldots, n_o n_i\) (\(k = (o - 1) n_i + i\)). Here, \(N_k(\omega)\) is the \(k\)-th element of the numerator matrix defined as: \[
N_k(\omega) = \mathbf{\Omega}(\omega)^T \boldsymbol{\beta}_k,
\] while \(d(\omega)\) is the common denominator polynomial defined as: \[
d(\omega) = \mathbf{\Omega}(\omega)^T \boldsymbol{\alpha}.
\]
In the previous equations, \(\mathbf{\Omega}(\omega)\) is the vector of polynomial basis functions defined as: \[
\mathbf{\Omega}(\omega)^T = [1, \Omega_1(\omega), \ldots, \Omega_n(\omega)] \text{ with } \Omega_n(\omega) = \exp(-i\omega n T_s),
\] where \(n\) is the order of the denominator polynomial (a.k.a the model order) and \(T_s\) is the sampling period. The vectors \(\boldsymbol{\beta}_k\) and \(\boldsymbol{\alpha}\) contain the coefficients of the numerator and denominator polynomials, respectively.
Based on the previous formulation of the FRFs, the LSCF method consists in estimating the coefficients of the numerator \(\boldsymbol{\beta}_k\) and denominator \(\boldsymbol{\alpha}\) polynomials by minimizing the following functional: \[
J(\boldsymbol{\theta}) = \sum_{k=1}^{n_o n_i} \sum_{j=1}^{n_f} \left| W_k(\omega_j)\left[H_k(\omega_j) d(\omega_j) - N_k(\omega_j)\right] \right|^2,
\] where \(\boldsymbol{\theta} = [\boldsymbol{\beta}_1^T, \ldots, \boldsymbol{\beta}_{n_o n_i}^T, \boldsymbol{\alpha}^T]^T\) is the vector of unknown parameters to be estimated, \(H_k(\omega_j)\) is the measured FRF at frequency \(\omega_j\), \(n_f\) is the number of frequency lines used for the identification and \(W_k(\omega_j)\) is an arbitrary weighting function. In this package, a unitary weighting function is used, i.e. \(W_k(\omega_j) = 1\).
The poles of the transfer function are obtained from the roots of the characteristic polynomial associated with the denominator polynomial \(d(\omega)\), which are computed after having estimated the coefficients \(\boldsymbol{\alpha}\) as the argument of the minimum of the functional \(J(\boldsymbol{\theta})\) (see Ref. [2] for details).
pLSCF method
The pLSCF method is an extension of the LSCF method that makes it possible to take into account multiple reference sensors simultaneously (see Ref. [2] for details). When using only one reference sensor, the pLSCF method is equivalent to the LSCF method.
In this approach, the transfer function \(\widehat{\mathbf{H}}_o(\omega) \in \mathbb{C}^{1 \times n_i}\) between an output \(o\) (\(o = 1, \ldots, n_o\)) and a set of \(n_i\) inputs can be modeled in the frequency domain as: \[
\widehat{\mathbf{H}}_o(\omega) = \mathbf{N}_o(\omega) \mathbf{D}(\omega)^{-1},
\] where \(\mathbf{N}_o(\omega) \in \mathbb{C}^{1 \times n_i}\) is the numerator matrix defined as: \[
\mathbf{N}_o(\omega) = \mathbf{\Omega}(\omega)^T \boldsymbol{\beta}_o,
\] while \(\mathbf{D}(\omega) \in \mathbb{C}^{n_i \times n_i}\) is the denominator matrix defined as: \[
\mathbf{D}(\omega) = \mathbf{\Omega}(\omega)^T \boldsymbol{\alpha}.
\]
Based on the previous formulation of the FRFs, the pLSCF method consists in estimating the coefficients of the numerator \(\boldsymbol{\beta}_o\) and denominator \(\boldsymbol{\alpha}\) polynomials. To this end, one defines the weighted error between the measured and modeled FRFs as: \[
\boldsymbol{\varepsilon}_o(\omega_j, \boldsymbol{\theta}) = W_o(\omega_j)\left[\mathbf{H}_o(\omega_j) \mathbf{D}(\omega_j) - \mathbf{N}_o(\omega_j)\right]
\]
Then, an error matrix \(\mathbf{E}_o(\boldsymbol{\theta}) \in \mathbb{C}^{n_f \times n_i}\) is constructed by stacking the weighted errors \(\boldsymbol{\varepsilon}_o(\omega_j, \boldsymbol{\theta}) \in \mathbb{C}^{1 \times n_i}\) for all frequency lines, that is: \[
\mathbf{E}_o(\boldsymbol{\theta}) = \begin{bmatrix}
\boldsymbol{\varepsilon}_o(\omega_1, \boldsymbol{\theta}) \\
\vdots \\
\boldsymbol{\varepsilon}_o(\omega_{n_f}, \boldsymbol{\theta})
\end{bmatrix}
\]
Finally, the coefficients of the numerator and denominator polynomials are estimated by minimizing the following functional: \[
J(\boldsymbol{\theta}) = \sum_{o=1}^{n_o} \Vert\mathbf{E}_o(\boldsymbol{\theta})\Vert^2.
\]
The poles of the transfer function are obtained from the characteristic polynomial associated with the denominator matrix \(\mathbf{D}(\omega)\), by computing the eigenvalues of the so-called companion matrix.
Stabilization diagram
When carrying out an EMA, the main question that arises is how many modes (or, more precisely, poles) are contained within the frequency band of interest. Analyzing the frequency response functions (transfer functions) provides an initial estimate. However, it is difficult to distinguish between multiple or nearby modes.
The stabilisation diagram is a fundamental EMA tool as it provides a clear and concise presentation of the poles of the system under study. The idea is to adjust the transfer function model by increasing the number of poles and displaying the results on a graph. Broadly speaking, increasing the order of the model (i.e. the number of poles) causes the true vibration modes to converge to a stable value, whereas the numerical modes appear inconsistently or randomly on the graph. Thus, as the poles converge to stable ones, markers appear on the graph to indicate the identified poles, as the order of the model is gradually increased. The selection of poles is generally based on the stability of natural frequencies, modal damping and/or mode shapes, moving from a model order \(n\) to an order \(n + 1\). In this package, the stability of the natural frequencies and damping ratios is used to identify the physical poles of the system. The stability criterion is set to 1% for the natural frequencies and 5% for the damping ratios, following common practice in the EMA field.
In StructuralVibration.jl, the stabilization diagram can be obtained using the stabilization function along with any of the Mdof methods described above. The identified stable poles can then be visualized using the stabilization_plot function. The stabilization function returns an EMAMdofStabilization structure that contains all the information related to the stabilization analysis, which allows the user to manually select the appropriate poles and to further analyze (e.g. by using quality and analysis indicators) and visualize the results.
