Random vibration of hysteretic systems under Poisson white noise excitations

Hysteresis widely exists in civil structures, and dissipates the mechanical energy of systems. Research on the random vibration of hysteretic systems, however, is still insufficient, particularly when the excitation is non-Gaussian. In this paper, the radial basis function (RBF) neural network (RBF-NN) method is adopted as a numerical method to investigate the random vibration of the Bouc-Wen hysteretic system under the Poisson white noise excitations. The solution to the reduced generalized Fokker-Planck-Kolmogorov (GFPK) equation is expressed in terms of the RBF-NNs with the Gaussian activation functions, whose weights are determined by minimizing the loss function of the reduced GFPK equation residual and constraint associated with the normalization condition. A steel fiber reinforced ceramsite concrete (SFRCC) column loaded by the Poisson white noise is studied as an example to illustrate the solution process. The effects of several important parameters of both the system and the excitation on the stochastic response are evaluated, and the obtained results are compared with those obtained by the Monte Carlo simulations (MCSs). The numerical results show that the RBF-NN method can accurately predict the stationary response with a considerable high computational efficiency.


Introduction
Hysteretic systems have been extensively used, e.g., to describe the constitutive relationship of magnetorheological elastomers [1] , the hysteretic energy-dissipation of structures [2] , seismic isolation element [3] , and shape memory alloy (SMA) [4] . Hysteresis exists in most civil structures, which are often under earthquake or other random excitations. The key structural members or dampers with hysteresis are usually served as the energy dissipator or moment resistant component, whose mechanical properties profoundly affect the seismic performance of the entire structure [5] . Hysteresis is characterized by the nonlinear restoring force of the material, and is a function of both instantaneous and historical deformations [6] . That is to say, the hysteretic force is a function of the displacement and velocity of the structure [7] . It is a challenge to calculate the response of structural systems with hysteresis for the reliability design of structures under random excitations. While the hysteretic system responses under Gaussian excitations have been extensively studied [8][9][10][11][12][13] , much less attention is paid to the responses of hysteretic systems under non-Gaussian random excitations. This paper presents a new radial basis function (RBF) neural network (RBF-NN) method for analyzing the responses of hysteretic systems under Poisson random excitations.
As a common non-Gaussian excitation, the Poisson process is often considered [14][15][16][17][18] . However, studies on the random vibration of hysteretic systems under Poisson excitations are still sparse. Iwankiewicz and Nielsen [19] presented such a study for the first time in 1992. In this work, the response moments of the Bouc-Wen hysteretic system are evaluated by the equivalent linearization technique and equivalent cubic expansion technique. Zeng and Li [20] derived the averaged generalized Fokker-Planck-Kolmogorov (GFPK) equation for Poisson-loaded bilinear hysteretic systems, and obtained the numerical solutions of the response probability density function (PDF). It is clear that more efficient methods for the response analysis of hysteretic systems under non-Gaussian excitations are needed.
Recently, neural networks have shown unprecedented performance in various applications, including solving partial differential equations (PDEs). Earlier studies of this topic can be found in Refs. [21]- [23]. The loss function of the neural network solution is defined as a combination of errors of the equation and the initial and boundary conditions. The minimization of the loss function by numerical optimization methods leads to the solution in the sense of the least square error. The physically informed neural networks (PINNs) are receiving lots of attention lately for solving nonlinear PDEs that obey certain physical laws [24][25][26][27] . Another novel method, graph convolutional networks (GCNs), is used to handle the PDEs with irregular domains [28] . The high-dimensional PDEs can be transformed to backward stochastic differential equations (BSDEs) by the Feynman-Ka formula, and the BSDEs are solved with deep neural networks [29][30][31][32][33][34][35][36][37][38][39] . A deep Ritz method [31] and the adversarial neural network methods [32] are also developed for the solution of general PDEs. High computational efficiency, solution accuracy, broad applicability to different physical systems in higher dimensional systems are the common goals constantly pursued by scholars using neural networks.
The RBF-NN has a simple structure with 3 layers, and is proved to be feasible for solving PDEs efficiently [21] . Since Kansa [33][34] proposed an RBF-NN method for solving PDEs, a number of similar methods have been developed [35][36][37][38] . Recently, the RBF-NNs with Gaussian activation functions have been successfully adopted to obtain the transient and steady state PDFs of the Fokker-Planck-Kolmogorov equation of nonlinear stochastic systems [39] and the reliability function of the first passage problem [40] . The RBF-NNs have a good balance of numerical performance and computational cost, which can be beneficial to high-dimensional problems [41] .
Inspired by the advantages of the RBF-NNs and the encouraging results reported in Ref. [39], we adopt it to solve the three-dimensional (3D) reduced GFPK equation of the Bouc-Wen hysteretic system under Poisson white noise excitations. The effects of intensity, average arrival rate, and other essential hysteretic parameters on the responses are studied. The selection of parameters for the RBF-NNs is also discussed, which is rare in the literature. The potential implementable solutions are also given for high-dimensional problems.
The paper is outlined as follows. The reduced GFPK equation of the Bouc-Wen hysteretic system forced by the Poisson excitations is derived in Section 2. The method for solving the reduced GFPK equation with the RBF-NNs is reviewed in Section 3, followed by a numerical example in Section 4. The parameter selection for the RBF-NNs is discussed in Section 5. Section 6 concludes the paper.
2 Bouc-Wen hysteretic system under a Poisson white noise excitation Consider a nonlinear hysteretic system under a Poisson white noise excitation. The response is governed byẍ where x is the nondimensional displacement. ξ is the viscous damping rate, and α denotes the post-yield to pre-yield stiffness ratio. z represents the hysteretic force. P (t) is the Poisson white noise excitation. The hysteretic force is usually highly nonlinear and history dependent.
Hence, x(t) is not a Markov process. Next, we shall consider a model of the hysteretic force and define an extended state vector, which is Markovian. There are various mathematical models to express the hysteresis characteristics of the structure. The Bouc-Wen model is a very popular one, and is adopted in this paper. The Bouc-Wen model is represented by the following nonlinear differential equation: where A, β, γ, and n are the parameters defining the shape of the hysterical curve. Define a new state vector as (x 1 , x 2 , x 3 ) = (x, y, z), where y =ẋ. Equations (1) and (2) can be written in the state space form as follows: The reduced GFPK equation for the stationary PDF p = p(x, y, z) reads where λ is the average arrival rate of the Poisson process. Y denotes the impulse amplitude of the Poisson white noise. E[·] is the mathematical expectation operator. No exact solution for Eq. (4) has been found up to now. It is common to truncate the equation and solve it numerically. The question is how many terms should be kept in the truncated equation. In the following, we explain a procedure for the truncation of the GFPK equation.
Assume that the impulse amplitude Y follows the Gaussian distribution with zero mean. Then, the coefficients of the items ∂ k p ∂y k in Eq. (4) can be obtained as These terms of Eq. (4) can be viewed as the Taylor series expansion of the function p = p(x, y, z) with respect to the variable y. The convergence of the Taylor series is the necessary condition for the existence of the solution p(x, y, z) in the state space. Hence, we can examine the trend of the coefficients in Eq. (6) as k → ∞ in order to investigate the convergence of the Taylor series.
Note that λ is independent of the index k and can be ignored in the discussion of convergence. Since k is even and The variation of E[Y 2 ] with k trunc for = 0.05 is shown in Fig. 1. It can be observed that for a given k trunc , there is a range of E[Y 2 ] values below the curve satisfying the inequality in Eq. (7). As k trunc increases, implying that more terms are kept in Eq. (4), a larger value of ] is usually small in engineering, a small number of terms to keep in Eq. (4) is sufficient. In the next section, a procedure with RBF-NNs to solve the truncated GFPK equation is presented.

