Approximate stochastic response of hysteretic system with fractional element and subjected to combined stochastic and periodic excitation

A method based on statistical linearization is proposed, for determining the response of single-degree-of-freedom hysteretic systems endowed with fractional derivative element and subjected to combined periodic and white/colored excitation. The method is developed by decomposing the system response into a combination of a periodic and of a zero-mean stochastic components. In this regard, first, the equation of motion is cast into two coupled fractional-order nonlinear differential equations with unknown deterministic and stochastic response components. Next, the harmonic balance method for the fractional-order deterministic equation and the statistical linearization for the stochastic equation are used, to obtain the Fourier coefficients of the deterministic response component and the variance of the stochastic response component, respectively. This yields two sets of coupled nonlinear algebraic equations which can be solved by appropriate standard numerical method. Pertinent numerical examples, including both softening and hardening Bouc–Wen hysteretic system endowed with different fractional-orders, are used to demonstrate the applicability and accuracy of the proposed method.


Introduction
Many mechanical and structural systems subjected to severe dynamic loads exhibit hysteretic behavior [1]. The hysteretic restoring force not only depends on the instantaneous system response, but also depends on the response history. Therefore, compared to the zero-memory nonlinear elements with elastic behavior, the hysteretic force-displacement relationships exhibit multi-valued loops [2]. To describe the historydependent nonlinear behavior of the hysteretic systems, various parametric models have been proposed. A dry friction-based hysteretic bilinear model is found in Ref. [3,4]. Proposed by Bouc [5] and improved by Wen [6], the Bouc-Wen hysteretic model and its extensions [7] has gained wide popularity among researchers, due to its versatility and simplicity; also see a comprehensive literature review in [8] and a book [9] devoted to various application of this celebrated model. A proper hysteretic model can represent various of engineering structural/material hysteresis behaviors. In this regard, the Bouc-Wen model belongs to a class of smooth analytical models described by differential equations. This yields a mathematical-based framework to develop analytical knowledge on relevant properties, which can be used for multi-purposes in engineering applications. In this regard, Chassiakos [10] et al. developed an online identification method for hysteretic systems using the Bouc-Wen model. Later on, the Bouc-Wen model was utilized for representing force-displacement relationship of MR dampers [11], lead-core rubber isolators [12] and steel dampers [13] and for describing constitutive law of smart materials [14].
Recently, a useful mathematical tool named fractionalorder calculus has been widely used in various engineering applications, especially for mechanical modeling of visco-elastic materials [15]. In this context, one may argue that if visco-elastic relaxation test of the material is well-fitted by a power law decay, then the fractional constitutive law in Caputo's derivative sense naturally appears [16]. Compared to the standard linear solid model, the fractional derivative model can describe the broadband of visco-elastic materials with small number of parameters [17]. Further, fractional derivatives has proved particularly useful in the fields of structural engineering, where visco-elastic materials are used in vibration mitigation devices for structures. Certain examples are related to visco-elastic dampers for connecting shear walls [18] and adjacent stories [19] of buildings. In this regard, the response determination of a structural dynamic systems endowed with fractional derivatives and subjected to random excitation arise naturally; see [20,21] for example.
For the engineering applications where the viscoelastic controlled/isolated structures or/and the involved visco-elastic dampers/isolators behave hysterically, one need to study a hysteretic dynamic model with fractional derivative. This kind of model has a ratedependent force-displacement relationship, which can be used for modeling visco-elastic controlled/isolated structures subjected to server external excitation. Very recently, Kang [22] et al. presented a new fractionalorder normalized Bouc-Wen model to describe the asymmetric and rate-dependent hysteresis nonlinearity of piezoelectric actuators.
From the preceding literature review, one may conclude that very few researchers focus on the response of a hysteretic dynamic system endowed with integer/fractional-damping element subjected to purely periodic or stochastic excitation. However, for some engineering applications, structures/systems are subjected to combined periodic and stochastic excitation. For example, aircrafts with rotation mechanism [23] are often subjected to the mixture of colored stochastic and harmonic excitation. Other examples include response of a gear system [24,25] and airfoil model [26] under stochastic excitation with uncertain disturbance, and so forth. More recently, nonlinear systems subjected to combined harmonic and stochastic excitation attract increasing interest of researchers to account a more realistic operation condition in energy harvest application [27,28]. The preceding literature reviews necessitate the following research regarding the response determination of a hysteretic system endowed with fractional element and subjected to combined periodic and stochastic excitation.
In this paper, a method based on statistical linearization(SL) [29][30][31] is proposed, for determining the response of single-degree-of-freedom (SDOF) hysteretic systems endowed with fractional derivative element and subjected to combined periodic and white/colored excitation. It can be regarded as an extension of the recent developed method for the response determination of a integer-order hysteretic dynamic system subjected to combined harmonic and white noise excitation [32]. This is achieved based on decomposing the system response into a summation of deterministic and stochastic component [33]. Accordingly, the original fractional-order equation of motion can be cast into two coupled fractional-order differential equations, in terms of the deterministic component and the stochastic component, respectively. The harmonic balance method and the SL method are then utilized for these sub-equations of motion, leading to two sets of nonlinear algebraic equations, which can be solved simultaneously using standard numerical schemes, such as Newton's iteration method. The applicability and accuracy of the proposed method is demonstrated by comparing with pertinent Monte Carlo simulations.

