Operational modal analysis using symbolic regression for a nonlinear vibration system

Abstract

The existence of nonlinearity is an inevitable frequent occurrence that should be considered to accurately identify the modal parameters of a vibration system using operational modal analysis. A problem is that the traditional operational modal analysis method based on the linear modal theory is not applicable to modal parameter identification of vibration systems with nonlinearity. A solution is as follows: this paper is aimed at solving the problem by proposing a new operational modal analysis method to carry out modal parameter identification for a nonlinear vibration system. The new operational modal analysis method, based on the forced response and symbolic regression method without assuming any pre-existing information and only using mathematical symbols, is introduced to solve the problem by automatically searching for the expression structure and modal parameters of a system in nonlinear normal modes. The simulation result of a three-degrees-of-freedom nonlinear system reveals the high accuracy of the proposed operational modal analysis method in extracting the modal parameters. Then, a rod fastening rotor model is considered, and the capability of the proposed operational modal analysis method to precisely extract its modal parameters is further evaluated.

Introduction

System identification is now aimed at the improvement of structural models built from input and output measurements for a real structure with vibration sensing devices. The modal parameters (namely, the natural frequencies, damping ratios, and mode shapes) are the basis of the model of the system to be expressed. In the past several decades, nonlinear normal modes (NNMs), which are a natural property of structural and mechanical systems, have been researched by many scholars. Moreover, NNMs are a good tool for performing modal analysis of nonlinear mechanical systems. The meaning of the NNM can be extended from the concept of the linear normal mode (LNM) for linear systems, but the NNM of structural and mechanical systems cannot be reflected in the decoupling of equations of motion. In other words, similar to the LNM for a linear system, the NNM for a nonlinear system is defined as any periodic motion of any undamped and autonomous system, but its principles and mechanisms are more complicated. A detailed description and explanation of NNMs can be found in the literature reviews.^1,2 To date, a few methods have been introduced to solve this problem.

In addition, as reported in Peeters et al., ³ no direct constructive method exists to determine the appropriate excitation of a given NNM. Thus, experimental modal analysis (EMA) is not suitable for NNM identification, but operational modal analysis (OMA) has great advantages in solving this problem.

The problem is as follows. Currently, OMA methods based on linear system theory are mainly suitable for linear vibration systems, and not when nonlinearity exists. For example, the frequency-domain modal identification method mainly extracts modal parameters by solving the frequency response function of the system. ⁴ In addition, modal parameters can also be identified using the response signal of the measured structure in the time domain. Ibrahim ⁵ identified the modal parameters of spacecraft using the stochastic reduction method combined with the Ibrahim time domain method. The random subspace method stochastic subspace identification ⁶ proposed by Peeters et al. in 1995 was also based on a linear system space model. Additionally, the blind source separation method^7,8 obtained modal parameters as unknown mixed channel parameters by obtaining observation signals through multiple sensors to estimate source signals.

For modal analysis of nonlinear vibration systems, Liao in 2014 proposed an approach to construct the relationship between the NNMs and the forced response of a nonlinear system. ⁹ Peeters et al. in 2011 proposed a nonlinear EMA method for modal testing of nonlinear vibrating structures based on NNMs. ¹⁰ The linearization approach was used to study nonlinear systems, and theories were developed in structural dynamics for the investigation of nonlinear systems.^11–13 He^14,15 and He and Sun ¹⁶ established a variational formulation in the thin film nonlinear vibration system, and the identification of the unknown function is given in detail.

There are no OMA methods for vibration system identification in the presence of nonlinearity at home or abroad. Therefore, we propose an OMA method to extract modal parameters using NNMs for vibration system identification in the presence of nonlinearity. Hence, based on forced response data, the symbolic regression method (SR) without assuming any pre-existing information and only using mathematical symbols is introduced to solve the problem by automatically searching for the expression structures and parameters of a nonlinear system in NNMs. Moreover, genetic programming has been successfully applied in the field of nonlinear function approximation, prediction, and modeling.

The paper is organized as follows: the upcoming section presents a theoretical analysis of identification methods to obtain the modal parameter of a system using NNMs. Next, a test of the proposed method is presented with a three-degrees-of-freedom nonlinear system. Subsequently, the application of the proposed method is presented to extract the NNMs of a rod fastening rotor model. Finally, a summary is provided.

