Power spectral density analysis for nonlinear systems based on Volterra series

A consequence of nonlinearities is a multi-harmonic response via a mono-harmonic excitation. A similar phenomenon also exists in random vibration. The power spectral density (PSD) analysis of random vibration for nonlinear systems is studied in this paper. The analytical formulation of output PSD subject to the zero-mean Gaussian random load is deduced by using the Volterra series expansion and the conception of generalized frequency response function (GFRF). For a class of nonlinear systems, the growing exponential method is used to determine the first 3rd-order GFRFs. The proposed approach is used to achieve the nonlinear system’s output PSD under a narrow-band stationary random input. The relationship between the peak of PSD and the parameters of the nonlinear system is discussed. By using the proposed method, the nonlinear characteristics of multi-band output via single-band input can be well predicted. The results reveal that changing nonlinear system parameters gives a one-of-a-kind change of the system’s output PSD. This paper provides a method for the research of random vibration prediction and control in real-world nonlinear systems.


Introduction
Nonlinear systems usually have more complex mechanical behavior than linear systems. One of the most powerful approaches is the Volterra series. Various theories and approaches have been developed to understand nonlinear systems and to tackle the related challenges. The Volterra series is utilized in a variety of fields, including mechanical engineering, electronic engineering, communication engineering, and biomedical engineering, for nonlinear system modeling and identification. Scholars have conducted extensive research on the Volterra series [1][2][3] . Peng and Cheng [4] and Cheng et al. [5] reviewed the achievements and progress of this research field. Under the action of simple harmonic load, a nonlinear system usually shows the influence of super-harmonic, sub-harmonic, and intermodulation interference in the frequency-domain analysis, which makes the analysis and design of nonlinear systems more complex. Jing et al. [3] studied the interaction mechanism of harmonics caused by different nonlinear components in the output spectrum of the Volterra system. Jing and Lang [6] jointly published a monograph, which is a good reference for studying the frequency-domain analysis methods of nonlinear systems.
The behavior mechanism of the nonlinear system under random load has a very wide engineering background. The linearized method based on moment equivalence is widely used [7][8] . Many studies show that the moment equivalence method can only guarantee integral equivalence, but improper results would be given in the prediction of power spectral density (PSD) [9] . Furthermore, the random average approach [10][11] , the Fokker-Planck-Kolmogorov (FPK) equation technique [12][13] , and the Wiener path integral approach [14][15] have been developed. In the field of practical application, the spectral analysis of nonlinear systems in the frequencydomain is very attractive. The theoretical framework of PSD analysis for linear systems has been outstanding. In contrast, due to the inherent difficulties of nonlinear systems, there are few studies on nonlinear systems' frequency-domain responses [16][17][18] . In many cases, the excitation of the actual system is a narrow-band stochastic process. Compared with the system under broadband random excitation, the dynamic behavior of a nonlinear system under narrow-band random excitation would produce more interesting phenomena [19] . But so far, the number of relevant publications on the analysis of output PSD characteristics of the nonlinear system under narrow-band random stimulation has been limited [20][21][22][23] .
Based on the Volterra series, the power spectrum analysis method of random vibration of the nonlinear system is studied in this paper. In the frequency-domain, the transfer relation between the input and output PSD is given by using the generalized frequency response function (GFRF). And the way to convert multi-dimensional PSD into physical PSD is also established. The effect of various parameters on the output PSD of a nonlinear system under a random load is discussed. The paper is organized as follows. In Section 2, the related theories of PSD analysis based on the Volterra series are briefly introduced. Section 3 introduces the growing exponential method and its flow chart. In Section 4, the numerical simulation of the output PSD of a class of nonlinear systems under narrow-band random excitation is introduced, and the causes of various phenomena are analyzed. Finally, the relevant conclusions of frequency-domain PSD analysis of nonlinear systems can be derived. The work of this paper has a valuable reference for the study of random vibration prediction and control of practical nonlinear systems.
2 PSD analysis method of nonlinear random vibration based on Volterra series 2.1 Volterra series for dynamic analysis of nonlinear systems If a system is linear and time-invariant, the system output y(t) subject to input x(t) can be expressed by Duhamel's integral, which is shown as follows: where h(τ ) is the impulse response function (IRF) that describes the response of a system with a unit impulse.
The IRF, h(τ ), and frequency response function (FRF), H(ω), constitute the Fourier transform pairs, which is widely used in the dynamic response analysis of linear systems.
Equation (1) is extended and applied to nonlinear systems. The response of the nonlinear system is expressed as an infinite series [24] , where h n (τ 1 , · · · , τ n ) is the nth-order Volterra kernel, which is the generalized function of IRF. The nth-order Volterra kernel in the frequency-domain is called GFRF [25] .
Similar to the linear system, the nth-order IRF h n (τ 1 , · · · , τ n ) and nth-order GFRF H n (ω 1 , · · · , ω n ) are defined as the multi-dimensional Fourier transform pairs, i.e., ωj τj The corresponding system output in the frequency-domain is where 2.2 PSD analysis of nonlinear random vibration based on Volterra series 2.2.1 Basic theory of PSD analysis The system under random load needs the description of system response statistics. For the frequency-domain analysis, the PSD of system response is most important. Based on the GFRF, the basic theory of PSD analysis for nonlinear systems will be given.
First of all, if m = 1 and n = 1, according to Eq. (12), we get If the input is a stationary random process, let Thus, use Eq. (16) as follows: By combining with Eq. (18) and using Eq. (14), one-dimensional PSD can be expressed as If m = 2 and n = 2, we can obtain The S For the higher-order expectation, such as E(x(t 1 ) · · · x(t 4 )), if x(t) is a Gaussian random process with zero-mean, then it has the following properties [27] If m + n is odd, If m + n is even, In Eq. (23), is the sum of the products of E( If m + n is even, where S (2) XX (ω i , ω j ) means the sum of the product of S (2) XX (ω i , ω j ) with the part (ω i , ω j ) being taken from ω 1 · · · ω m+n (any two combinations). It has (m + n)!/ m+n cases. Then, Eq. (20) can be expressed as Thus, (27) If ω > 0 and δ(ω) = 0, Eq. (27) can be simplified as If m = 3 and n = 3, we obtain There are two cases. In Case 1, The one-dimensional PSD can be derived as The number of times of Case 1 is In Case 2, The one-dimensional PSD can be derived as The number of times of Case 2 is Thus, Form Eq. (24), the following equation can be derived: For the case of m = 1, n = 3 or m = 3, n = 1, the number of times is Thus, Hence, the first 3rd-order PSD of the nonlinear system can be given as Here, o(ω) is the Peano remainder. Note that the Volterra series is infinite, and its convergence discrimination is an urgent problem. This question is still difficult to answer [28] . The latest related research is due to Zhu and Lang [29] , who have deduced a new criterion to judge the convergence of the Volterra series of nonlinear systems. In this paper, the problem of convergence is not studied. The truncation Volterra series is used for approximate expression. Some literature studies show that the first 3rd-order of weak nonlinearity has a good accuracy [30] . For Eq. (15), considering the first 3rd-order of the series and ignoring the cross term (m = n), we get the Volterra series of PSD as follows: In Eq. (42), S Y2Y2 (ω) and S Y3Y3 (ω) can be calculated by residue theorem or numerical integration. In this paper, numerical integration is used.

Calculation of GFRFs
How to obtain the GFRF from the nonlinear system is the difficulty of the research. When the differential equation is known, the harmonic probing method [31] , growing exponential method [32] , and Carleman linearized method [33] can be used to solve the GFRF. If the equation is unknown, the least square approach [17] , neural network [34] , and other approaches based on the system's test data can be used to calculate the GFRF.
The growing exponential approach is the simplest and most logical of the preceding approaches. And the approach can simplify the operation according to the parity of the system equation. It is a simple and useful mathematical approach for determining a nonlinear system's GFRF. The assumption that the input signal is the sum of several independent growth exponents is the key of the algorithm. The desired GFRF of each order can be generated by solving the system's dynamic response and categorizing the output signal with the same exponent. Compared with the direct integration method based on the time-domain analysis, this method can significantly improve the calculation efficiency while ensuring calculation accuracy.
The power series is used to replace the original input and output functions. The assumed input is And the output is Substituting Eq. (43) and Eq. (44) into the nonlinear systems gives the G-function, By making the same exponential coefficient equal, and letting λ N = iω, the linear equations for obtaining the GFRF are obtained as The symmetry is represented by the 'sym' subscript in the formula. The following properties are met by the symmetric GFRF, Because all of the theories discussed in this study are based on the symmetric Volterra kernel theory, the subscript 'sym' will be dropped to save space [28] .

Numerical examples 4.1 Description of problem
Consider the following nonlinear system [3] : The output PSD of the nonlinear system is given by the Volterra series. And the GFRFs need to be solved by using the growing exponential method.
First, suppose that the input is Then, the output corresponds to y(t) = G 100 e λ1t + G 010 e λ2t + G 001 e λ3t + G 200 e 2λ1t + G 020 e 2λ2t + G 002 e 2λ3t Substituting Eq. (50) into Eq. (48), and making the coefficients of the following items at both ends of the equation equal: we obtain the following equations: These problems are solved, and the first 3rd-order GFRFs are obtained according to Eq. (46), The PSD analysis of nonlinear random vibration based on the Volterra series subject to narrow-band stationary random excitation is considered. The calculation flow is shown in Fig. 1. ...

Fig. 1 Calculation flow of PSD analysis based on Volterra series
The input x(t) in this study uses a narrow-band stationary Gaussian random process since the bandwidth of the random input is limited in some realistic engineering applications. For the input x(t), the frequency bandwidth is 5.5 rad · s −1 to 7.5 rad · s −1 , and the value of PSD is 10.
The direct integration method is utilized for computation comparison to verify the proposed method. The trigonometric series formula [35] is used to imitate the random input.
where ω k and ω l are the frequency's lower and upper limits, S XX (ω k ) is the PSD of the corresponding frequency point, N is a sufficiently high integer (in this work, N = 1 000), and ϕ k is a random number obeying the uniform distribution rate in the range of [0, 2π]. Figure 2 depicts the time-domain input signal. The sampling rate is 3.3 times the maximum angular frequency of ω max (ω max = 30 rad · s −1 ). The periodogram code in MATLAB is used to estimate the PSD of the input signal x(t). The PSD of the input signal is obtained as shown in Fig. 3. Several simple cases of a nonlinear system (i.e., Eq. (48)), will be explored to clearly explain the properties of this system. The proposed method is compared with the direct integration method (using the MATLAB code ode45) in this paper. For the direct integration method, the output PSD is calculated on the mathematical software MATLAB combined with the periodogram method. The following are the PC configurations. GPU: NVIDIA GeForce GT 1030; CPU: Intel Core (TM) i-9400f, RAM: 32.0 GB.

Discussion on different nonlinear parameters 4.2.1 Case 1 (d = 0)
First, the case c = 1, a = b = d = 0 is discussed. In this case, the nonlinear system degenerates into a linear system, and the output PSD is shown in Fig. 4. The results of the three methods are consistent. It can be seen that the frequency span of output PSD is from 5.5 rad · s −1 to 7.5 rad · s −1 , which is consistent with the input PSD. For a linear system, it is easy to estimate the output PSD by using the linear FRF of the system, which has been widely studied and applied.
Consider the instance a = 1, b = c = d = 0 as illustrated in Fig. 5, where 2 peaks emerge from the output PSD. The output PSDs from the three methods are consistent. Due to the intermodulation of the input signal, the span of the first peak is from 0 to 2 rad·s −1 (7.5−5.5 = 2 rad · s −1 ). The second peak appears in the range of (5.5 to 7.5) × 2 = (11 to 15) rad · s −1 for the square nonlinear parameters. Furthermore, the second peak changes from flat to spike. The reason is that the frequency band contains not only double frequency but also intermodulation from the first peak (2 rad · s −1 ). Consider the case b = 1, a = c = d = 0 as illustrated in Fig. 6, where the output PSDs from the three methods are consistent. The output PSD has 2 peaks. The first peak span is from 3.5 rad · s −1 to 9.5 rad · s −1 . The reason is that the signal interferes with intermodulation. The left boundary of the first peak is 5.5 × 2 − 7.5 = 3.5 rad · s −1 , and the right boundary is 7.5 × 2 − 5.5 = 9.5 rad · s −1 . The frequency width of the output PSD is three times that of the input (6 rad · s −1 ). The second peak span is from 16.5 rad · s −1 to 22.5 rad · s −1 ((5.5 to 7.5)× 3 rad · s −1 ), which is due to the cubic nonlinear parameters. In addition, this frequency band also contains the effect of signal intermodulation interference.
In the case a = b = 1, c = d = 0, the output PSD is shown in Fig. 7. There are four peaks, which are equivalent to the sum of Fig. 5 and Fig. 6. The output PSDs from the three methods are consistent. In the case a = b = c = 1, d = 0, the PSD of the output signal is shown in Fig. 8. There are also four peaks, which are equivalent to the accumulation of the values in Fig. 4, Fig. 5, and Fig. 6. There are some differences in the output PSDs from the three methods. From Fig. 8, ignoring the cross-terms would make the output PSD smaller. In this case, ignoring the cross term would bring some errors. The PSD of each order based on the expansion of the Volterra series is shown in Fig. 9.  In the case a = b = c = 1, d = 0, the values of output PSD generated by the three methods are shown in Table 1. The results are slightly different. Moreover, this work pays more attention to the nonlinear phenomenon of single-band input and multi-band output.

Case 2 (d = 0)
Consider the case a = b = 0, c = 1, d = 0.005, which is similar to the 'nonlinear stiffness' in a dynamic system. The output PSDs are shown in Fig. 10. And there are 2 peaks at 5.5 rad · s −1 to 7.5 rad · s −1 and 17 rad · s −1 to 21 rad · s −1 , respectively. The second peak span is from 17 rad · s −1 to 21 rad · s −1 , which is due to the third nonlinear term in the system parameters. When d = 0.01, the output PSDs are shown in Fig. 11. The prediction for the first peak is well., while the error at the second peak is larger. The phenomenon indicates that with the increase of cubic coefficient, the nonlinearity would be strengthened, and the fitting accuracy of the Volterra series method would be weakened. Table 2 shows the PSD obtained by the Volterra series method and the direct integration method. The difference of the results increases with the increase in d. With the enhancement of the input PSD, the errors are increasing, which probably requires more terms in the expansion of the Volterra series [36] .

peaks input PSD
In this section, an example of 2 peaks input PSD is discussed. The selected nonlinear case is c = 1 and d = 0.01. The input and output PSDs with 2 peaks are shown in Fig. 12. The output PSD is calculated by four methods. It can be seen from Fig. 12 that the presented method is consistent with the direct integration method.

Spike input PSD
In this section, an example of spike input is discussed. The selected nonlinear case is b = 1. The input and output PSDs are shown in Fig. 13. It can be seen from Fig. 13 that the presented method is consistent with the direct integration method. Volterra series can make a good prediction for this case. But with the complexity of the input PSD shape, the relationship between the peak of PSD and the parameters of the nonlinear system would become complex and not easily explained.

The influence of bandwidth of input PSD
One consequence of nonlinearities is a multi-harmonic response via a mono-harmonic excitation [37] . A similar phenomenon also exists in random vibration. In this paper, more attention has been paid to this phenomenon. Compared with the system under broadband random excitation, the behavior of a nonlinear system under narrow-band random excitation would be more obvious. In Fig. 14, subgraphs 14 (a) Input with bandwidth 3 rad · s −1 to 5 rad · s −1 (b) Output with bandwidth 3 rad · s −1 to 5 rad · s −1 (c) Input with bandwidth 3 rad · s −1 to 10 rad · s −1 (d) Output with bandwidth 3 rad · s −1 to 10 rad · s −1 (e) Input with bandwidth 3 rad · s −1 to 15 rad · s −1 (f) Output with bandwidth 3 rad · s −1 to 15 rad · s −1 S XX (ω) S YY (ω) linear phenomenon of single-band input and multi-band output becomes less obvious with the increase of bandwidth.

Conclusions
Using the Volterra series approach, the output PSD characteristics of the nonlinear system under the narrow-band random process are obtained. The output PSD and the system's first 3rd-order GFRFs are given as expressions. The results show that changing the nonlinear system parameters causes a complicated phenomenon in the system's output PSD. This influence must be considered in the PSD analysis of the nonlinear system. By comparing different methods, the nonlinear characteristic of multi-band output via single-band input is well predicted.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.