Mathematical formulation
Consider a SDOF hysteretic system endowed with fractional element and subjected to a combined periodic and stochastic excitation where x (t) andẍ (t) is the displacement and acceleration of the system, respectively; m, c and k is the coefficient for mass, damping and stiffness, respectively; D C (·) represents the fractional order in Caputo's sense Note the damping coefficient c = 2ζ mω 2−q n to ensure a consistent dimension of the fractional-order damping coefficient with the one of the integer-order derivative, where ω n = √ k/m is the natural frequency and ζ is the damping ratio. The superscript q denotes the order of fractional derivative; α is usually referred to as the 'rigidity ratio'; w (t) is a white/colored stochastic excitation with power spectrum density S w (ω, t); f (t) is a deterministic periodic excitation, written in a form of Fourier series with A n and B n being the Fourier coefficients of the nth harmonic term; z (t) is the hysteretic displacement governed in a differential forṁ As an archetypal hysteretic model, the Bouc-Wen model has been widely used in structural engineering applications. In this context, consideṙ where A, n, γ and β are the Bouc-Wen hysteretic parameters. In particular, when n = 1, Eq. (5) becomeṡ Assume next the steady-state response of the system, governed by Eq. (1), can be cast into a combination of a deterministic and of a stochastic component. That is wherex (t) andẑ (t) are the zero-mean stochastic processes; μ x (t) and μ z (t) are the deterministic mean processes which can be further written as a Fourier series where C n , D n , V n , D n are the Fourier coefficients of response. Substituting Eqs. (7a)-(7b) into the equation of motion shown in Eq. (1) yields Taking mathematical expection on both sides of Eq. (9) leads to Substracting Eq. (10) into Eq. (9) yields Similary, substituting Eq. (7) into the auxiliary firstorder differential equation governing the hysteretic term shown in Eq. (6) yieldṡ Taking mathematical expectation on both sides of Eq. (12) leads to Substituting Eq. (13) from Eq. (12) yieldṡ Note that these two sets of differential equations (Eqs. 10-11 and 13-14) are coupled, because the expectations of the nonlinear functions include both μ x , μ z and σ˙x , σˆz. Therefore, a simultaneous solution of them should be adopted. In the following part, the harmonic balance method is utilized for solving the Fourier coefficients of the deterministic component, while the statistical linearization method is used for the variance/co-variance of the stochastic component.