Theoretical formulation

Nonlinear system model

All systems introduced in this report are assumed to be deterministic and time-invariant. For instance, under the condition of a given excitation, the system response is always the same without any uncertainty. The generalized multi-degrees-of-freedom nonlinear vibration system is shown in Figure 1, and its generalized differential equation is as follows (n degrees of freedom) ⁹

$$\mathbf{M}\ddot{\mathit{x}}+\mathit{C}\dot{\mathit{x}}+\mathit{Kx}+{\mathit{f}}_{\mathit{nl}}\left(\dot{\mathit{x}},\mathit{x},t\right)=\mathit{p}\left(t\right)$$

(1)

Figure 1.

Structure schematic of the nonlinear system.

Figure 2.

Calculation and verification procedure of the OMA method based on SR.

$\ddot{\mathit{x}}$ , $\dot{\mathit{x}}$ , and $\mathit{x}$ are, respectively, the acceleration, velocity, and displacement arrays of order n described by the physical coordinates; $\mathit{p}(t)$ is an external excitation array of order n; M and K are, respectively, the system mass matrix and stiffness matrix, which are n-order square arrays (real symmetric ones); C is the viscous damping matrix, which is a positive definite or semidefinite symmetric matrix and an n-order square array; and ${\mathit{f}}_{\mathit{nl}}\left(\dot{\mathit{x}},\mathit{x},t\right)$ is the nonlinear internal force. The nonlinear factor ${\mathit{f}}_{\mathit{nl}}(\dot{\mathit{x}},\mathit{x},t)$ will produce a self-excitation phenomenon under single-frequency excitation and induce nonlinear modes.

By means of mathematics, general forces and responses can undoubtedly be written as sine and cosine series, respectively. In other words, the response or excitation is expressed as a linear combination of contributions from resonant NNMs. ¹⁷

Nonlinear response $\mathit{x}\left(t\right)$

$$\mathit{x}\left(t\right)={\displaystyle \sum }_{k=1}^{\infty }{\mathit{X}}_k\text{cos}\left(k\omega t\right)$$

(2)

Excitation $\mathit{p}\left(t\right)$

$$\mathit{p}\left(t\right)={\displaystyle \sum }_{k=1}^{\infty }{\mathit{P}}_k\text{sin}\left(k\omega t\right)$$

(3)

A number of scholars have contributed to the development of the HBM (harmonic balance method).^18,19 The HBM expands the periodic response of a nonlinear system as a truncated Fourier series. The unknown time response $\mathit{x}\left(t\right)$ can be written through the practical application of the Fourier series of up to N_H

$$\mathit{x}\left(t\right)={\displaystyle \sum }_{k=0}^{N_H}{\mathit{U}}_k^c\cos \left(k\omega t\right)+{\mathit{U}}_k^s\sin \left(k\omega t\right)+\epsilon \left(t\right)$$

(4)where ${\mathit{U}}_k^c$ and ${\mathit{U}}_k^s$ are vectors of cosinusoidal and sinusoidal harmonic coefficients; $\epsilon (t)$ is the error, which can be as small as desired, and this is an optimization problem determined by the method selected.

For most events in mechanical engineering, the response is related to the positions of the mechanical structure, so the equations of motions (1) of the structure can be written more clearly as

$$\left[\begin{array}{cc}{\mathbf{M}}_{\mathit{UU}}& {\mathbf{M}}_{\mathit{VU}}\\ {\mathbf{M}}_{\mathit{UV}}& {\mathbf{M}}_{\mathit{VV}}\end{array}\right]\left[\begin{array}{c}{\ddot{\mathit{w}}}_{\mathit{U}}\\ {\ddot{\mathit{w}}}_{\mathit{V}}\end{array}\right]+\left[\begin{array}{cc}{\mathbf{C}}_{\mathit{UU}}& {\mathbf{C}}_{\mathit{VU}}\\ {\mathbf{C}}_{\mathit{UV}}& {\mathbf{C}}_{\mathit{VV}}\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{w}}}_{\mathit{U}}\\ {\dot{\mathbf{w}}}_{\mathit{V}}\end{array}\right]+\left[\begin{array}{cc}{\mathbf{K}}_{\mathit{UU}}& {\mathbf{K}}_{\mathit{VU}}\\ {\mathbf{K}}_{\mathit{UV}}& {\mathbf{K}}_{\mathit{VV}}\end{array}\right]\left[\begin{array}{c}{\mathbf{w}}_{\mathit{U}}\\ {\mathbf{w}}_{\mathit{V}}\end{array}\right]+\left[\begin{array}{c}{\mathit{f}}_{\mathit{nl},\mathit{U}}\left({\dot{\mathit{w}}}_{\mathit{U}},{\mathit{w}}_{\mathit{U}},{\dot{\mathit{w}}}_{\mathit{V}},{\mathit{w}}_{\mathit{V}},t\right)\\ {\mathit{f}}_{\mathit{nl},\mathit{V}}\left({\dot{\mathit{w}}}_{\mathit{U}},{\mathit{w}}_{\mathit{U}},{\dot{\mathit{w}}}_{\mathit{V}},{\mathit{w}}_{\mathit{V}},t\right)\end{array}\right]=\mathit{p}\left(t\right)$$