RBF-NN solution of the GFPK equation
Suppose that the RBF-NN solution to Eq. (4) readŝ where N is the total number of the radial basis Gaussian functions covering a domain of interest, and w = [w 1 , w 2 , · · · , w N ] T denotes the weights of the functions. r represents the state vector, µ j denotes the mean of the jth Gaussian function, and σ j is the corresponding standard deviation. The radial form of the Gaussian function reads where m denotes the dimension of the state space, m = 3 for the hysteretic system, x i is the ith state variable, µ ij is the jth mean location of the Gaussian activation function in the ith state direction, and σ ij is the jth standard deviation of the Gaussian activation function in the ith state direction. C j is the normalization constant such that the integration of G(r, µ j , σ j ) over the entire state space is equal to unity. Becausep(r, w) is also normalized such that its integration over the state space is one, we arrive at a constraint on the coefficients w j : For stable stochastic systems, the PDF and its finite order of partial derivatives all vanish as r → ∞. The RBF-NN solution in Eq. (8) satisfies this condition when all the Gaussian functions are distributed in a finite region around the origin.
Very recently, the RBF-NNs in the form of Eq. (8) have been used to approximate the transient solution of the two-dimensional (2D) Fokker-Planck-Kolmogorov equation [39] . The study showed a powerful computing ability of the method. More importantly, it has the potential to solving higher dimensional GFPK. In the following, we apply the solution in Eq. (8) to the 3D stationary GFPK equation (4). We first obtain the residual error of the equation as where m i (i = 1, 2, 3) are the drift terms of the system, and Define a loss function as where r i ∈ D in which D denotes the domain of interest where the center points of Gaussian functions are located, and S ∈ R N ×N is symmetric and is known as the data matrix defined as in which N s is the number of sample points per dimension. The weights w j can be obtained by minimizing the loss function via a proper optimization method subject to the constraint of Eq. (10). To this end, we apply the method of Lagrange multipliers, and consider the following extended loss function: The necessary conditions for minimizing the extended loss function are given by Then, an algebraic equation for the optimal coefficients w is obtained, which can be uniquely determined by the matrix inversion when the data matrix is non-singular. For systems in highdimensional state space, various gradient descent algorithms can be adopted to compute the optimal weights in the framework of neural network training. This is a topic for another study. We point out that when all the weight coefficients w j are positive, the RBF-NN solution in Eq. (8) is guaranteed to be positive. This is a sufficient condition. In most cases, w j is positive, and therefore imposing the constraint w j > 0 has little influence on the PDF solution.
After the solutionp(r, w) =p(x, y, z, w) is obtained, different marginal probability densities of displacement, velocity, and hysteretic force can be computed. Three examples are given as follows: where p 3 (z) describes the steady-state probabilistic distribution of the hysteretic force, which is not commonly available.

