Bayesian Model Identiﬁcation Through Harmonic Balance Method For Hysteresis Prediction in Bolted Joints

Hysteresis is a nonlinear dissipative phenomenon present in many structures, such as those assembled by bolted joints. However, the approaches for identiﬁcation are still limited in the literature due to their complexity. Addition-ally, there are also uncertainties held in bolted structures caused to ﬂuctuations in tightening torque and pressure distribution along the contact surface. Thus, this paper yields a methodology for identifying a stochastic Bouc-Wen model for bolted joints based on the harmonic balance method. Un-fortunately, some challenges are encountered when applying conventional series approximation to hysteresis, caused mainly by non-smooth behavior, which induces abrupt transitions between diﬀerent motion regimes. In this work, previous adaptations were made to split the hysteresis loop in smooth paths and then use a piecewise harmonic balance approach. In this way, it was possible to deal with a deterministic identiﬁcation problem based on minimizing the error between the Fourier amplitudes of an experimental signal and those obtained through harmonic balance applying the Cross-Entropy optimization method. So, the results were extended to a stochastic model by applying the Bayesian paradigm, in which the maximized likelihood function was also based on the harmonic balance amplitudes. This methodology was demonstrated to identify Bouc-Wen parameters capable of predicting hysteresis in the BERT benchmark, composed of two aluminum beams jointed by a bolted joint in a cantilever boundary condition. Evaluating the results in the time and frequency domain and the nonlinear behavior through the hysteresis loop, it can be concluded that the method was able to identify an accurate stochastic Bouc-Wen model in pre-L. dicting the dynamics of bolted structures even taking into account the probable uncertainties of the system.