(5)where p (t) is the external excitation vector and U and V are two positions.

On the one hand, linear modal theory is not suitable for the modal analysis of vibration systems with nonlinearity, but the concept of the NNM helps us solve the problem well by extending linear modal theory to the generalized nonlinear mode.^1,17 On the other hand, the response is a binary function of time and position, and the internal interactions at different positions are coupled. After the response is subjected to the variable separation method, ²⁰ its expression can be updated to

$$\begin{array}{l}w\left(x,t\right)={\displaystyle \sum _{m=1}^M{\displaystyle \sum _{k=1}^N\left\{A_{m,k}\cos \left(k{\omega }_{k}t\right)+B_{m,k}\sin \left(k{\omega }_{k}t\right)\right\}\cdot C_{m,k}\sin \left({\lambda }_{m}x\right)}}\\ +\left\{A_{m,k}\cos \left(k{\omega }_{k}t\right)+B_{m,k}\sin \left(k{\omega }_{k}t\right)\right\}\cdot D_{m,k}\cos \left({\lambda }_{m}x\right)\end{array}$$

(6)where N and M are the number of generalized nonlinear normal modes and the number of positions, respectively; A_m,k and B_m,k are the coefficients in generalized mode coordinates of each order; and ${\omega }_k$ is the generalized nonlinear modal frequency.

Proposed OMA method

There are no OMA methods for vibration system identification in the presence of nonlinearity at home or abroad. Therefore, we propose an OMA method to extract modal parameters using NNMs for vibration system identification in the presence of nonlinearity. Hence, based on forced response data, the symbolic regression method without assuming any pre-existing information and only using mathematical symbols is introduced to solve the problem by automatically searching for the expression structures and parameters of a nonlinear system in NNMs. Moreover, genetic programming has been successfully applied in the field of nonlinear function approximation, prediction and modeling.

Symbolic regression has been conducted in extensive research by many domestic and foreign scholars. The natural laws implied in experimental data can be excavated with the usage of symbolic regression. In addition, the symbolic regression algorithms mainly include genetic programming (GP), gene expression programming, and stepwise regression algorithms for randomly generating candidate factor sets. Of the several methods, research on GP is the most widespread. The GP algorithm can identify the model’s expressions and parameters without assuming or relying on any pre-existing information and only using mathematical symbols before identification. Furthermore, GP does not assume the form of a function and automatically searches for the form of the mathematical formula and its parameters using the data and only mathematical symbols.^21–26 In particular, the process of GP is an adaptive nonlinear search process under the guidance of the fitness and a generalized hierarchical computer program to describe the problem. ²⁷

To demonstrate how the proposed nonlinear OMA method works, a simple NNM is taken as an example and shown in Figure 2. Its function expression $f\left(x\right)=\mathit{A}\cos \left(\omega t+\varphi \right)$ is selected for illustrating the process. The GP algorithm (pictured in Figure 3) was presented by Professor Koza of Stanford University in the early 1990s. Based on the idea of biological evolution, it involves a function search process in which genetic operators, such as copy, crossover, and mutation, are applied on a randomly generated population, as shown in Figure 3(a). The measure of the predictive ability in extracting the modal parameters is determined using the fitness function