Harmonic balance method for the deterministic component
An assumption often used in the statistical linearization is the Gaussian distribution of response. In this regard, the expectation in Eqs. (13)- (14) can be evaluated in an approximate closed-form [32]. That is where σ˙x and σˆz are standard deviation ofẋ andẑ, respectively; ρ is the correlation coefficient ofẋ and z.
Taking fractional-, first-and second-order derivative on both sides of Eq. (8a) yields and −C n ω 2 n cos ω n t + D n ω 2 n sin ω n t , respectively. Substituting Eqs (3), (8a), (8b), (16a) and (16c) into Eq. (10) and taking harmonic balance yield − mC n ω 2 n + cω q n C n cos Substituting Eqs. (16b) and (8b) into Eq. (18) leads to an equation with both sides being summation of weighted harmonics. For multi-harmonic representation of μẋ and μ z , one may refer to [34] for an analytical solution of the weighting factors. For a special case when the periodic excitation is mono-harmonic, one may adopt an approximation that the response μẋ and μ z is also mono-harmonic. In this regard, consider and for an illustrative example, where U 0 and V 0 , and C 0 and D 0 are Fourier coefficients for μ z and μ x , respectively. Equation (18) reduces to The weighting factors and are derived from the Fourier series of μ z μ 2 x and μẋ μ 2 z . That is where coefficients of frequency 3ω 0 are omitted. Further, taking harmonic balance on both sides of Eq. (20) gives Combining Eqs. (17a)-(17b) (n = 0) and Eqs. (23a)-(23b) yields a solution of U 0 , V 0 , C 0 and D 0 in terms of σ˙x , σˆz and ρ. However, σ˙x , σˆz and ρ are unknowns. Therefore, more relationships between the unknown Fourier coefficients and response second moments are needed to supplement the set of nonlinear algebraic equations shown in Eq. (17a)-(17b) and (23a)-(23b). This can be achieved by the statistical linearization method, shown in the following section.

Statistical linearization for the stochastic component
Let q = {x,ẑ} T , then the equation of motion of the zero mean responsex andẑ shown in Eqs. (11) and (14) can be cast into where are the mass, integer-order damping, fractional damping and stiffness matrix, respectively; is the nonlinear term.
Next, according to the statistical linearization method, Eq. (24) can be equivalently written as [35] where Further manipulation yields The expectation of nonlinear functions in Eq. (27a)-(27b) can be approximated as [32] E Combining Eqs. (28a)-(28d) with Eqs. (27a)-(27b), one may find that the equivalent parameters of the system are time varying, and for the mean processes, μẋ and μ z are included. Note that the standard deviations σ z and σ˙x are cyclostationary when t → ∞. Therefore, the fast varying content can be eliminated by averaging over the fundamental period T 0 = 2π/ω 0 . That is where the deterministic averaging on the right-hand side of Eq. (29a)-(29b) yields the following approxi-mations: [32] E zsgn (ẋ) Note that the approximate equivalent parametersc e and k e depend on seven unknowns C 0 , D 0 , U 0 , V 0 , σ˙x , σˆz and ρ. Further, the supplemented relationship between the unknown response statistics σ˙x , σˆz and ρ, and the Fourier coefficients of the deterministic response component can be obtained by the standard linear random vibration theory for dynamic systems with fractional elements [36]. In this regard, where is the two-sided PSD matrix of Q and q, respectively; S w (ω) is the power spectral density (PSD) of the stochastic excitation ; Sxx and Sˆzˆz is the two-sided auto-PSD ofx andẑ, respectively; Sxˆz = S * zx is the twosided cross-PSD of (x,ẑ), where the symbol * denotes complex conjugate. The frequency response matrix in Eq. (31) can be written as The covariance σ 2 x and σ 2 z and the correlation coefficient ρ can be derived by

