A Volterra-PEM approach for random vibration spectrum analysis of nonlinear systems

In this paper, the Volterra series and the pseudo-excitation method (PEM) are combined to establish a frequency domain method for the power spectral density (PSD) analysis of random vibration of nonlinear systems. The explicit expression of the multi-dimensional power spectral density (MPSD) of the random vibration response is derived analytically. Furthermore, a fast calculation strategy from MPSD to physical PSD is given. The PSD characteristics analysis of the random vibration response of nonlinear systems is effectively achieved. First, within the framework of Volterra series theory, an improved PEM is established for MPSD analysis of nonlinear systems. As a generalized PEM for nonlinear random vibration analysis, the Volterra-PEM is used to analyse the response MPSD, which also has a very concise expression. Second, in the case of computation difficulties with multi-dimensional integration from MPSD to PSD, the computational efficiency is improved by converting the multi-dimensional integral into a matrix operation. Finally, as numerical examples, the Volterra-PEM is used to estimate the response PSD for stationary random vibration of a nonlinear spring-damped oscillator and a non-ideal boundary beam with geometrical nonlinearity. Compared with Monte Carlo simulation, the results show that by constructing generalized pseudo-excitation and matrix operation methods, Volterra-PEM can be used for input PSD with arbitrary energy distribution, not only restricted to broadband white noise excitation, and accurately predict the secondary resonance phenomenon of the random vibration response of nonlinear systems in the frequency domain.