$${\mathit{J}}_j={\displaystyle \sum }_{i=1}^m\left|\frac{{\mathit{Y}}_j\left(i\right)-{\widehat{\mathit{Y}}}_j\left(i\right)}{{\mathit{Y}}_j\left(i\right)}\right|$$

(7)where Y is the searched function expression vector for the vibration system model, Yj is the jth element of Y , and i is the ith value of the response time series. The purpose of the GP algorithm is to find a symbolic expression that best satisfies a given set of system response time series produced for an example NNM according to the fitness criterion. First, we must determine the function symbol set F and the variable symbol set T considering the NNMs used for operational modal analysis. The characteristics of the expression structure of the example NNM for nonlinear OMA determine the choice of function set: +, –, *, /, sin, cos and end points t and const. Then, the structure of the chromosome is constructed. For an individual representing an NNM in GP, a hierarchical tree shape is used (Figure 3(c)). Finally, the example NNM will be searched for by the GP procedure.

Figure 3.

(a) Schematic of the automatic evolution of GP for the OMA method in searching for the modal parameter expressions; (b) GP algorithm workflow; (c) hierarchical structure of the function expression $f\left(x\right)=A\text{cos}\left(\omega t+\varphi \right)$ .

Faced with the complicated expression structure in equation (4), we have to select the most suitable expression from all the expressions searched by GP, as shown in equation (8). Generally, the proposed nonlinear OMA method is carried out by obtaining the nonlinear modes and their modal parameters directly from the response time series. The method uses only mathematical symbols, including +, –, *, /, sin, cos and the end points t and const, without assuming or relying on any pre-existing information, especially the linear modal theory. In other words, we normally well consider and adapt to the presence of nonlinear components in the vibration system using GP.

$${\mathit{x}}_{\mathit{Symbol} \mathit{Regression}}\left(t\right)={\displaystyle \sum }_{k=0}^{N_{}}\left({\mathit{A}}_k^c\cos \left({\omega }_{k}t\right)+{\mathit{A}}_k^s\sin \left({\omega }_{k}t\right)\right)$$

(8)where N is the number of generalized nonlinear normal modes; ${\mathit{A}}_k^c$ and ${\mathit{A}}_k^s$ are the coefficients in generalized mode coordinates of each order; and ${\omega }_k$ is the generalized nonlinear modal frequency.

For the general nonlinear vibration system, when its forced response is obtained, GP without assuming or relying on any pre-existing information and only using mathematical symbols is used before identification. In addition, the NNM as a function expression and its coefficients can be automatically extracted by performing NNM identification with the advantages of high accuracy and relatively simple operation. Finally, GP can be well applied to OMA of nonlinear vibration.

A test model

To avoid the situation in which the selected test system is too weakly nonlinear to influence the general results, a strongly nonlinear system is needed. A strongly nonlinear system means a system for which nonlinear terms are of the same level as linear terms. Taking a three-degrees-of-freedom nonlinear vibration system as an example, the external force selected is harmonic. The single NNM and its natural frequency w_nl can be identified using the OMA based on GP under the condition that the system’s response time series are obtained. Finally, the performance of the proposed OMA method is compared with that of the theoretical reference estimation calculated by the Runge–Kutta method ²⁸ based on several criteria.

Model equations

The three vibrators in Figure 4 are connected to three nonlinear springs and dampers. We write out the equations for this nonlinear vibration system and introduce a transformation matrix p such that $\mathit{x}=\mathit{pv}$ , which are decoupled according to the linear part and can be generalized as ²⁹

$$\left[\begin{array}{c}{\ddot{u}}_1\\ {\ddot{u}}_2\\ {\ddot{u}}_3\end{array}\right]+\left[\begin{array}{l}{\mu }_{1}00\hfill \\ 0{\mu }_{2}0\hfill \\ 00{\mu }_3\hfill \end{array}\right]\left[\begin{array}{c}{\dot{u}}_1\\ {\dot{u}}_2\\ {\dot{u}}_3\end{array}\right]+\left[\begin{array}{ccc}\frac{k_1}{m_1}& 0& 0\\ 0& \frac{k_2}{m_2}& 0\\ 0& 0& \frac{k_3}{m_3}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]+\left[\begin{array}{c}{f}_1\left(u_1,{\dot{u}}_1,u_2,{\dot{u}}_2,u_3,{\dot{u}}_3\right)\\ f_2\left(u_1,{\dot{u}}_1,u_2,{\dot{u}}_2,u_3,{\dot{u}}_3\right)\\ f_3\left(u_1,{\dot{u}}_1,u_2,{\dot{u}}_2,u_3,{\dot{u}}_3\right)\end{array}\right]=\left[\begin{array}{c}{F}_1\text{cos}({\text{Ω}}_{1}t)\\ F_2\text{cos}({\text{Ω}}_{2}t)\\ F_3\text{cos}({\text{Ω}}_{3}t)\end{array}\right]$$