Introduction
A representative number of assembled structures are connected by bolted joints, which leads these elements to have a substantial presence in the industry. The presence of these joints reflects on the dynamic behavior of the structure, in particular on its global damping and stiffness characteristics, due to nonlinear contact interactions, e.g., friction, that take place at the micro-scale of the joint area. The behavior of these structures is commonly associated with the presence of hysteresis, a nonlinear effect that relates inputs and outputs in a non-smooth way, induces delays, rate-dependent or independent memory effects, and multiple solutions [35]. Moreover, due to slipping effects that occur when the joint is subjected to pronounced vibration amplitudes, nonlinear softening effects take place in the vicinity of resonance frequencies.
There are approaches available in the literature to address the dynamics of jointed structures. Finite element (FE) analysis offers models with high fidelity in describing the contact interface (when hundreds of degrees of freedom are taken into account). However, the computational cost required to simulate these models during a nonlinear regime of motion, especially in obtaining responses in the time domain, is sometimes prohibitive when a dynamic response that spans several oscillation cycles is desired. To reduce the efforts related to numerical simulations of FE models, several works in literature have proposed the use of simplified joint models, such as the Iwan model. This model stands out for its ability to describe typical nonlinearities of bolted joints, such as softening stiffness and hysteresis contact forces, by using a distribution of friction sliders. The popularization of such a model is also directly related to its simplicity of implementation, offering the possibility to calibrate its parameters based on experimental observations. It is argued that any Iwan model can be represented by an equivalent Bouc-Wen model, which is described by a set of nonlinear differential equations and, due to its versatility, it is capable of representing several hysteresis loops in the restoring force × displacement plane just by adjusting its parameters properly. Applications of this model include approximations of magnetorheological fluid dampers [18], modeling of piezoelectric actuators [30], reconstruction of hysteretic forces of seismic isolators [24], and others [14]. Concerning engineering problems involving joints, Oldfield et al. [25] proposed using the Bouc-Wen model to adjust the hysteresis loops yielded by a finite element model of bolted joints. Fantetti et al. [11], in turn, replicated hysteresis loops from a fretting test rig using a modified Bouc-Wen model, which takes wear into account, to simulate the evolution of contact parameters. Teloli et al. [34] put forward the Bouc-Wen model to fit the experimental response of the BERT benchmark, an experimental setup that consists of assembling two beams by a double bolted joint. However, compared to the Iwan model, there is a lack of contributions involving the use of the Bouc-Wen model to characterize contact interactions in bolted joints. One of the reasons behind this fact is related to the difficulty in describing the response of this phenomenological model through analytical expressions since closed-form solutions are a useful tool in parameter estimation and calibration procedures.
Based on the context presented, the paper's main contribution lies in using closed-form solutions to identify parameters that predict the response of bolted structures considering a reduced-order stochastic Bouc-Wen model. The proposed procedure is mainly based on analytically approximating the system's response and its nonlinear restoring force via Fourier series using the harmonic balance method (HBM). The specific literature has already brought different ways to overcome the challenging issues faced when approximating the output of hysteretical systems (which are generally non-smooth systems) by smooth series, as HBM aims. Among them, we have, for example, some works that had done it applying the Galerkin method [37,23] or even an incremental version of HBM [16]. Another way proposed by Cameron and Griffin [6] is the Alternating Frequency-Time (AFT) method, which aims to take advantage of working in each domain inconvenient steps and thus iteratively moves between them. In other works, variations of the AFT method have also been applied to frictional systems [19] or even combined with the HBM [40,39]. A new approach has also been presented by Miguel et al. [21], in which a smoothing procedure was previously carried out to the hysteresis force, then achieving a piecewise HBM, which is employed in the present work. Although the HBM-based procedure for identifying a Bouc-Wen model using this smoothing approach has already been proposed in Miguel et al. [22], this paper goes further and expands the early deterministic model proposed to a stochastic one through the Bayesian paradigm. It was also employed the harmonic amplitudes from the Fourier series as identification metrics, as done in the deterministic methodology. Such model improvement proves to be quite relevant because it predicts responses with statistical confidence for bolted joints, which are known as a great source of uncertainties due to possible variability in the assembly and the contact surface dynamics [3]. Also, an updated connection model with a reduced dimension that admits an analytical solution with statistical confidence could be a feasible alternative to satisfactorily circumvent the high computational cost required by the refined FE meshes currently necessary to obtain the dynamic behavior of bolted structures.
The paper is organized as follows: in addition to the introduction, the following section shows an overview of the smoothing procedure for the Bouc-Wen model and the HBM approach applied to describe the hysteresis effect [21]. Then. The HBM-based procedure is presented in Section 3, which contains two main subsections: first, in Section 3.2, the deterministic identification technique of the Bouc-Wen model proposed by Miguel et al. [22] is presented. The methodology is based on minimizing the error between mono-harmonic responses from the experimental data and the analytical approximation through the Cross-Entropy optimization procedure [8,17,9]. It is also assumed the adjustment of a transient velocity is an auxiliary identification metric to achieve better results. So, in Section 3.3, the procedure to identify a stochastic model using Bayesian inference is described. This step, which brings the main contributions of this work, considers the system's parameters as random variables with their respective probability distributions. Starting from a priori knowledge about the structural dynamics based on the previous deterministic results, a posteriori statistical model employing the Markov Chain Monte Carlo/Metropolis-Hastings (MCMC) algorithm is identified through the maximization of the likelihood of the harmonic amplitudes obtained via HBM. In Section 4 the efficacy of this approach is demonstrated applying it on the BoltEd stRucTure (BERT ) benchmark [22,34], which is composed of two aluminum beams jointed by a symmetric double-bolted joint in a cantilever boundary condition. The preliminary deterministic model is briefly presented to illustrate the preliminary knowledge acquired and then expanded to the stochastic model using the Bayesian paradigm formulation, allowing to generate a model that can accommodate the uncertainties existing in the system through a mean and a confidence band. Also, it is For more information access: https://github.com/ shm-unesp/DATASET_BOLTEDBEAM essential to point out that the BERT benchmark was chosen because it has all the nonlinear behavior expected to bolted structures, besides it had already been statistically identified by Teloli et al. [34] using higher-order FRFs estimated through Volterra series expansion, enabling a methodological comparison that is carried out at the end of this section. Finally, the concluding remarks and trends for future research directions in this topic are presented. Finally, the concluding remarks and trends for future research directions in this topic are presented in Section 5.

The HBM Applied to Hysteretical Systems
Among several methods to approximate output signals in nonlinear systems, harmonic balance is an analytical method in the frequency domain that carries a simple proposal and exact physical meaning. In opposite to what happens in linear systems, when a pure sinusoidal input is applied to a nonlinear system, the output does not follow the input's mono-harmonic frequency but is distorted by the presence of high-order harmonics [38]. Based on this feature, the harmonic balance's fundamental premise is to capture both the fundamental response and the higher-order harmonic terms through the sum of sinusoidal functions with a truncated Fourier series, according to the required precision. A complete discussion with simulated examples of the method, including the particular case of hysteretic systems, can be found in Miguel et al. [21], where the formulation applied to this work is presented. This section aims to summarize the key ideas of the methodology.

General HBM formulation
To briefly illustrate the main ideas of HBM, a general nonlinear mechanical system is considered, whose differential equation of motion is given by: in which , and represent the displacement, velocity and acceleration, respectively; is the mass constant, and are the viscous damping and linear stiffness, respectively; Z( ) is a general representation of the nonlinear restoring force element. For the application of HBM it is assumed that the system is subject to a mono-harmonic sinusoidal excitation ( ) = sin ( ), where [N] is the input amplitude, and the frequency. For this input signal it is assumed a steadystate harmonic output that generates a restoring force with these same features, allowing them to be approximated by an expansion in Fourier series: where is the number of harmonic terms considered in the approximation, and , , A and B are Fourier coefficients of displacement and nonlinear force, respectively. In this second one, the coefficients are conveniently expressed through the integrals of classical Fourier analysis: Replacing the Eqs. 3 to 5 on Eq. 1 and balancing the harmonic terms of the equality, we arrive at a system of nonlinear equations that in general is solved numerically. It is worth noting that in some frequency ranges, the system can present multiple solutions due to physical effects such as the jumping phenomena [4] and even numerical instabilities that can occur mainly when high order harmonic terms are considered. This second case can manifest itself in situations where the method does not reach the convergence or reach a false convergence.