Implementation procedures
To summarize, the harmonic balance method for the deterministic response component leads to four indeterminate nonlinear algebraic Eq.
where the superscript (i) denotes the ith iteration step and ε is a predefined error.
For applying Newton's iterative method in Step 2, one needs to solve the following matrix equation for the unknown Fourier coefficients of the harmonic compo-nent. That is The column vector K with four entries can be written as the Jacobin matrix is Pertinent Monte Carlo simulation is used to validate the proposed method. The specific procedures of Monte Carlo simulation is the following: 1. Use the spectral representation method [37,38] to generate n white/colored noise samples according to the target spectrum. Substitute every white/ colored noise sample into the right-hand side of Eq. (1) to obtain n combined excitation samples. 2. Use Newmark-β method of the incremental form for nonlinear systems endowed with fractional derivative element to solve a sample response of the considered system under a combined excitation sample obtained in Step 1. Repeat the preceding numerical algorithm for all the combined excitation samples to get n response samples. 3. Find the mean value and standard deviation of these responses, which is the deterministic response component and the standard deviation of the stochastic response component, respectively.
To the authors' best knowledge, very few attention has been devoted to the numerical algorithms for the deterministic response of a hysteretic dynamic system endowed with fractional derivative elements. For numerical algorithms of linear fractional-order dynamic systems, one may see Ref. [39]. Therefore, the authors adopt a discretization scheme for the Riemann-Liouville definition of fractional derivative [39] and the Newmark-β discretization scheme for integer-order derivatives [40]. The adopted algorithm for the sample response calculation is developed in an incremental form, and the Newton-Raphson iteration is used in every time-step for hysteretic nonlinearity is considered in the following numerical example. The work regarding the numerical scheme has been recorded as an independent manuscript and will be submitted soon to an international journal.

Numerical example
When α = 1, F 0 = 0, q = 1, the fractional-order hysteretic system reduces to a special case of the linear system with integer-order derivative and subjected to stochastic excitation only. The displacement standard deviation of the reduced linear system is σx = σ x = 1.
In this illustrative case, to demonstrate the applicability of the proposed method, the amplitude and the frequency of the harmonic excitation is chosen as F 0 = 1 and ω 0 = 1.2, respectively. Both softening and hardening Bouc-Wen systems with fractional order q = 0.5 are considered herein. For the softening Bouc-Wen system, the hysteretic parameters are A = 1, γ = 0.5, β = 0.5, n = 1, α = 0.1, whereas for the hardening Bouc-Wen system, all the other hysteretic parameters are the same as those in the softening case except β = −0.35, γ = 0.65. The response comparison shown in Fig. 1a-b for the softening Bouc-Wen model verifies the accuracy of the proposed method. Specifically, the deterministic component time history and the stochastic component standard deviation of the displacement x (t) obtained by the proposed method is in good agreement with the pertinent Monte Carlo estimates over 10,000 samples. Besides, the proposed method yields a stationary standard deviation of the stochastic component, vis-a-vis the MC estimates gives a cyclostationary standard deviation. This is caused by the coupling effect between the deterministic and random response component, shown in Eqs. (28a)-(28d). These equations indicate that the equivalent parameters, and thus, the standard deviation of the stochastic component, depend on deterministic component. However, in the proposed analytical method, the oscillation of the equivalent parameters, and thus, of the cyclostationary variance, is eliminated by averaging over one period for further simplification (see Eqs. (29a)-(29b)). Therefore, the response variance during the cyclostationary phase obtained by the Monte Carlo data is averaged for the ensuing error analysis. Specifically, for the deterministic response μ x (t) shown in Fig. 1a, the response amplitudes obtained by the proposed method is around 8.89% less than the Monte Carlo estimates, whereas for the stationary standard deviation shown in Fig. 1b, the analytical result is about 11.68% less the MC data. Further MC investigations show that the amplitude of harmonic-like cyclostationary part of the standard deviation depends on the amplitude and the frequency of the harmonic excitation. Besides, it seems that the oscillation frequency of the cyclostationary standard deviation is twice as higher as the harmonic excitation frequency. This conclusion is reasonable because the equivalent stiffness, shown in Eq. (27b), depends on quadratic terms of the deterministic response component.
Numerical investigations, as shown in Fig. 2a-2b, demonstrate the applicability of the proposed method for the considered hardening Bouc-Wen model endowed with fractional element. Error analysis sug- Specifically, for the displacement response x (t), the amplitude of the deterministic component and the standard deviation of the stochastic component, obtained by the proposed method, is 16.55% and 10.5% lower than the corresponding MC estimates, respectively. All the above errors are within the reasonable range of the classic statistical linearization method reported by previous researches [35]. Sample hysteretic loops of the considered softening and hardening Bouc-Wen system endowed with fractional element under combined harmonic and white noise excitation is shown in Fig. 3a-b, respectively. It seems that for both circumstances, the hysteretic loops exhibit history-dependent strong nonlinearity. Besides, compared to the hysteretic loops under harmonic excitation, the white noise disturbance renders the loop center move around the equilibrium position. The ensuing analyses further investigate the applicability of the proposed method under other harmonic excitation with different frequencies. Consider the softening Bouc-Wen model with the same parameters as in Fig. 3a. Figure 4 shows the response standard deviation σx , σ˙x , σˆz versus excitation frequency curve when The amplitude and frequency of harmonic excitation are selected to be F 0 = 1 and ω 0 = 1.2, while both softening and hardening system with the same hysteretic parameters as in Figs. 3a-b are considered herein for illustration. Figure 6a and b shows the response standard deviation σx , σ˙x , σˆz versus fractional order curve for the softening and for the hardening system, respectively. It seems that for the considered conditions with fractional order q ∈ [0.1, 0.9], the response statistical characteristics obtained by the proposed method  Figs. 7a-b, comparisons regarding the fractional-order versus deterministic response (μ x , μẋ , μ z ) amplitude curve obtained by the two method demonstrate that for both softening and hardening system, the proposed method agree well with Monte Carlo simulation.