(9)where $u_1,u_2,\text{ and }{u}_3$ are the responses; $m_{1,2,3}$ and $k_{1,2,3}$ are the system mass and stiffness, respectively; ${\mu }_{1,2,3}$ are the viscous damping values; $f_{1,2,3}\left(u_1,{\dot{u}}_1,u_2,{\dot{u}}_2,u_3,{\dot{u}}_3\right)$ are the nonlinear internal forces, including all possible relationships between the parameters in the system, and $F_1,F_2,F_3,{\text{Ω}}_1,{\text{Ω}}_2,{\text{and Ω}}_3$ are the amplitudes and circular frequencies of external excitations.

Figure 4.

A three-DOF coupled nonlinear oscillator system.

Considering the object of this study, the system excitation equation used is as follows

$$\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]=\left[\begin{array}{c}{f}_{11}{u}_1^2{\dot{u}}_1+f_{12}{u}_1^3+f_{13}{u}_1^2{u}_2+f_{14}{u}_1^{}{u}_2^2+f_{15}{u}_2^3\\ f_{21}{u}_1^2{u}_2+f_{22}{u}_1^{}{u}_2^2+f_{23}{u}_2^3+f_{24}{u}_1^{}{u}_3^2+f_{25}{u}_2^2{\dot{u}}_2\\ f_{31}{u}_1^2{u}_2+f_{32}{u}_1^{}{u}_2^2+f_{33}{u}_2{u}_3^2+f_{34}{u}_3^3+f_{35}{u}_3^2{\dot{u}}_3\end{array}\right]$$

(10)

At the same time, the harmonic balance method ¹⁶ is adopted to solve the above equation to obtain the analytic solution equation (11). According to the expression structure of equation (11), the OMA method proposed in an earlier section is suitable for this work to perform OMA and extract the modal parameters.

$$\left[\begin{array}{l}{u}_1\\ u_2\\ u_3\end{array}\right]=\left[\begin{array}{c}{\displaystyle \sum _{i=0}^{\infty }{\text{α}}_{\text{i}}\cos \left({\text{τ}}_{\text{i}}t+{\phi }_i\right)-{\text{μ}}_1\cos \left({\text{Ω}}_{1}t\right)-{\text{δ}}_1\text{sin}\left({\text{Ω}}_{1}t\right)}\\ {\displaystyle \sum _{j=0}^{\infty }{\text{α}}_{\text{j}}\cos \left({\text{τ}}_{\text{j}}t+{\phi }_j\right)-{\text{μ}}_2\cos \left({\text{Ω}}_{2}t\right)-{\text{δ}}_2\text{sin}({\text{Ω}}_{2}t)}\\ {\displaystyle \sum _{k=0}^{\infty }{\text{α}}_{\text{k}}\cos \left({\text{τ}}_{\text{k}}t+{\phi }_k\right)-{\text{μ}}_3\cos \left({\text{Ω}}_{3}t\right)-{\text{δ}}_3\text{sin}({\text{Ω}}_{3}t)}\end{array}\right]$$

(11)

Numerical results

To extract the single NNM and further obtain the modal parameters as a theoretical reference, the forced responses of the above system must be calculated first. After choosing reasonable coefficients, ${\mu }_1=2$ , ${\mu }_2=8$ , ${\mu }_3=4, \frac{k_1}{m_1}=6, \frac{k_2}{m_2}=2, \frac{k_3}{m_3}=2, F_{1, 2, 3}=20$ , and ${\Omega }_{1,2,3}=20$ , with the external force given as

$$f_{\mathit{select}}=\left[\begin{array}{c}6u_1^2{\dot{u}}_1+2u_1^2{u}_2\\ 0.5u_1^2{u}_2+10u_2^3+u_1^{}{u}_3^2\\ u_2{u}_3^2+5u_3^2{\dot{u}}_3\end{array}\right]$$