Particularities of HBM applied to hysteresis systems
Unfortunately, this simple and direct application presented earlier is only possible when dealing with smooth nonlinearities, since the restoring force is continuous and does not undergo abrupt transitions between different operating regimes, differently from what occurs in non-smooth systems such as the hysteretic ones. Another factor that limits its use is that the hysteresis forces are often described in differential equations with no single and trivial solution in terms purely related to the output. These features are exemplified here to the Bouc-Wen model, which is the chosen one for modeling applications further in this work. A dynamical system coupled to a hysteretic dissipation element can be modeled as a Bouc-Wen oscillator, whose the nonlinear restoring force follows the first order differential equation of Z( , ) [36]: in which , , and are responsible for determining the shape of the obtained hysteresis loop, and are called Bouc-Wen parameters. As pointed out by Jalali [15], the term Z( , ) does not offer an explicit expansion in terms of the output or even its derivatives and , restricting the applicability of the HBM along with other complications associated with the non-smoothness features of hysteresis models. In order to overcome such issues, a smoothing procedure was developed based on the integration of the differential equation that governs the hysteresis force [33]. For this purpose Eq. 6 considering the case in which = 1 (without loss of generality) is initially divided by , resulting in the following equation which in turn is a differential equation in . The definite integral of the equation at convenient intervals that ensure all the different combinations of the signal function's argument results in four intervals for Z, within which smooth and explicit functions can represent the force in terms of displacement, given by: (ii) path: 0, Z 0 (iii) path: 0, Z 0 (iv) path: 0, Z 0 in which paths (iii) and (iv), characterized by 0, make up the loading cycle, while paths (i) and (ii), characterized by 0, make up the unloading cycle. Additionally, 0 was also defined as the displacement to Z = 0, as shown in Fig. 2.2. Besides that, the smooth equations obtained allow a finite expansion in Taylor series around the term 0 : With the hysteresis force wholly described in the polynomial form, the challenges related to it are solved; however, it is necessary to adapt the Fourier coefficients' integrals. Considering that the restoring force was divided into four smooth intervals alternating every quarter of the excitation period, the integrals of the Fourier analysis, contained in Eqs. 16 and 5, are also split into a quarter of the period, resulting in the following integrals: that replaced in Eq. 3 allow finding an average Fourier series considering all regimes of motion. Some aspects of the displacement signal produced by the Bouc-Wen model must also be assumed, such as the fact that it has symmetric motion between its maximum and minimum points. Due to this symmetrical aspect, the contribution of even-order harmonic and constant displacement terms can be considered nulls [38]. Once the series has been calculated, obtaining and solving the system usually proceeds as described in the general case.

Problem Definition
This paper aims to propose a parameter estimation procedure applied to problems involving bolted joints that admit approximation by hysteresis models. The Bouc-Wen model is used to represent the set of observed data D of the first vibration mode of the BERT beam, a benchmark that contains a fully symmetric double-bolted joint. This data exhibit fluctuation due to environmental variation and uncertainties in the measurement process since the experimental tests were conducted on different days. First, a reference model is calibrated considering deterministic Bouc-Wen parameters. In this step, the cross-entropy (CE) method is used to estimate the set of feasible parameters = * ∈ R n that best fits a single experimental realization D ∈ D, as illustrated by Miguel et al. [22].
The construction of the stochastic Bouc-Wen model is then performed. The calibrated parameter values for the deterministic model are used to propose a priori uniform dis- is the randomized version of the vector . These probability density distributions are then updated to posteriori ones ( |D) by the MCMC algorithm based on the learning inferred from experimental observations.
As discussed by Ikhouane and Rodellar [13], there are several combinations of the Bouc-Wen parameters that could reproduce the same input-output behavior, making the procedure of optimizing its parameters computationally laborious. To circumvent this technical issue and reduce the set of possible solutions S to the problem, this work divides the parameter estimation procedure into two steps: 1. Linear Regime of Motion: This step takes advantage of the knowledge acquired by the underlying physics of the system of interest by estimating its modal parameters. These parameters are calibrated at specific low displacement conditions such that the system is still operating in a linear vibration regime and subject to a broadband power input spectrum with very low excitation amplitude; 2. Nonlinear Regime of Motion: Having the modal parameters to constraint the nonlinear optimization problem, the parameters responsible for the nonlinear dynamics of the Bouc-Wen oscillator are estimated via analytical expressions derived by the HBM-based formulation to approximate experimental measurements from steppedsine tests. Under these controlled periodic excitations, the system of interest bevahes nonlinearly. Figure 2 shows the framework of the parameter estimation procedure proposed by this work.