Example of a steel fiber reinforced ceramsite concrete (SFRCC) column
The random vibration of an SFRCC column under Poisson white noise excitations is investigated to illustrate the proposed solution method. The parameters of the Bouc-Wen model are obtained from the data of low-cyclic reciprocating loading experiments with the aid of the least square regression [13] as follows: n = 1, A = 1.103 7, β = 0.202, γ = 0.314 7.
The comparison between the theoretical and experimental curves is shown in Fig. 2. The other parameters are ξ = 0.025 and α = 0.2.
For all the cases reported below, we set the number of partitions in each coordinate as N c = 10, resulting in a total number of Gaussian functions centered at the grids N = N 3 c = 1 000. It has been reported that when N s 2N , it leads to sufficiently accurate results [39] . In this paper, we select wherep i and p i represent the probability densities at the ith data point r i obtained by the RBF-NN method and MCSs, respectively. N s denotes the total number of the sample points. We first study the effects of the impulse amplitude E[Y 2 ] and average arrival rate λ on the system response. Figures 3 and 4 show the comparison between the stationary PDFs of displacement, velocity, and hysteresis obtained by the RBF-NN method and MCSs. In the MCSs, the 4th-order stochastic Runge-Kutta algorithm is used, and 10 7 sample points are generated for all cases. It can be observed from these figures that the RBF-NN solutions agree well with the MCSs in almost all the cases. For the hysteretic restoring force, the difference of PDFs between the RBF-NN solutions and MCSs is small and acceptable despite the nonsmoothness of the system. Furthermore, the figures indicate that the PDFs of displacement and velocity display a heavier tail as E[Y 2 ] or λ increases. That is to say, the response of the system has higher probabilities to achieve larger amplitudes with the increase in the excitation intensity. Next, the effect of the post-yield to pre-yield stiffness ratio α is studied. Different values of α are selected, and the numerical results of the stationary PDFs are shown in Fig. 5. Excellent agreement is achieved between the MCSs and RBF-NN solutions. Note that the post-yield stiffness of structural members made of steel or concrete is generally 10% of the initial stiffness [42] , indicating that α can be set to be around 0.1 in most cases. Thus, the proposed RBF-NN method can accurately predict the response of common structures with hysteresis.
The differences of the stationary responses   Fig. 7. Apparently, the solution accuracy of the 4th-order truncation arrives a larger RMS error than that of the 8th-order truncation, impling that a higher-order truncation is indeed necessary in some situations when E[Y 2 ] is not small, as indicated by Fig. 1. However, it is common to truncate the GFPK equation to the 4th-order because a higher order truncation would lead to much higher computational cost [20,43] . It is found that remarkable computational efficiency is achieved with the RBF-NN method even when the higher order truncation is considered. It takes about 3.4 s to solve the 4th-order truncated GFPK and 5 s for the 8th-order one on the computer with Intel(R) Xeon(R) Silver 4210R CPU 2.40 GHz. One of the reasons for such high computational efficiency is that we do not use the Autodiff technique to compute the derivatives [43] and that we obtain the explicit analytical expression of s j in Eq. (14).