Introduction
In practice, structural vibration is always nonlinear with different levels. Linear dynamic models can be used only for small deformations with small amplitude vibrations or inapparent nonlinear characteristics. More often, nonlinear dynamic models must be used to describe many practical situations, such as thin plate buckling problems [1], hysteresis effects from nonlinear constitutive relations [2], contact nonlinearity in bolted structures [3,4], and geometric nonlinearity from large deflection beams [5]. Nonlinear systems typically exhibit complex mechanical behaviour, such as the dependency of the natural frequency and dissipation on the input excitation amplitude, discontinuities in the frequency response function, and a multi-harmonic response via a single harmonic excitation [6]. It is critical to study and explain the special phenomena that occur in nonlinear systems in order to accurately predict and evaluate the dynamical behaviour of nonlinear systems.
Extensive research work has been carried out for the analysis of the dynamic response of nonlinear systems [7][8][9][10][11][12]. The most widely used method is the direct integration method in the time domain, which has the significant advantage of being able to handle the dynamics of systems with different load types and different strength nonlinearities. However, computational resource consumption is a major concern in practical applications, especially for complex structures. In addition, the cumulative error in time-domain integration cannot be ignored [11]. In the frequency domain, the theory of nonlinear solutions based on series is of interest to many scholars. For example, the harmonic balance method (HBM) was effectively developed using a Fourier series expansion. This is a frequency domain method with high efficiency and accuracy. Nonlinear phenomena such as softening and hardening can be effectively described using the HBM. The HBM is the most feasible method for the steady-state solution of nonlinear systems with bolted complex engineering equipment under periodic loads [12]. In addition, the Volterra series has shown promising applications by extending impulse response functions and frequency response functions in linear systems to nonlinear systems [13][14][15][16][17][18]. In related studies of the Volterra series, Cheng et al. [19] reviewed the research progress of Volterra series theory and its application in mechanical engineering. For the hysteresis nonlinearity of the material, Ran et al. [20] derived the frequency response function (FRF) of the Bouc-Wen model using the Volterra series. They proposed twice transformation method to reduce the difficulty of identifying and gave a Bode diagram of the displacement FRF of the simplified bridge model under three different levels of excitation. Based on a nonlinear time-domain parametric model, Peyton and Yaser [21] proposed a new algorithm for computing high-order multi-input and multi-output Volterra general frequency response function (GFRF).
In recent work, Volterra series theory has also been applied to studies related to the identification of nonlinear systems [22]. It is notable that all the abovementioned extensive and fruitful works constitute dynamical analyses of systems subjected to deterministic excitation. At present, the vibration theory of homogeneous media materials is relatively well developed, but the vibration problem of porous materials, which are widely used in modern industry, has been extremely challenging due to the complex mechanical characteristics. In this field, Abdel-Wahab et al. [23][24][25] first creatively proposed an isogeometric analysis (IGA) method for the dynamic characterization of porous materials (including porous rock materials, and functional gradient nanoplates) and obtained high-precision numerical results using quartic NURBS basis functions.
In practice, many engineering structures with nonlinear characteristics are subjected to random loads, such as earthquakes, waves, and jet noise [26]. Due to the complexity of dynamic problems, the safety requirements of structures under these random loads may bring greater challenges to engineers. To solve the random vibration problem of nonlinear systems, many methods have been established. One of the methods is the Fokker-Planck-Kolmogorov (FPK) approach. However, the FPK equation is difficult to solve analytically, and only a few equations can be solved accurately under some strict assumptions. In this context, the relevant numerical methods were further developed to solve the approximate solution of the stable FPK equation [27,28]. Another class of methods is derived from deterministic nonlinear vibration theory, including the equivalent linearization method [29], Monte Carlo simulation (MCS) [30] and the series expansion method [15][16][17]. One of the more widely used methods is the equivalent linearization method (ELM) based on statistical moments [31]. For example, Oliva et al. [32] used an equivalent linearization technique to analyse the time-domain responses and their PSDs for a system with a nonlinear energy sink (NES) under white noise excitation. Boussaa and Bouc [33] proposed a new linearization technique for analysing PSDs of responses of hysteretic nonlinear systems. Elliott et al. [34] studied the equivalent linearization modelling technique for nonlinear damped systems and analysed the velocity PSD of the system under white noise. That study considered only the case where the random excitation is white noise. The phenomenon of multiple resonance peaks, which is unique to nonlinear systems, has not been studied much. In recent years, Roncen et al. [6,35] creatively developed the HBM for nonlinear random vibration problems. A significant feature of this method is the introduction of the alternate frequency time (AFT) method and nonlinear iterative techniques in the nonlinear solution process and the ensemble averaging of sample PSD as an estimate of the theoretical PSD. In terms of the series approach, Belousov et al. [36] studied the time-domain response of a Duffing oscillator under white noise based on the Volterra series. In the study of random vibrations, the PSD analysis method has always had an important position due to the second-order statistical characteristics of a steady random process can be characterized by the PSD. For the random vibration of linear systems, the response PSD is easily calculated using the FRF. However, the FRF is not directly applicable to the random vibrations analysis of nonlinear systems, which becomes an obstacle to PSD analysis in the frequency domain. For this, a frequency domain method of PSD analysis is developed for a nonlinear random vibration based on the concept of generalized FRF, combined with the Volterra series and the PEM.
For linear systems, the response PSD analysis method has been very well established. As a typical method, the PEM [37,38] transforms stationary random vibration analysis into harmonic vibration analysis and transforms non-stationary random vibration analysis into conventional step-by-step integration analysis. The PEM achieves the same accuracy as complete quadratic combination (CQC) with a significant increase in computational efficiency. In recent years, the PEM has been widely used in random vibration problems in earthquake and wind resistance analysis, and vehicle engineering [39][40][41][42][43]. If the PEM can be used for the response PSD analysis of nonlinear random vibrations, it will be a very interesting topic.
With the above motivation, this paper extends the PEM to a nonlinear random vibration PSD analysis based on the Volterra series and solves two key problems during this process: (1) MPSD analysis: general pseudo-excitation is constructed by the Volterra series, and the MPSD of the random vibration response of nonlinear systems is computed from the general pseudo-response; (2) PSD analysis: the associated multi-dimensional integral is transformed into a matrix operation by the characteristics of the problem, which effectively overcomes the numerical computation difficulties from MPSD to PSD. This also further makes the present method applicable to arbitrarily shaped input PSDs and not limited to white noise or band-limited white noise excitation. The paper is organized as follows: In Sect. 2, based on the Volterra series, a generalized PEM method is established for MPSD analysis. In Sect. 3, the computational procedure of the GFRF is given, which is based on the growing exponential method. In Sect. 4, the PSD of random vibrations for a nonlinear spring-damped oscillator and non-ideal boundary beam considering geometrical nonlinearity is studied. Comparison with the MCS shows that the proposed method has good accuracy and validity for predicting the response PSD of random vibration. The Volterra series develops the concepts of the IRF and FRF of linear systems and applies the concepts to nonlinear systems [13]. With these concepts, a nonlinear system response can be expressed as where h n s 1 ; . . .; s n ð Þis the Volterra kernel function or the general impulse response function (GIRF). Similar to the theoretical framework of linear systems, the norder GIRF h n s 1 ; . . .; s n ð Þ and the n-order GFRF H n x 1 ; . . .; x n ð Þhave the multi-dimensional Fourier transform relationship [44]. The corresponding frequency domain expression of Eq. (1) is 2.2 Generalized PEM for MPSD analysis Statistical information must be used to describe the dynamical behaviour of a system under random loads. For the frequency domain analysis, the important response statistic information is the PSD. In the following section, the generalized PEM is derived for PSD analysis of the random vibration response of nonlinear systems based on the GFRF. The response y t 1 ð Þ, y t 2 ð Þ are expanded in Volterra series according to Eq. (1), respectively, and the autocorrelation function of the response can be obtained [14]: R y m y n t 1 ; t 2 ð Þ ð4Þ R y m y n t 1 ; t 2 ð Þis defined as the partial autocorrelation function; consequently, R y m y n t 1 ; t 2 ð Þ¼R ðmþnÞ y m y n t 1 ; . . .; t mþn ð Þ in which R ðmþnÞ y m y n t 1 ; . . .; t mþn ð Þ is defined as the multi-dimensional partial autocorrelation function [14].
Since the multi-dimensional partial autocorrelation function R ðmþnÞ y m y n and the multi-dimensional partial selfpower spectral density function S ðmþnÞ Y m Y n also have the Wiener-Khinchin relation [14], then the multi-dimensional partial self-power spectral density function can be rewritten as When the input load is a zero-mean Gaussian random process [45], the multi-dimensional partial input autocorrelation function can be written as where is the sum of the product of any two combinations of x t 1 ð Þ Á Á Á x t n ð Þ.
Applying the Fourier transform to Eq. (7), the multi-dimensional partial input PSD can be written as where x i ; x j À Á are any two combinations in When m þ n is even, Eq. (6) can be transformed as When the random process is stationary, let t j ¼ t i þ s, and the 2-dimensional spectrum is as follows: By combining Eqs. (9) and (10), we obtain Using the filtering property of the d function, the number of variables in the m þ n dimensional PSD can be cut in half, i.e., reduced from m þ n to mþn 2 (defining mþn 2 as n), and the integral variable x is changed to X. We obtain S ðnÞ Y m Y n X 1 ; . . .; X n where w 1 ; . . .; w m f g[ c 1 ; . . .; c n f g¼ ÀX n Á Á Á ; ÀX 1 ; È X 1 ; . . .; X n g and X n ¼ X 1 ; . . .; X n È É . To deduce the formula conveniently, the following variables are defined: . . .; c n f g , and X AEn ¼ ÀX n Á Á Á ; ÀX 1 ; X 1 ; . . .; X n È É . Then, Eq. (12) can be expressed as If the following n-order pseudo-excitation is constructed then the corresponding n-order pseudo-response can be expressed as e y n X n ; t Using the property of the m-order GFRF [14] and setting H ¼ ÀW gives where H Ã m is the complex conjugate of H m . Then, the m-order conjugate pseudo-response becomes e y Ã m X n ; t From Eqs. (15) and (17), the n-dimensional PSD can be expressed as where the total number of combinations is ð2nÞ! n ð Þ!2 n . Equation (18) can also be represented by the commonly integral variable x, such as

Calculation from MPSD to PSD
Using multi-dimensional integration, Eq. (19) can be expressed as The response PSD is produced by It should be noted that because the Volterra series is an infinite series, determining its convergence remains a difficult problem [46]. Zhu and Lang [8] derived a new criterion to determine the convergence of the Volterra series representation of the nonlinear system response. This problem is not studied in this paper, and only the truncated Volterra series is used for approximate representation, which already has good accuracy for the first three orders for weak nonlinearity [47].
Based on the first three orders of Volterra series expansions, Eq. (18) is used to compute the self-PSD of the response. When m ¼ 1; n ¼ 1 and n ¼ 1, then Equation (22) is the classic formulation for the response PSD of linear systems.
When x [ 0, by letting x ¼ x and eliminating the d function, Eq. (23) can be transformed into the matrix operational form Here, Dx 1 ¼ 2x u 2Nþ1 , N is the number of sampling points for the matrix operation, and x u is the upper limit of frequency.
Analogously, eliminating the d function, Eq. (26) can be transformed into matrix operational form as and For the cases (e.g., m ¼ 3; n ¼ 1 or m ¼ 1; n ¼ 3), the corresponding matrix operation is not given here, which is similar to Eqs. (24) and (27). For the S Y 3 Y 3 x ð Þ, that consumes a lot of calculation cost, Volterra-PEM only requires a single summation calculation using Eq. (27). If the calculation is performed by the traditional numerical integration method, the double integration will need to be converted into a double summation.
In some practical engineering, velocity and acceleration are the quantities with which people are most concerned. Based on the displacement response PSD, the velocity and acceleration response PSDs can be solved using the following equations: Thus, The independent variable of the PSD function can be represented in terms of the frequency f or circular frequency x. The conversion relationship is

Computational algorithm of the Volterra-PEM
This part gives the details of the computing algorithm of the Volterra-PEM. First, the growth exponential method is used to compute each order of GFRFs. Second, for a given input PSD, the output PSD of a nonlinear system is computed at each frequency point. At each determined frequency, the multi-dimensional integral of the response MPSD can be transformed into the matrix operation. Finally, the response PSD of the nonlinear system at the corresponding frequency point is computed by summing the integration matrix of each order. The algorithm is described in Algorithm 1.

Growth exponential method for calculating GFRFs
In the Volterra series theory, how to solve GFRFs is a difficult and important problem. Because of the dimension disaster (the number of model parameters increases exponentially with the complexity of the input-output relation), the identification of GFRFs is very challenging [48]. When the system equation is known, it can be solved by the harmonic probing method [22], growth exponential method [49], Carleman linearization method [50], etc. If the equation is unknown, it can be identified from the experimental data of the system using the least squares method [51], neural network method [52,53], Bayesian method [54], etc. The growth exponential method is the simplest and most logical of the preceding approaches. It is a very convenient and effective tool for determining GFRFs. The basic principle of the growth exponential method is to assume that the input signal is the sum of a specified number of different growth exponents, and by solving the dynamic response of the system and then categorizing the output signal with the same exponent, the GFRFs of each order can be obtained.
The input and output functions can be replaced by exponential sequences [55]. We can assume that the input is The assumed output is Substituting Eqs. (34) and (35) By making the same exponential coefficients equal and by making k n ¼ ix n , the linear equations about the symmetric GFRFs to be solved are obtained as The subscript ''sym'' in Eq. (37) represents symmetry. For a symmetric GFRF, the properties are as follows: Because the theories in this paper are based on the symmetric Volterra kernel theory, the subscript ''sym'' is typically omitted for brevity [46].

Numerical examples
In this section, numerical examples are given to validate the accuracy and effectiveness of the proposed method. The response PSD is studied for a nonlinear spring-damped oscillator and a non-ideal boundary beam under a stationary random load. In the case of the nonlinear spring-damped oscillator, the results and the computational efficiency are compared and discussed. In the case of the non-ideal boundary beam considering geometrical nonlinearity, the nonlinear phenomena generated by the system and the factors influencing the deviation of the results are analysed and discussed.

Random vibration spectrum analysis
for a nonlinear spring-damped oscillator

Model description
The nonlinear spring-damped oscillator is shown in Fig. 1, and the equation of motion is The characteristic parameters of the nonlinear spring-damped oscillator [51] are shown in Table 1.

GFRFs of the nonlinear spring-damped oscillator
In the analysis of the random vibration of nonlinear systems, the pseudo-response is computed using Eq. (15), which requires the solution of a GFRF. The assumed third-order input is Then, the corresponding output is wðtÞ ¼ G 100 e k 1 t þ G 010 e k 2 t þ G 001 e k 3 t þ G 200 e 2k 1 t þG 020 e 2k 2 t þ G 002 e 2k 3 t þ G 110 e k 1 þk 2 Substituting Eqs. (40) and (41) into Eq. (39), and setting the coefficients of the following items (Eq. (42)) at both ends of the equation equal, The first three orders of GFRFs for the nonlinear spring-damped oscillator can be obtained as Figure 2 shows the first-order GFRF H 1 x ð Þ computed using Eq. (43). We set Fig. 3.

Characterization of the output PSD of the nonlinear spring-damped oscillator
This section discusses the random vibration behaviour of a nonlinear spring-damped oscillator under noise.
To validate the accuracy and effectiveness of the Volterra-PEM, the MCS is used as a comparison method. The MCS uses a triangular series [56] to generate a sample of a random process. When the number of terms of the triangular series is sufficiently large, the sample accurately reflects the probability characteristics of the random process (autocorrelation function or PSD function). In the MCS approach, the expression of the input random process xðtÞ is and x u and x l are the upper and lower limits of the analysed frequency domain, respectively. S XX x k ð Þ is the PSD at the corresponding frequency. N D is a sufficiently large integer, and u k is a random phase obeying a uniform distribution in the range of [0, 2p].
When the time-domain load is applied to the nonlinear spring-damped oscillator, the time-domain response of the system (displacement, velocity, etc.) is solved by the Runge-Kutta numerical integration method. Subsequently, for the time-domain response (displacement, velocity, etc.), the PSD estimation method is used to obtain a sample response PSD. The response PSD of a nonlinear system is obtained by the average of sample response PSDs. The computational algorithm of the MCS for the response PSD computation of nonlinear systems is described in Algorithm 2. The Volterra-PEM (Algorithm 1) and MCS  (Algorithm 2) were used to compute the output PSD of the nonlinear spring-damped oscillator with Eq. (39). Figure 4 shows the three different input PSDs and the corresponding output PSDs. According to the peak intensity, Fig. 4a, c, and e are defined as Input PSD-1 (with a maximum peak of 5 9 10 -3 N 2 /Hz), Input PSD-2 (with a maximum peak of 2 9 10 -2 N 2 /Hz) and Input PSD-3 (with a maximum peak of 2 9 10 -1 N 2 / Hz), respectively. A high level of excitation (bandwidth from 18 to 22 Hz) was designed at approximately 20 Hz (the natural frequency of the nonlinear spring-damped oscillator is approximately 20 Hz) in Input PSD-1, Input PSD-2 and Input PSD-3. The aim of this design was to excite a more significant nonlinear behaviour by applying a stronger excitation at the natural frequency. Figure 4b, d, and f shows the corresponding outputs with three different inputs, defined similarly as Output PSD-1, Output PSD-2 and Output PSD-3. Output PSD-1, Output PSD-2 and Output PSD-3 contain the results from the linear response (ignoring nonlinear terms), Volterra-PEM and MCS (with 200 samples).
For Output PSD-1 (see Fig. 4b), the computed results from the linear response, Volterra-PEM and MCS remain consistent. A very slight peak appears at approximately 60 Hz. Local magnification reveals that the Volterra-PEM predicts a very slight nonlinear response.
In the case of Output PSD-2 (see Fig. 4d), the system output PSD shows the second peak at 60 Hz, which is caused by the nonlinearity of the system. This nonlinear peak is more prominent than that of Output PSD-1. The nonlinear phenomenon of ''a multifrequency band response via a mono-frequency band excitation'' becomes more pronounced as the excitation intensity increases. The Volterra-PEM and the MCS produce very consistent results at these excitation intensities. This confirms the validity of the Volterra-PEM. However, the linear response does not exhibit this phenomenon.
In the case of Output PSD-3 (see Fig. 4f), S Y 3 Y 3 is the source of the phenomena of a multi-frequency band response via a mono-frequency band excitation. In this case, the MCS is closer to the exact solution. The output PSD computed from the MCS shows a decrease in the peak near the main resonance peak at 20 Hz, compared to the linear result.
For computational efficiency, the Volterra-PEM has a significant advantage. Table 2 presents a comparison of the computational times of the Volterra-PEM and MCS. The following personal computer configurations were used during the numerical computations: CPU: Intel Core (TM) i5-9400F, RAM: 32.0 GB. In Algorithm 1, the number of integral points for matrix operation N is 200. In Algorithm 2, the integration method for the nonlinear dynamical equations is the Runge-Kutta method, the sampling frequency Fs is 1000 Hz (periodogram analysis), the integration step Dt is 1 Fs (time-domain analysis), and the total number of computational steps N step is 2 12 .

Computational efficiency analysis of the Volterra-PEM
This part discusses relevant issues for the Volterra-PEM in terms of computational efficiency, which is used to illustrate the significant advantages of the Volterra-PEM in computational efficiency. With Input PSD-2 selected as an example, Fig

Model description
In this part, random vibration PSD analysis of a nonideal boundary beam with geometrical nonlinearity is briefly described. The beam is subjected to a random acceleration load applied by the base. The beam model equations proposed by Nayfeh [57] are used. The material properties of the beam are assumed to be constant along its length. The mechanical model of the beam is shown in Fig. 7. Taking into account the dissipation of the dynamical system, a linear viscous damping term is added. The equation of motion can be expressed as where q is the density, A is the cross-sectional area, l is the linear viscous damping coefficient, E is Young's modulus, I is the cross-sectional moment of inertia, a e is the acceleration at the ends of the clamped beam, t is time, x is the abscissa coordinate, w is the lateral displacement in the midpoint of the beam, and u is the axial displacement. TðtÞ is the tension force in the beam, which is written by Since the structure and its excitation are symmetric, the response w is also symmetric. Only one half of the beam must be studied. The constrained boundary springs k r and k b , the viscous damping coefficient l, and the half-length l are confirmed in experiments by model updating [6]. The values of the model-related characteristic parameters [58] are shown in Table 3. More modelling details can be found in reference [6].
By introducing boundary conditions, transforming to modal coordinates, and defining the set of modes under consideration as M 2 1; 2; . . . f g , Eq. (46) can be expressed as Here, k 2 i is the i-th solution of the equation [59] 1 À Y i can be written as If the frequency range under consideration is within 500 Hz, the first-order mode of the structure will be dominant; i.e., M ¼ 1 f g. Equation (48) can be simplified as For a random vibration problem, the random process can be used to replace the harmonic excitation a e cosðXtÞ.

GFRFs of a non-ideal boundary beam with geometric nonlinearity
Using the growth exponential method described in Sect. 3, the first three-order GFRFs of Eq. (52) are written by Correspondingly, Fig. 8 shows the first-order GFRF H 1 x ð Þ of the system. Letting Fig. 9.

Characterization of the output PSD of a nonideal boundary beam with geometric nonlinearity
For the output PSD analysis of a non-ideal boundary beam with geometrical nonlinearity, two random excitation cases (RMS levels of 2 m/s 2 and 4 m/s 2 ) are considered. Figure 10 presents the input PSDs and the output PSDs of the structure. In the input PSDs (see Fig. 10a, c), from 210 to 230 Hz, two peaks are designed: a primary peak and a secondary peak. This is used to show that the method can be applied to arbitrarily shapes of input PSD.
The output PSD was computed for the non-ideal boundary beam described by Eq. (53) using both Volterra-PEM and MCS (100 samples), respectively. Figure 10a shows the level of excitation of 2 m/s 2 RMS, and the corresponding output PSD is shown in Fig. 10b. The results obtained by the Volterra-PEM and MCS are very consistent. For the input PSD with 2 peaks, the output PSD has 3 peaks. This phenomenon is not present in the corresponding linear response (see the green curve in Fig. 10b). The third peak is the ''triple resonance peak'' generated by the input PSD through the nonlinear system (the frequency corresponding to the nonlinear peak is approximately 380 Hz, which is approximately three times the natural frequency of 126 Hz).
When the excitation RMS level is increased to 4 m/ s 2 (see Fig. 10c), the corresponding output PSD is as shown in Fig. 10d. The generation of ''three nonlinear peaks'' can also be observed. The computed results from the Volterra-PEM and MCS are still in good agreement. However, there is some deviation for the third peak. This is closely related to the level of excitation and is discussed further in a later section for higher excitation strengths.
In Algorithm 1, the number of integration points for matrix operation N is 400. In Algorithm 2, the sampling frequency Fs is 3000 Hz (periodogram analysis), and the total number of integration steps N step is 2 13 . Table 4 gives a comparison of the computational times for the Volterra-PEM and MCS. As can be seen from the table, the computational time for the Volterra-PEM is approximately 9 s for both load excitation cases. However, the MCS  computational time is 47 s. The Volterra-PEM has a very good computational efficiency.

Comparative analysis of the Volterra-PEM and equivalent linearization method
To show the advantages of the Volterra-PEM in PSD prediction of random vibration, the Volterra-PEM and ELM are compared below. For nonlinear systems, the ELM can be used to approximate the system response to a certain extent. Various methods have been developed to determine the equivalent linear coefficients. In this paper, an error minimization scheme is used [34].
For the beam with a non-ideal boundary studied in this section, the system can be modelled as a system with nonlinear stiffness. By introducing the equivalent stiffness coefficients of the nonlinear system, Eq. (52) can be replaced by the following equivalent linear system: where k e is the equivalent stiffness coefficient and the other variables are the same as in Eq. (52). For a random vibration problem, the random process can be used to replace harmonic excitation a e cosðXtÞ. We construct the residuals of Eqs. (52) and (54): For Eq. (54) to be the optimal equivalent linear system of Eq. (52) in the least squares sense, it is necessary that the mean square value E e 2 ½ be minimized, with the necessary condition: Substituting Eq. (55) into Eq. (56), we can obtain: For a zero-mean stochastic process, we obtain: where S ww x ð Þ is the self-PSD of the response w. The equivalent stiffness coefficient k e in Eq. (57) needs to be solved iteratively since it contains an unknown E w 2 ½ . Here, two conditions with RMS levels of 2 m/s 2 and 4 m/s 2 are selected to compare the calculation results of the Volterra-PEM and ELM. During the iterative calculation of the equivalent stiffness coefficient, k e is calculated for each iteration, as shown in Fig. 11.
For the input PSDs with RMS levels of 2 m/s 2 and 4 m/s 2 , the respective output PSDs of the beam with a non-ideal boundary are calculated and shown in Fig. 12. It can be found in Fig. 12 that near the peak of the main resonance, ELM can also exhibit the hardening effect specific to nonlinear beams (i.e., in the output PSD, a rightward shift of the curve is observed). However, the ELM does not predict the nonlinear peaks of the system at triple frequencies. This is due to the fact that ELM is essentially a conversion of a nonlinear system into a linear system with only variance equivalence, which is highly deficient in the prediction of the output PSD. As shown in Fig. 12, the Volterra-PEM can predict the multi-peak phenomenon of the nonlinear random vibration PSD of the system very well. This shows  Fig. 11 Equivalent stiffness for various steps in the iterative method: a input PSD with RMS = 2 m/s 2 ; b input PSD with RMS = 4 m/s 2 that the method in this paper has a considerable advantage in predicting the output PSD. In addition, the response variance of the system obtained by the different methods is given in Fig. 13 according to Eq. (58). At both RMS levels, the Volterra-PEM and ELM are in relatively good agreement. The value of the variance of the linear response in Fig. 13b is slightly smaller, and this part of the deviation is related to the variation of the main resonance peak due to nonlinearity.

Effect of the window function on the PSD estimation
The periodogram method is used to estimate the output PSD of nonlinear systems with MCS. The selection of window functions has a significant impact on the estimation of output PSD. Window functions are an important part of the implementation of spectrum analysis. The window function both corrects for signal non-periodicity and reduces measurement inaccuracies caused by leakage. The fast Fourier transform (FFT) is based on the assumption that the time signal has an infinite period. During the PSD analysis, often only a portion of the time signal is intercepted for analysis, and a window function is added to reduce leakage. Rectangular windows, Hanning windows, and Hamming windows are the most commonly used window functions [60]. The most common one is the rectangular window. Usually, not adding a window function is considered to make the signal pass a rectangular window. The rectangular window is not an ideal one because it is less discriminating of the magnitude of PSD. Both Hanning and Hamming windows are cosine-type windows, which have lower frequency resolution but less leakage of PSD estimates than rectangular windows. Furthermore, the Hamming  Fig. 12 Output PSD obtained by different methods: a input PSD with RMS = 2 m/s 2 ; b input PSD with RMS = 4 m/s 2 window decays at a slightly slower rate than the Hanning window. If the test signal has multiple frequency components and the spectrum is very complex, the Hanning window is usually a good choice.
Taking the level of excitation RMS = 2 m/s 2 as an example, Fig. 14 shows the output PSD estimates obtained using different window functions (using the same set of 100 samples). Good predictions can be achieved at the main resonance peak (Peak 1) with different window functions. However, in the peak regions of Peak 2 and Peak 3, it can be found that an inappropriate choice of window function gives an inappropriate spectral estimation result. Furthermore, it can also be seen from Fig. 14 that the Volterra-PEM can give a good PSD prediction result.

Computational results analysis of the Volterra-PEM with different excitation levels
The random vibration response of the nonlinear system is computed for excitation RMS levels of 0.4 m/s 2 , 2 m/s 2 , 4 m/s 2 and 6 m/s 2 . The output PSD of the nonlinear system under different excitation levels is studied. Figure 15a shows the input PSD at different RMS levels. The corresponding peak values at 114-134 Hz are 0.005 (m/s 2 ) 2 /Hz, 0.2 (m/s 2 ) 2 /Hz, 0.8(m/s 2 ) 2 /Hz and 1.8 (m/s 2 ) 2 /Hz. The peak values at 210-230 Hz are all 6.43 9 10 -4 (m/s 2 ) 2 /Hz. Figure 15b shows the output PSD of the system in the range of 50-500 Hz. The resonant peak (114-134 Hz) and the nonlinear peak (320-440 Hz) increase gradually with increasing excitation intensity. Figure 15c, d shows the local detail amplification of the resonant and nonlinear peaks, respectively. As shown in Fig. 15c, the predictions of the output PSD on the left side of the natural frequency by the Volterra-PEM and MCS are consistent with the increases in the RMS levels. However, the resonant and nonlinear peaks obtained by the MCS are ''widening'' in Fig. 15c, d. This ''PSD widening'' is a particular phenomenon with nonlinear random vibration. The reason for this phenomenon is that the natural frequency of nonlinear systems can change with the excitation, and the strength of the excitation influences the level of ''hardening'' of the nonlinear system.
For RMS = 6 m/s 2 , a very extreme level, the proposed Volterra-PEM has a large deviation at the nonlinear peak of the output PSD at 350 Hz. The Volterra-PEM approach has limitations in this case. Volterra series are more appropriate for weakly nonlinear systems, and the MCS based on timefrequency analysis is considered more effective for strong nonlinearities. For the application of the Volterra series, one of the key points may be the  Fig. 15 Output PSD at different excitation RMS levels: a input PSD; b output PSD; c the local detail at the resonance peak for (b); d the local detail of (b) at the nonlinear peak; the solid line is the result of the MCS, the '' ? '' is the result of the Volterra-PEM, the colours represent different input cases and correspond to (a) GFRF. This may require further development of other methods different from the growth exponential method (e.g., identification methods [54,61]) to overcome to some extent the shortcomings of the Volterra series for strongly nonlinear systems.

Conclusion
In this work, a Volterra-PEM approach was developed to compute the response PSD for nonlinear systems under random excitation. In the framework of Volterra series theory, the response MPSD of a nonlinear system response is transformed into the form of multiplying pseudo-responses. A numerical approach of matrix operation is presented to address the numerical difficulties of multi-dimensional integral from MPSD to PSD. As numerical examples, the Volterra-PEM is used to estimate the response PSD for stationary random vibration of a nonlinear spring-damped oscillator and a non-ideal boundary beam with geometrical nonlinearity. The first three orders GFRFs are obtained by the growth exponential method. The proposed method is a direct theoretical PSD analysis method in the frequency domain, which does not involve the time-domain calculation process and does not require statistical averaging of the sample PSD. Compared with the traditional Volterra series method, by constructing pseudo-excitation and matrix operation method, it can be applied to the arbitrary shape of the input PSD, not only limited to white noise excitation.
PSD analysis of a nonlinear system under a stationary random load is performed using Volterra-PEM and MCS. The results show that the Volterra-PEM, as a frequency domain method, can directly compute the output PSD from the input PSD and does not need any samples. This shows the advantages of Volterra-PEM in efficiency. The proposed Volterra-PEM is a direct theoretical PSD analysis method in the frequency domain, and statistical averaging of the sample PSD is not required in the analysis process. Compared with the more widely used ELM, it can predict the phenomenon of multiple resonance peaks with inherent nonlinearity very well. In addition, compared with the traditional Volterra series method, the present method can be applied to the input PSD of an arbitrary shape by constructing pseudo-excitation and matrix operations and is not restricted to white noise excitation. Under a range of excitation intensities, the proposed Volterra-PEM can accurately predict the response PSD of a nonlinear system, and the Volterra-PEM computational consumption time also has a significant advantage over the MCS method. An interesting phenomenon can also be found in the study: ''a multi-frequency band response via a monofrequency band excitation'' is observed for both nonlinear spring-damped oscillator and non-ideal boundary beam with geometrical nonlinearity under narrow-band random excitation. This nonlinear phenomenon can be understood from a signal analysis point of view as the generation of a frequency component in the response signal that is different from the primary resonance. For such cases, the influence of the secondary resonance frequency component on the structure should be considered, which is of some reference value for structural design and damage diagnosis.
Currently, Volterra series is widely used in different types of engineering nonlinear systems, such as foam vibration isolation systems [10], structures with bolted connections [54], damped nonlinear beams [16], nonlinear oscillators [53], etc. For many specific situations, the nonlinear behaviour of a structure is locally characterized. Modelling engineering nonlinear systems using simplified nonlinear systems, that reflect the nature of dynamics, not only simplifies the problem but also allows for accurate prediction of the dominant nonlinear dynamical behaviour of the system. The above works have conducted relevant studies on this idea. In the numerical example of this paper, the simplified nonlinear system modelling is used, and the analysis of the PSD characteristics and influencing factors of the random vibration of the system is carried out. The engineering background (nonlinear system modelling problem) represented by this simplified model is not covered much. Moreover, at present, the analysis of nonlinear systems using the Volterra series is still at the numerical or laboratory stage and is far from complex industrial cases [54]. If it is used for the design of a practical system, more factors need to be considered [19] and a lot of research needs to be carried out later. Data availability Data generated or analysed during this study are included in this article.