Deterministic identification procedure
The Bouc-Wen model from Eqs. (1) and (6) is rewritten in the mass normalized form, yielding: where is the linear resonance frequency defined as = √︁˜+˜a nd is the damping ratio; = ˜,˜,˜ is the set of parameters to be estimated through the CE method when the structure is under nonlinear regime of motion.
Firstly, for low vibration amplitude, it is assumed that the structure behaves linearly, in which the pair ( , ) can be adjusted by any classical modal analysis method, or curve fitting strategy, such as the Complex-Exponential Method applied in this work, among others [20].
These parameters are then fixed as constraints on the nonlinear optimization problem, which is formulated to find a set of feasible parameters that minimizes an objective function ∈ S ↦ → R ( ) [22]: where the set of admissible parameters is limited to the interval S = [ min , max ], and R ( ) is the total residue between features inside D built of different vibration data and the predicted ones, yielding: Assuming that the experimental response admits a firstorder harmonic solution considering the complex representation: where exp 1 ( ) and exp 1 ( ) are the harmonic coefficients calculated directly from the Fourier transform exp and its complex conjugate exp measured at a frequency . Note that these harmonic coefficients will be extracted at each frequency increment that composes stepped sine tests. Thus, the HBM terms are given by the distance between the experimental Fourier amplitudes and the analytical ones assuming = 1 harmonic term and = 0,1,2,3 on Taylor series: Fig. 2 Flowchart of the identification procedure. As can be seen, the procedure has two main branches: deterministic and later stochastic identification. Each uses the data acquired and processed and has a linear and nonlinear identification step, as is proposed throughout the work. At the end of the deterministic identification step, the knowledge obtained is used as an initial guess for the iterations of the MCMC/Metropolis-Hastings algorithm that performs the Bayesian identification, in the form of a uniform distribution called prior, to obtain the posteriori distribution of the parameters.
in which ) are the amplitude vectors over the frequency spectrum.
Considering that there is a negligible delay between the sinusoidal input and the hysteresis force, the time instant = 0 when the nonlinear restoring force is zero coincides with the zero input value, one can conclude that the displacement term 0 on Eqs. (8)-(15) results in 0 = exp 1 [22]. This enables one to evaluate the analytical amplitude at each excitation frequency using the HBM-based formulation.
The transient term, in turn, is given by: resulting from the difference between the experimentally measured velocityˆ exp ( ) and its predictive counterpart ( ; ) numerically integrated through the Runge-Kutta 4th order scheme applied on Eqs. (18)- (19). This term is included in the calibration of the reference model to take into account transient effects in system response, since the HBM only considers steady-state response. It will be seen in Section 4 that the velocity responses correspond to data generated by swept sine excitation tests in the vicinity of the first resonant frequency. In Eqs.
Based on the nonlinear optimization problem addressed in Eq. (20), the CE method [29] is considered to find a feasible solution * . For the case of parameter identification, it consists of firstly selecting a probability distribution ( , ) defined by hyperparameters = [ , ], where is the mean and is the standard deviation of the random vector . In this work, the realizations ( = 1, . . . , ) are sampled according to a Gaussian truncated distribution within the S domain.
After sampling the random vector , the objective function is evaluated and sorted in ascending order, such that R (1) ≤ · · · ≤ R ( ) . At this step, an elite sample E is composed of those { 1 , . . . , } for which the performances R (1) ≤ · · · ≤ R ( ) present their smallest values, since the values that minimize the objective function are of utmost importance. The number of candidates that compose an elite sample is defined by selecting a subset of = samples, in which 0 < < 1 is the rarity parameter and it represents the -quantile of performances [17].
Based on the subset E and following the analytical formulas for the mean and standard deviation for truncated Gaussian distributions, the hyperparameters then are updated: The CE method is a sampling technique, which refines the solution candidates at each iteration. Instead of considering the hyperparameter values from equations (26)- (27) to sample in the next iteration, a heuristic smoothing procedure is applied to prevent components of from being zero or one at the first few iterations [17]: where 0 < < 1 is a fixed smoothing parameter that weights the hyperparameter's update on the learning phase.
The main idea of the CE method lies on using the information obtained by evaluating R ( ) for all the independent and identically distributed samples to iteratively drive the prior assumed distribution in direction to the optimal set of parameters (by moving its mean) and to concentrate it around * (by decreasing the deviation). In other words, the method makes: In which ( − * ) is a multivariable Dirac delta [8]. Thus, an usual stopping criterion is do the iterations while ∞ > , in which is a small tolerance related to an ideal Dirac delta, in which = 0, and • ∞ is the L ∞ -norm. Algorithm 1 summarizes the cross-entropy optimization procedure employed [9,17,8].