Parameter selection for the RBF-NNs
Like in other studies, the parameters for RBF-NNs have a considerable impact on the results. In this section, a guideline for the parameter selection for the RBF-NNs is presented.  Fig. 8 for λE[Y 2 ] = 0.2 and N c = 10. The error above the line R c = R s , i.e., R c > R s , is clearly much larger than that below the line, i.e., R c < R s , suggesting that the domain of the sample points should include the domain of the center points. This observation is in agreement with Ref. [39]. Furthermore, there appears to be a region in the R c R s -plane where the RMS errors are the smallest. This region is the dark area below the line R c = R s , which is approximately described by the inequalities 0.2R s < R c < 0.5R s .

Selection of the standard deviations
The standard deviation σ of the Gaussian functions also plays an important role. To investigate its effect, various values of σ are considered. log(e rms ) of the RBF solution is shown in Fig. 9, which suggests that the best accuracy is obtained when σ = 2R c /N c , i.e., σ is equal to the grid interval for the centers of the Gaussian functions.

Application of multiprocessing
The computational time and solution accuracy of the RBF-NN method as a function of N c and σ are calculated and shown in Fig. 9. A few observations can be made from the figure. First, there exists a minimum RMS error when σ = 2R c /N c = 0.4, as pointed out earlier. Second, as N c increases further, both the RMS error and the computational time start to increase, which is an often-observed phenomenon.
To improve the computational efficiency, we make use of the multiprocessing module provided by Python, which starts multiple subprocesses simultaneously, and executes customized tasks in the subprocesses. We find that the most time-consuming step is the process to create the data matrix S. We divide it into blocks. Each sub-block matrix of S can be regarded as a subprocess. With multiprocessing, several sub-blocks can be computed in parallel and then be assembled to form S. The results in Fig. 9 indicate that multiprocessing can be more beneficial for higher dimensional GFPK equations of multi-degree-of-freedom systems.
In summary, when RBF-NNs are used to solve the GFPK equation, R c can be estimated.
Then, select R s = kR c , where k ∈ (2,5). The domains of the center points and sample points in each dimension are determined as [−R c , R c ] and [−R s , R s ], respectively. The standard deviation of the Gaussian functions σ can be set equal to the grid interval, i.e., σ = 2R c /N c . N c is set to be 10 empirically for the 3D GFPK equation. N s can be selected to be no less than 2N .

Concluding remarks
In this paper, the RBF-NN method is adopted to analyze the stationary response of the structures with the Bouc-Wen hysteresis under the Poisson white noise excitations. An SFRCC column is taken as an example. The RBF-NN solution to the 3D GFPK equation is obtained and validated with MCSs. The influence of various parameters is then investigated to provide a guidance to the application of the proposed method. We have found that the RBF-NN method can solve higher order truncated GFPK equations with high computational efficiency and accuracy. A multiprocessing technology is also adopted that helps the RBF-NN method to improve the computational speed and has a potential to deal with multiple degrees of freedom systems.

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/.