Colored noise as the stochastic excitation
The proposed method can be extended for response determination of a hysteretic systems endowed with where ξ g = 0.2 and ω g = 2 is the damping ratio and natural frequency of the overlying soil of the local seismic site ; S 0 = 2ζ /π is the power spectral density of the white noise generated by the bedrock to be filtered. In this illustrative case, to demonstrate the applicability of the proposed method, the amplitude and the frequency of the harmonic excitation is chosen as F 0 = 1 and ω 0 = 1, respectively. Consider both a softening and a hardening Bouc-Wen system. For the softening Bouc-Wen system, the hysteretic parameters are A = 1, γ = 0.5, β = 0.5, n = 1, α = 0.1, whereas for hardening Bouc-Wen system, all the other hysteretic parameters are the same as those in the softening case except β = −0.1, γ = 0.8. The response comparison shown in Figs. 8a-b for the softening Bouc-Wen system verifies the accuracy of the proposed method. Specifically, for the deterministic harmonic response x (t) shown in Fig. 8a, the amplitude of the deterministic response component, obtained by the proposed method is around 3.83% less than the Monte Carlo estimates. The stationary standard deviation of the stochastic response component calculated by the proposed method, shown in Fig. 8b, is about 11.63% less than the time average of the cyclostationary standard deviation obtained by the Monte Carlo simulation. Further MC investigations show that the amplitude of the harmonic-like cyclostationary part of the standard deviation depends on the amplitude and the frequency of the harmonic excitation.
Numerical investigations, as shown in Fig. 9a-b, demonstrate the applicability of the proposed method for the considered hardening Bouc-Wen system endowed with fractional elements. Error analysis suggests satisfactory accuracy of the proposed method. Specifically, for the displacement response x (t), the amplitude of the deterministic component and the standard deviation of the stochastic component, obtained by the proposed method, is 11.32% and 5.91% less than the corresponding MC estimates, respectively. All the above errors are within the reasonable error range of the classic statistical linearization method.
The same conclusion as in the white noise case can be obtained that the sample hysteretic loops of the considered hysteretic system under combined harmonic and colored noise excitation, exhibit history depended strong nonlinearity. Besides, compared to the hysteretic loops under harmonic excitation, the colored noise disturbance renders the loop center move around the equilibrium position.
The ensuing analysis further investigates the applicability of the proposed method under other harmonic excitations with different frequencies. Consider the softening Bouc-Wen model with the same parameters as in Fig. 8a-b. Figure 10 shows the response standard deviation σx , σ˙x , σˆz versus excitation frequency curve, when the harmonic excitation amplitude is F 0 = 1. It seems that the proposed linearization method exhibits reasonable accuracy (< 18%) over all the considered excitation frequencies for the derived displacement standard deviation compared with the Monte Carlo estimates. The amplitudes of the deterministic response component versus the harmonic excitation with different frequencies, obtained by the proposed method and estimated by the pertinent Monte Carlo simulation, are  Fig. 11 when F 0 = 1. It seems that under the considered conditions, the amplitude of the deterministic response component versus harmonic excitation frequency curve can be well-predicted by the proposed analytical method. The last analysis investigates the applicability of the proposed method for the considered Bouc-Wen systems endowed with different fractional orders. The amplitude and frequency of harmonic excitation are selected to be F 0 = 1 and ω 0 = 1, while both softening and hardening system with the same parameters as in Figs. 8a-b and 9a-b , respectively, are considered herein for illustration. Figure 12a-b shows the response standard deviation σx , σ˙x , σˆz versus fractional order curve for softening and hardening system, respectively. Figure 13a-b shows the deterministic response (μ x , μẋ , μ z ) versus fractional order curve for softening and hardening system, respectively. It seems that for the considered conditions, the response statistical characteristics obtained by the proposed method show reasonable accuracy when compared with the pertinent Monte Carlo estimates. It should be noted that all the preceding numerical examples yield results with reasonable accuracy when the relevant assumptions [32] are not violated.   Fig. 2b. Similary, the proposed method cannot predict the slowly varying nonstationarity of the stochastic response component due to the non-stationarity of the stochastic excitation. Another limitation is that the increase in the harmonic excitation amplitude reduces the accuracy of the proposed method. This is due to the assumption adopted in Eqs. (15a)-(15b) and (28a)-(28d) that the deterministic response μẋ and μ z is smaller compared to the stochastic component standard deviation σ˙x and σˆz , respectively; see also Ref. [32] for more details.