(12)

The Runge–Kutta method ²⁸ is adopted for solving equation (9), and then, periodic solution expressions of responses $u_1$ , $u_2$ , and $u_3$ are acquired. The calculated solutions are taken as the theoretical references, and two NNMs are obtained, for which the natural frequencies ${\text{ω}}_{1,0}=0.2251$ and ${\text{ω}}_{2,0}=0.3898$ are obtained. Based on the response time series u_i, the single NNMs given in Table 1 are extracted by the OMA method based on GP. Generally, these natural frequencies are too close to tell apart, but the estimated responses are searched for by the OMA method based on GP, as shown in equations (13) to (15). In addition, in Figure 5, the estimated curves fit the reference curves well, and the magnitude of the absolute error is less than 1e⁻³ except for the initial segment as a shock. Thus, the proposed OMA method has the resolution ability to extract NNMs and their modal parameters easily. Figure 6 presents the evolution process of OMA based on SR for responses u₁, u₂, and u₃ (plot a is for response u₁, plot b is for response u₂, and plot c is for response u₃).

$$\begin{array}{l}{\mathbf{u}}_{\mathbf{SR}1}=0.000460527304115933\text{cos}\left(2.36139156091025t\right)\\ +0.0504553957227678\text{cos}\left(3.25509022847277+19.9983152901949\text{t}\right)\end{array}$$

(13)

$${\mathbf{u}}_{\mathbf{SR}2}=-0.0464585117289493\text{cos}\left(5.88836982310503-19.9987243900362t\right)$$

(14)

$$\begin{array}{l}{\mathbf{u}}_{\mathbf{SR}3}=0.00280208085575998\text{cos}(0.296798346254101+0.83469375104736t)\\ +0.00234681861505726\text{cos}\left(2.3228387246822+t\right)\\ +0.0492492754329264\text{cos}\left(3.34162788055375+19.9996416596751t\right)\\ -0.000783966126954913\end{array}$$

(15)

Table 1.

NNMs of the system extracted by the OMA method based on symbolic regression.

Item	NNM 1	NNM 2	NNM 3	Forced element
u1	0	0	0.0004cos(2.3613t)	0.0504cos(3.2550 + 19.9983t)
u2	0	0	0	−0.0464cos(5.8883–19.9987t)
u3	0.0028cos(0.2967+ 0.8346t)	0.0023cos(2.3228+t)	0	0.0028cos(0.2967 + 0.8346t)+0.0023cos(2.3228+t)+0.0492cos(3.3416+ 19.9996t)−0.0007

NNM: nonlinear normal modes.

Figure 5.

Responses $u_1,u_2,\text{and }{u}_3$ for the theoretical reference and estimation by OMA based on SR (symbolic regression).

Figure 6.

Evolution process of OMA based on SR for responses $u_1,u_2,\text{and }{u}_3$ (plot a is for response $u_1$ , plot b is for response $u_2$ , and plot c is for response $u_3$ ).

Comparison

The expression structure of the system response shows that the amplitude of the part of the response signal expression whose frequency is the same as the external force frequency is relatively small compared to the other parts, which have obvious amplitudes. The above performance indicates that the system undergoes a self-excitation phenomenon and that the energy of the system is reassigned to other frequency parts. It is very advantageous to identify the function structure with the OMA method based on GP. Additionally, the response time series are obtained by GP on Eureqa. The function expressions of the response time series are given in equations (13) to (15), and the natural frequency of a single NNM can be easily identified using the coefficients in the function structure.

Three criteria, including the absolute error and phase figure, are selected for comparison of the performance of the proposed nonlinear OMA method, and the results are listed in Table 2. According to the absolute error criterion, GP is suitable for the extraction of an NNM and its modal parameters. Considering the phase figure criterion, the overall NNM identification performance clearly indicates that the GP with symbolic regression method is more suitable than the reference method.

Table 2.

Nonlinear modal parameter estimation with the OMA method using SR.

Modal parameters	Linear theoretical reference	OMA method based on SR
Natural frequency 1 1 w1	0.2251	0.8346
Natural frequency 2 w2	0.3898	1.0000
Natural frequency 3	Null	2.3613