Bayesian inference in stochastic identification using HBM amplitudes
The Bayesian inference for identification purposes is based on selecting random parameters that maximize a likelihood function (D | ) that best represents the set of observed data [7]. Estimation of this function allows us to update posterior probability density functions (PDF) of the numerical model parameters ( | D) based on the information inferred about the system of interest. Besides, due to Bayesian paradigm, the amount of samples used for model learning is not exhaustive to measure, which makes this technique widely used to identify stochastic models. From the Bayes rule of conditional probability, the posterior PDF is given by [10]: where ( ) is an a prior probability distribution introduced from preliminary assumptions, constraints, or parameters of a deterministic reference model, as done by this work; (D) is the PDF of the observed data that ensures ( | D) is a probability density function with integral equal to unity. Thus, note that if ( ) ∼ U (1 − Δ) * , (1 + Δ) * , then ( | D) ∝ (D | ). The generation of the posterior distribution can be understood as a process of increasing statistical information about the system's response.
The analytical expression of the likelihood function can be derived assuming the following: where D M ( ) is the model M prediction given a set of parameters , ∼ N (0, I 2 ) is an additive decorrelated Gaussian noise ∈ R with covariance I 2 , and I is the identity matrix. Thus, (D | ) results in: To evaluate the likelihood function in this work, the model prediction is taken also as mono-harmonic Fourier amplitudes of a Bouc-Wen model defined by and computed through HBM. Therefore, the Bayesian inference is based on minimizing the residue: where H ( ) is the absolute value vector composed by monoharmonic analytical responses along the frequency range of interest, with elements given by: In a similar way, H (D ) is the absolute value vector of experimental responses from a realization D ∈ D along the same frequency range. Also, the summation indicates that R ( ) is composed by the sum of the residue evaluated in all realizations of D. Thus, the likelihood function can be rewritten as: As can be seen on Eq. 35, as both residues reduces, the likelihood increases. In other words, minimizing the residue yields maximizing the likelihood, which in turn becomes a measure of how well a set of parameters fits the entire dataset. The most representative parameters = * are those that reach the maximum a posteriori probability (MAP): where K denotes the domain constrained by the uniform prior ( ), i.e., K = (1 − Δ) * , (1 + Δ) * . Although the Bayesian paradigm has a solid theoretical basis for evaluating the likelihood function, it does not specify by itself methods for obtaining the posterior efficiently, primarily when the model is defined by a combination of parameters with different distributions sampled simultaneously. Due to this fact, there is an auxiliary algorithm in the literature from the family of Monte Carlo methods that are useful to avoid high-dimensional integration, such as the Markov Chain Monte Carlo/Metropolis-Hastings [12,28], illustrated in Algorithm 2 [27].
From the Algorithm 2, X ∈ R × is the matrix that contains the sets of parameters accepted by the algorithm to compose the posterior distributions. The operator L { } evaluates the likelihood function of Eq. 35 for the set , originated from the interpolation of the vector Y, which each element belongs to the range [0,1], between the boundaries of the domain K. Each sample Y is a candidate to compose the set X, and is generated from a normal distribution with standard deviation adjusted to obtain ∼ 50 % acceptance rate of candidates. Finally, must be chosen to ensure that the model reaches the convergence, evaluated by the function:
which takes into account the mean of the squared euclidean norm of all the lines of the matrix H * ∈ R × , that stores the vector H ( ) to each accepted candidate ( = 1, . . . , ), and which the dimension is increased at each iteration up to . It is important to note that there are other alternative methodologies to reach a stochastic model. A pretty common way is to perform repeatedly deterministic identification procedures to a dataset and evaluate the parameters identified to each one. The main drawback of this approach is that you will need a different realization to each set of parameters identified, increasing the amount of data required. Further, as the procedure will be repeated several times, the identification should be as simple as possible to do not increase the computational cost too much. An example of this kind of methodology was carried out as a preliminary investigation on the linear parameters of this work and can be seen at the beginning of Section 4.2.2. It was feasible because it was performed only through classical modal analysis on white noise inputs.