Concluding remarks
A statistical linearization method has been proposed for the response determination of a SDOF hysteretic system endowed with fractional element and subjected to combined periodic and white/colored excitation. This has been achieved by decomposing the total stochastic response into a summation of deterministic and zeromean stochastic component. By doing so, the differential equation governing the total response has been decomposed into two sets of differential equations, governing the deterministic response component and the zero-mean stochastic response component, respectively. Next, the harmonic balance method has been utilized to establish a relationship between the Fourier coefficients of the periodic excitation and those of the deterministic response component, whereas the statistical linearization has been used for the second moments of the stochastic response component. However, the derived nonlinear algebraic equations obtained by the two techniques are coupled, and thus standard numerical iteration method has been adopted to obtain the response quantities numerically and simultaneously. Finally, the pertinent Monte Carlo simulation has been used to demonstrate the applicability and accuracy of the proposed method. Specifically, a softening and hardening Bouc-Wen hysteretic system endowed with fractional element and subjected to combined excitation has been used as numerical examples for illustration. It has been found that the accuracy of the proposed method is comparable with the well-established statistical linearization method for a integer-order Bouc-Wen system subjected to stochastic excitation only. It should be noted that the proposed method can be served as an beneficial attempt toward efficient and versatile solution for the considered problem. Further investigation should focus on the extension of the method for multi-degree-of-freedom fractional-order systems with other type of hysteretic models and subjected to combined periodic and non-stationary or/and non-Gaussian stochastic excitation.