OMA: operational modal analysis.

Because the three natural frequencies are too near to distinguish, the FRF is not selected, and the phase portraits are used instead for characterizing the accuracy and effectiveness of the symbolic regression method. In Figure 7, portraits of the theoretical values and estimations from symbolic regression for u_i (i = 1, 2, 3) are shown. We can also obtain the phase portraits between the estimate ions and the forced response solutions. It can be seen in Figure 7 that the phase portraits for the estimations approach the forced response solutions. Moreover, it is obvious that the portrait from symbolic regression shows a nearly 80% goodness of fit, and the shapes better coincide with the corresponding forced responses compared with the reference.

Figure 7.

Portraits of displacements of response u_i (i = 1, 2, 3) for the theoretical reference and estimations from the OMA method based on SR.

Table 2 shows that the OMA method based on SR can extract the three NNMs well, but the traditional method extracts only two LNMs. Finally, symbolic regression with GP can be applied to structural search and parameter identification very well and is more suitable for nonlinear OMA of nonlinear vibration systems.

The example problem

Rod fastening rotor model

The mechanical model of a simplified posterior rotor ³⁰ is shown in Figure 8. The masses of the two rigid disks are m₁ and m₂, respectively, and their eccentricities are e₁ and e₂. The angle of the two mass eccentricity vectors is $\phi$ . The tensile bar structure has the characteristic of a nonlinear stiffness, and its elastic recovery force can be expressed as $f\left(x\right)=k_{3}x+k_4{x}^3$ . Moreover, the bending stiffnesses of the shafts are k₁ and k₂. The damping coefficients of the shafts and the bending spring are c₁, c₂, and c₃. Only the bending vibration of the rotor model in the x–y plane is considered.

Figure 8.

The simplified mechanical model.

The motion of the disk mass center is expressed as