Experimental setup
The BERT benchmark is based on some experimental setups with bolted joints present in the literature [31,1,16,2,32] to experimentally illustrate the hysteretic behavior in bolted joints dynamics. This benchmark is proposed to test identification frameworks, and due to this reason, it was used in this work to exemplify the application of the methodology described in Section 3. The beam is composed of two aluminum beams whose dimensions are 270 mm x 25.4 mm x 6.35 mm, and are connected by a double-bolted joint with a contact area of 40 mm x 25.4 mm and 5 Nm tightening torque applied to M5 metric hex bolts. The relative motion of micro-slip between the two surfaces in contact induces a global hysteresis behavior that the Bouc-Wen model can fit. Figure 3 shows the complete experimental setup of the BERT beam. Figure 3(b) details the structure and the instrumentation employed: to measure the output acceleration of the system, four accelerometers span the length of the beams; to measure the velocity at the free end of the beam, a laser vibrometer Polytec OFV-525/5000S is used. As this work focuses on the first bending mode, a reduced-order model is constructed, taking only the laser signal into account once it provides better observability of this mode that has more pronounced displacements at the free end of the beam. Besides that, the first resonance range was at relatively low input frequencies (around 5 and 30 Hz), and the laser was capable of measuring low-frequency signals with higher quality than the accelerometers. The tests performed on the test bench for this work include white noise, swept-sine, and stepped-sine excitations, and the application of each one will be presented along with the development of the results in Section 4.  These results indicate that the peak amplitude decreases as the excitation level increases, i.e., the structure does not hold the superposition principle, which is the main feature of a linear system. Moreover, Fig. 4(b) shows the BERT's response to the stepped-sine test performed in a frequency range from 13 up to 23 Hz with an incremental step 0.10 Hz and a sampling rate of 1600 Hz. For this test, only the low and medium amplitudes were considered to ensure structural health. It concentrates more energy at fixed frequencies (mainly at resonance range), generating mono-harmonic inputs that induce large displacements. Note the decreasing of resonance frequencies as the excitation increases, e. g., from 18.5 Hz to 18.3 Hz for the amplitudes of 0.05 V and 0.10 V, respectively. These aspects characterize the occurrence of the so-called softening behavior, which is a nonlinear ef-fect that generates a loss on the structural stiffness in regimes with more significant displacement (such as in the vicinity of resonance for a high level of input) that shifts the resonance frequency to the left.
Both nonlinear mechanisms of amplitude attenuation and loss of stiffness depending on the excitation amplitude are compatible with the presence of hysteresis behavior in structures assembled by bolted joints, which paves the way for the proposition of the Bouc-Wen oscillator, which is a very versatile hysteretic model, for modeling the BERT benchmark. Therefore, its parameters can be calibrated by the proposed optimization procedure in Section 2.

Deterministic identification step
As previously stated in Section 3.2, the first step of the deterministic identification consists in obtaining the modal parameters that best fit the experimental FFR of the first bending mode, applying any classical modal analysis procedure [20]. A white noise test was carried out with a low amplitude level supplied in the shaker amplifier (0.05 V) to achieve this. This is a broadband excitation that, for low input amplitudes, generates responses close to the linear behavior. Figure 5 shows the experimental FRF compared to the linear model identified, which has a resonance frequency of = 18.8 Hz and damping ratio of = 0.4634 %.
To cover a large frequency band with a sinusoidal input (for which the HBM is formulated), a stepped-sine test D ∈ D was performed to the nonlinear identification step. The input was set to excite each frequency for 32 seconds to ensure a steady-state response. The test conditions were the same from Fig. 4. From the measured response, it was possible to estimate experimental Fourier coefficients exp 1 ( ) and exp 1 ( ) for each frequency increment. Then, these values were compared to the respective HBM amplitudes 1 ( ; ) and 1 ( ; ) in function of a sampled candidate to obtain the residues R ( ) and R ∫ ( ) from Eq. 23 and 24. A swept-sine test was also carried out to evaluate the transient behavior through the residue R ( ) from Eq. 25. Finally, having the residues, the objective function R ( ) from Eq. 21 is minimized by the CE method. Table 1 presents the lower and upper limits min and max , respectively, of the feasible solutions S. They were selected through preliminary assumptions. For instance, the observation of a softening behavior reveals the requirement of 0 <˜≤ [ 34].
Hereupon    Table 1 Interval S of feasible solutions considered for the optimization procedure.
To illustrate the validity of the updated parameters, Fig.  6 compares the hysteresis loops on the displacement × nonlinear restoring force plane predicted by the deterministic Bouc-Wen model and the ones obtained from the experi-mental data. Such verification was performed considering the low (0.05 V), medium (0.10 V), and high (0.15 V) levels of excitation amplitude. Figure 6 indicates that the Bouc-Wen restoring force identified is capable of predicting the nonlinear behavior observed in the experimental bolted structure properly. Notwithstanding, note that the energy dissipated by the hysteresis during each oscillation cycle, which is proportional to the loop's internal area, increases according to the input level. These results are consistent with the ones observed in Fig.  4(a), in which the peak amplitudes are presumably reduced  Table 2 Bouc-Wen parameters of the model identified. due to nonlinear damping. The predicted hysteresis loops fit the experimental observations in terms of shape and scale, which indicates that the model proposed has a good agreement in assessing the energy dissipation by the hysteresis effect.