$$\{\begin{array}{l}{m}_1{\ddot{x}}_1+\left(c_1+c_3\right){\dot{x}}_1-c_3{\dot{x}}_2+\left(k_1+k_3\right)x_1-k_3{x}_2+k_4{\left(x_1-x_2\right)}^3=m_1{e}_1{\omega }^2\text{cos}\left(\omega t\right)\\ m_2{\ddot{x}}_2+\left(c_2+c_3\right){\dot{x}}_2-c_3{\dot{x}}_1+\left(k_2+k_3\right)x_2-k_3{x}_1-k_4{\left(x_1-x_2\right)}^3=m_2{e}_2{\omega }^2\text{cos}\left(\omega t+\phi \right)\\ m_1{\ddot{y}}_1+\left(c_1+c_3\right){\dot{y}}_1-c_3{\dot{y}}_2+\left(k_1+k_3\right)y_1-k_3{y}_2+k_4{\left(y_1-y_2\right)}^3=m_1{e}_1{\omega }^2\text{sin}\left(\omega t\right)\\ m_2{\ddot{y}}_2+\left(c_2+c_3\right){\dot{y}}_2-c_3{\dot{y}}_1+\left(k_2+k_3\right)y_2-k_3{y}_2+k_4{\left(y_1-y_2\right)}^3=m_2{e}_2{\omega }^2\sin \left(\omega t+\phi \right)\end{array}$$

(16)and related parameters of the system are shown in Table 3. Additionally, the predicted frequency response curve is obtained by the proposed OMA method. The response time series are not collected with sensors but replaced with the solution in the time domain obtained using the Runge–Kutta method. The solution is subjected to the proposed OMA method, and the predicted frequency response curve is shown in Figure 9. The absolute errors are all less than 0.01 except for the initial segment as a shock.

Table 3.

Mechanical properties of the rod fastening rotor model.

m1 (kg)	m2 (kg)	Angle $\text{φ}$	$e_1$ (mm)	$e_2$ (mm)
30	33	0	0.15	0.1
Stiffness $k_1$ (N/m)	Stiffness $k_2$ (N/m)	Stiffness $k_{3,4}$ (N/m)	Damping $c_{1,2,3}$ $\left(\text{N}\cdot \text{s}/\text{m}\right)$	Frequency ω (r/min)
0.2e4	4e4	1e4	0.2e4	2π·60

Figure 9.

Response curves of the theoretical reference and absolute error between the reference and the results from OMA with symbolic regression (plot a is the displacement for x₁, x₂, y₁, and y₂; plot b is the velocity for x₁, x₂, y₁, and y₂; and plot c is the acceleration for x₁, x₂, y₁, and y₂).

Results and analyses

The frequency response curve acquired is used as a reference to verify the performance in investigating the response of the system. Hence, the above model, which can accurately capture the structural motion of the system, consists of two NNMs. To further show the performance in the estimation of the NNMs, parameters including the displacement, velocity and acceleration are selected for analysis. Considering that the time interval is sufficiently small, we use the Euler stepping method described in equation (17) to acquire the velocity and acceleration from the displacement. Based on forced responses x_1,2 and y_1,2, symbolic regression is adopted to optimize the estimation of the NNMs for the extraction of modal parameters by OMA. The searched results are as follows: the displacements of responses x₁, x₂, y₁, and y₂ are, respectively, x₁ = 0.144cos(3.36 + 377t), x₂ = −0.103cos(0.101 + 377t), y₁ = 0.144sin(6.06–377t), and y₂ = −0.103cos(1.47–377t); the velocities of responses x₁, x₂, y₁, and y₂ are, respectively, x₁ = 54.4cos(4.93 + 377t), x₂ = 38.9cos(4.81 + 377t), y₁= –54.4cos(0.221 + 377t), and y₂ = −38.9cos(0.102 + 377t); and the accelerations of responses x₁, x₂, y₁, and y₂ are, respectively, x₁ = 54.4cos(0.938–174t)–2.04e4cos(3.22 + 377t), x₂ = 1.46e4cos(6.24 + 377t), y₁ = 2.04e4cos(4.79 + 377t), and y₂ = −1.46e4cos(1.53 + 377t). From the acceleration of response x₁, we can see that the term “–2.04e4cos(3.22 + 377t)” is the forced response, while the term “54.4cos(0.938–174t)” is the nonlinear part extracted by symbolic regression. The forced response appears in all responses above. Because the nonlinear part of the response is weak, we only extract the nonlinear part well from the acceleration of response x₁. Finally, the results are shown in Figures 10 and 11. Moreover, in Figure 10, the good performance of the OMA method with SR, which recovers all the modal amplitudes with good accuracy, is obvious. The OMA method with SR can extract the single NNM well between the left and right edges because the FRF curves bend, and the natural frequency of the nonlinear system changes. In Figure 11, we exhibit the displacement, velocity, and acceleration relationships for the theoretical reference and estimation from the OMA method with SR. The phase diagrams show that the estimated curves from the OMA method with SR are in good agreement with the theoretical reference. The shape is oval. The example shows that OMA using symbolic regression is suitable.

$$\{\begin{array}{l}\mathit{\upsilon }=\left[\mathit{u}\left(\mathit{k}+1\right)-\mathit{u}\left(\mathit{k}\right)\right]/\mathit{h}\\ \mathit{a}=\left[\mathit{\upsilon }\left(\mathit{k}+1\right)-\mathit{\upsilon }\left(\mathit{k}\right)\right]/\mathit{h}\end{array}$$

(17)where h is the sampling time interval.

Figure 10.

Spectral analysis of the responses of x_1,2 and y_1,2 for the reference and their estimation by the OMA method based on SR with harmonic terms (plots a₁, b₁, c₁, and d₁ are, respectively, the displacement of responses x₁, x₂, y₁, and y₂; plots a₂, b₂, c₂, and d₂ are, respectively, the velocity of responses x₁, x₂, y₁, and y₂; plots a₃, b₃, c₃, and d₃ are, respectively, the acceleration of responses x₁, x₂, y₁, and y₂).

Figure 11.

Phase diagrams of the rod fastening rotor model.

Conclusion

One of the greatest challenges in mechanical dynamics is to accurately model the nonlinearity and optimize the modal parameters in realistic vibrating structures using the OMA method. In this paper, a nonlinear OMA method based on symbolic regression is proposed. The natural frequencies of NNMs of the system are identified by using symbolic regression based on the forced response under harmonic excitation. The OMA method using only response time series is an advantageous modal analysis method for parameter identification of nonlinear vibration systems. The existence of nonlinearity is a huge challenge for accurately identifying the modal parameters of vibration systems using OMA. Additionally, there are many deficiencies in the research, and the current research is far from sufficient.