Stochastic identification step
Based on the deterministic model, the stochastic identification procedure is then carried out. Firstly, the modal parameters were identified through frequency response curves estimated from 50 experimental runs considering a white noise excitation at a low input level (0.05 V). From these tests it was preliminary analyzed the coefficient of variation, that is, a percentage that evaluates the dispersion of a random variable in a dimensionless way and is defined as For the natural frequency, it was computed a coefficient of variation = 0.06% and for the damping ratio = 9.37%. Since the resonance frequency does not vary considerably between the experimentally estimated samples, this parameter is fixed at its mean value = 18.8 Hz during the stochastic identification procedure. Due to the higher computational cost and fewer experiments available, this preliminary finding was not feasible for nonlinear parameters.
Thus, the vector of parameters considered in the stochastic identification is: = [˜,˜,˜, ] T .
Defined the set to be identified, another important step that deserves attention is adjusting the MCMC/Metropolis-Hastings algorithm. Based on the inferred information about the deterministic model, and a priori distribution ( ) is proposed. A natural consideration is to define the a priori PDF as a uniform distribution centered on the deterministic parameter set * . However, this first suggestion may not be the most representative when considering the whole D dataset and needs to be updated when appropriate. For example, the posteriori distribution of some of the parameters can tend to concentrate its values near the upper or lower limits of the proposed Uniform distribution, indicating that it may be appropriate to update the Uniform distribution to a new mean value. Thus, after adjusting the Uniform distribution, Tab. 3 presents its maximum and minimum values by considering a Δ = 30 % for each side of a central value.  Table 3 Limits of the uniform prior distribution.
As an estimate of experimental variance 2 , the sum of the main diagonal of the covariance matrix evaluated to D was considered, resulting on 2 = 3.5842 × 10 −8 . It is worth mentioning here that the data considered as experimental was series of synthetic data generated from really experimental stepped-sine signal from the test bench in Fig.  4. Thus, it was considered 10 realizations with medium excitation level generated by simply adding directly on the experimental stepped-sine curve a Gaussian noise with standard deviation = 1.0% of the RMS value of the curve. Lastly, using 2 = 8.6 × 10 −3 adjusted to ensure the acceptance rate ≈ 50%, the convergence of the method were evaluated. Figure 7 presents the convergence function calculated after each sampling.
Note that the convergence is reached around the 1200 th sample. A burn-in of 200 samples was assumed, then only 1000 samples were available to estimate the posterior PDFs for the parameter set , which are depicted in Fig 8. Some statistical aspects about the PDFs can be readily noted: on the one hand, for , and unimodal PDFs were identified,   being much more concentrated around this unique mode for the last two. On the other hand, resulted in a bimodal distri- bution covering all the assumed prior range, indicating quite a difficulty in identifying the linear damping ratio. Despite this, the methodology was able to predict in which regions there are more suitable parameters. . Table 4 shows the mean and MAP values for each PDF from Fig. 8, and also allows comparing them with the first deterministic result * As suggested by the PDFs, Tab. 4 confirms that the final parameters identified considering the entire dataset tended to slightly distance from the deterministic value. However, when considering distributions, we see that in the case of the parameter˜,˜, the value identified initially remains ac-  cepted within the distribution. In opposite, for˜, the values of the admitted parameters were lower, causing the deterministic values to be out of the distribution. These results indicate that when it is intended to adjust the model to represent a broader set of experiments and not just a specific one, the parameters may differ slightly. As the Bouc-Wen model admits more than one set of parameters to reproduce the same physical system, one can conclude that the model depends on specific combinations of them. Because of that, different sets do not necessarily imply any model fault.
To evaluate the results, the following comparisons consider the uncertainty propagation of the Bouc-Wen model generated by series of 1000 Monte-Carlo simulations carried out along the parameter distributions. They were done for swept-sine inputs considering low, medium, and high amplitude with the same test conditions as the responses shown previously in Fig. 4. For comparing the results it was used 10 swept-sine outputs for each excitation level, that were all generated through additive Gaussian noise in experimental signals adjusted to obtain a similar variability than previously computed on Bayesian inference, evaluated by the variance 2 . The main idea was to ensure testing with statistically similar data to learning. Further, it is essential to note that swept-sine data was selected to illustrate the identified model due to the high computational cost in performing Monte Carlo simulations by integrating stepped-sine signals over a long experimental time. From the responses obtained were computed an average response and the limits of 0.995 and 0.005 quantiles, which compose the confidence band of 99 %. Figure 9 shows the comparison between the experimental temporal response and that of the model, with the mean value and confidence intervals with 99 % of the amount responses.
It is noted in the temporal responses that the stochastic model proved to be quite efficient to accommodate the system's response and possible variations, even in transient conditions, once the swept-sine test is primarily a transient input. It is essential to highlight that despite the stochastic model identified through the stepped-sine test, it also fitted the transient behavior of the swept-sine, illustrating the method's robustness. Also, it was possible to reasonably predict responses to low and high input levels; even the model was only identified considering the medium one. Notwith-  standing, it was noticed the presence of relatively worse prediction in the model's output concerning the experimental response just after the resonance region, which results in a large dispersion in the confidence band. This also indicates that it is a region more sensitive to variations in parameters due to the larger envelope when compared to the other areas simulated by the same model and conditions. Figures 10  (b), (d) and (f) reflect the same pattern of uncertainty in the model but regarding the frequency response curves, i.e., a similar effect of more significant uncertainty and prediction difficulty in the region around 20 Hz when analyzing the FRF for the same amplitudes. And finally, Figs. 10 (a), (c), and (e) also illustrates the hysteresis loops for the stochastic model. A notable feature in this stochastic model's loops is that the confidence bands tend to decrease near the plastic regime, i.e., as the nonlinear restoring force goes asymptotically to aa horizontal line opposite to what occurs on the sloping lines of the elastic regime. As most of the uncertainties were inserted in the nonlinear system's parameters, this behavior illustrates the increased influence of hysteresis in the response as the system reaches more significant displacements, accentuating the nonlinearity. However, the model was able to predict the experimental response for all the input levels with quite an accuracy, in particular to the medium amplitude used to identify it. Thus, it can be stated that the model is helpful in predicting the energy dissipated by the hysteresis loop with a reasonable confidence band, demonstrating the adequate functioning of the methodology.
Although the Bouc-Wen model equations considered here to describe the BERT benchmark's response are the same as those considered by Teloli et al. [34], the set of experimental realizations used for each work is different. However, a comparative analysis between both identification frameworks can be drawn, as summarized on the Tab. 5.  Considering the deterministic procedure, there are no considerable differences in the experimental tests needed once both methodologies use white noise, transient sweptsine, and sinusoidal signals in a frequency range; however, while the method based on the higher-order FRFs estimated through the Volterra series expansion considers the sequential quadratic programming as an optimization algorithm, this work uses the CE-method. Regarding the stochastic identification, the framework proposed by Teloli et al. [34] identifies, firstly, the parameters˜and˜considering sine inputs with excitation frequency around one-third of the resonance frequency ( /3) to ensure a closed hysteresis loop, in which there is no dependence of the parameter˜. Depending on the mode of interest, this excitation condition may be close to another vibrating mode, causing undesirable coupling effects. In this sense, the advantage of the HBM-based methodology lies in the possibility of isolating nonlinear modes through stepped-sine tests around the resonance frequency. Still, in stochastic identification, the HBM requires fewer signal classes than the Volterra series formulation, representing a smaller number of experimental tests. Also, the theoretical basis of HBM has its formulation based basically on the Fourier series, which are much simpler than the multidimensional Fourier transforms present on the Volterra Series formulation. However, both have a known physical meaning [26]. Nevertheless, it is important to highlight that both methodologies are helpful and present promising results to identify systems with a hysteresis.

Final Remarks
This work presented a new methodology to identify a stochastic Bouc-Wen model capable of predicting the hysteretic behavior in bolted joint structures. The identification procedure was applied in a BERT benchmark with two beams joined by a bolted connection. Preliminary tests showed that there was the presence of nonlinear behavior induced by frictional effects in the contact surface between them.
The paper brought all the theoretical and practical formulation of the procedure, including previous steps to adjust a deterministic first guess for the identification algorithm. A Bayesian inference so updated the stochastic model over displacement amplitude of the fundamental harmonic component obtained in the function of the Bouc-Wen parameters through a piecewise HBM approach recently presented in the literature [21]. The methodology was proved to be quite efficient in represent the nonlinear behavior of the structure statistically, showing how practical and robust the proper formulation of a simple analytical tool such as HBM can be to identification purposes, even when dealing with more complex nonlinearities than those commonly identified through it in the literature, such as exemplified in the review made by Busseta et al. [5].
Finally, it is also important to highlight that the procedure is in line with current trends in the specific literature to identify hysteretic systems, as the work presented by Teloli et  al. [34], over which this present work even has some methodological advantages. As the HBM approach used is capable to adequately predict the response in function of all the Bouc-Wen parameters without any additional assumptions about the nonlinear behavior presented by the hysteresis loop. It can update all the nonlinear parameters in the same Bayesian step, making the methodology more straightforward to implement and requires a dataset with less variability from different experimental tests.
The presence of methods on the literature that allow not only analytically approximate hysteresis by series of smooth functions but also propose its use as features to adjust hysteretic models to experimental data paves the way for future works on how to overcome challenging issues in the simulation bolted connections of more complex structures that nowadays need FE models with a large number of degrees of freedom and high computational cost to represent the frictional effects on them. Following the contributions presented by this work, it would be possible to adjust a model such as Bouc-Wen with statistical confidence to the transmissibility through the joint and even then use the HBM approach to calculate the response from it, avoiding costly numerical integrations.