Bayesian model identification through harmonic balance method for hysteresis prediction in bolted joints

Hysteresis is a nonlinear phenomenon present in many structures, such as those assembled by bolted joints. Despite a large number of recent findings related to identification techniques for these systems, the problem is still challenging and open to contributions. To fill the gap concerning the proposition of identification algorithms based on closed-form solutions, this work introduces the use of the harmonic balance method to identify a stochastic Bouc–Wen model for predicting the nonlinear behavior of bolted structures. A piecewise smooth procedure is applied on the hysteretic restoring force to become possible to derive an analytical approximation of the response based on the Fourier series. Firstly, the analytical approximation is used to calibrate deterministic Bouc–Wen parameters by minimizing the error between the Fourier amplitudes of the numerical model and those extracted from experimental data using the cross-entropy optimization method. Since the experimental data investigated here contain variability due to the measurement process (aleatoric uncertainties), the deterministic parameters are then used as a priori conditions to update their probability density functions via the Bayesian inference. Having the model parameters as random variables, the stochastic Bouc–Wen model is obtained. This methodology was illustrated in a bolted structure benchmark. The results indicate that the method proposed can identify an accurate stochastic Bouc–Wen model for predicting the dynamics of bolted structures, even taking into account data variability.


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 action of these joints exhibits the dynamic behavior of the structure, particularly 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 operation of these structures is commonly associated with the presence of hysteresis, a nonlinear effect that relates inputs and outputs in a nonsmooth way, induces delays, rate-dependent or independent memory effects, and multiple solutions [1]. Moreover, due to slipping effects that occur when the joint is subjected to pronounced vibration amplitudes, nonlinear softening impacts occur in the vicinity of resonance frequencies.
There are some approaches in the literature to address the dynamics of jointed structures. Finite element (FE) analysis offers high-fidelity models in which hundreds of degrees of freedom are taken into account to describe the contact interface. However, the computational cost required to simulate these models during a nonlinear motion regime is sometimes prohibitive, especially in obtaining responses in the time domain. In this scenario and for applications where modeling the contact area is not of utmost importance, lumped models that take the form of single-degree-of-freedom oscillators driven by hysteretic terms stand out as a solution to reconstruct global dynamics of bolted structures. Reduced-order representations of these structures based on the Iwan model have been notably addressed in the literature (see, for instance, Mathis et al. [2]). This model can describe typical nonlinearities of bolted joints [3], such as softening stiffness and hysteresis contact forces, by using a distribution of friction sliders [4,5]. The popularization of such a model in the research community that investigates jointed structures is also directly related to its simplicity of implementation, offering the possibility to calibrate its parameters based on experimental observations.
Hysteretic terms can also be represented by the Bouc-Wen model. Unlike the Iwan model that considers sliders (Jenkins elements) associated in series or parallel, the Bouc-Wen model assumes a constitutive hysteresis relation in a differential form to capture the evolution of the restoring force. Due to its general mathematical properties, the model is versatile and can describe the hysteresis cycles of a wide variety of applications, ranging from magnetorheological fluid dampers [6], and piezoelectric actuators [7], to the reconstruction of hysteretic forces of seismic isolators [8]. Notwithstanding, it is argued that any Iwan model can be represented by an equivalent Bouc-Wen model; in the context of modeling bolted structures, there exist only a few works that further explore the use of the latter model to this end. Oldfield et al. [9] proposed a simplified Bouc-Wen model to adjust the hysteresis loops yielded by contact interactions from an FE-based model of a structure assembled by bolts. More recently, Teloli et al. [10] put forward an equivalent Bouc-Wen model to fit whole nonlinear restoring force that actuates on the experimental response of the BoltEd stRucTure (BERT) benchmark, which is composed of two aluminum beams jointed by a sym-metric double-bolted joint in a cantilever boundary condition.
However, there is still a lack of contributions related to the calibration process of Bouc-Wen parameters from experimental data of jointed structures since the model constitutive equations restrict the use of conventional nonlinear identification techniques. For example, the non-smoothness and multiple solutions present in the restoring force limit the proposition of solutions for the response through series expansion if no smoothing procedure is employed. Note that closed-form solutions through analytical approximations are helpful for parameter calibration purposes. Recently, to avoid this technical limitation and render it possible to expand the response of the Bouc-Wen model using the Fourier series, Miguel et al. [11] have presented an alternative way to obtain these analytical solutions by dividing the non-smooth restoring force into smooth paths and applying a piecewise approach of the harmonic balance method (HBM) to estimate the response's Fourier coefficients.
It is worth noting that the HBM has already been used to numerically calculate solutions to non-smooth problems, including approximations of the Bouc-Wen model response. In this context of applications, the HBM framework is used as a numeric solver of the system response by using continuation methods (prediction step followed by a correction step) that can calculate the progress of the solution given a frequency interval. This strategy does not necessarily depend on prior knowledge of the constitutive equations of the system under analysis, but only on the estimation of the Fourier coefficients of its response, either in the time or frequency domain. Thus, for parameter calibration purposes, analytical approximations such as the one proposed by Miguel et al. [11] are more appropriate (see [12], in which the analytical HBM framework was used to calibrate parameters from a Duffing-van der Pol oscillator).
Toward this background, this paper's main contribution lies in using closed-form solutions to identify parameters that predict the response of uncertain bolted structures considering a reduced-order stochastic Bouc-Wen model. The identification procedure is mainly based on analytically approximating the system's response and its nonlinear restoring force via Fourier series using the piecewise HBM approach proposed by Miguel et al. [11]. Although this HBM framework has already proved to be helpful to identification purposes on bolted structures by Miguel et al. [13], their work was limited to the assumption of a deterministic system. Therefore, this paper goes further and expands the early deterministic model to a stochastic one through the Bayesian paradigm [14]. Such model improvement proves to be quite relevant since joints are well known for being a great source of uncertainties due to possible variability in the assembly and the contact surface dynamics [15]. Also, an updated hysteresis model with a reduced dimension that admits analytical solution and statistical confidence could be a feasible alternative to satisfactorily circumvent the high computational cost required when assuming a representation based on matrices or meshes extracted from finite element geometries.
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 [11]. Then, the identification procedure is presented in Sect. 3. Firstly, the section addresses the deterministic identification technique that is based on minimizing the error between Fourier amplitudes from the experimental data and the analytical approximation through the cross-entropy (CE) optimization method, which is a global search algorithm based on a sampling technique from the family of Monte Carlo methods [16][17][18]. Then, the section goes through the procedure to identify a stochastic Bouc-Wen model using the Bayesian inference. This step, which brings the main contributions of this work, considers a priori information of the model parameters acquired through the CE framework to update a posteriori probability density functions of the Bouc-Wen parameters by using the Markov chain Monte Carlo/Metropolis-Hastings (MCMC) algorithm [19][20][21]. In Sect. 4, the effectiveness of this approach is demonstrated in an experimental application involving the BERT 1 benchmark [10,13]. The BERT benchmark was chosen because it exhibits the expected nonlinear behavior in bolted structures (e.g., softening effects and hysteresis behavior). Finally, the concluding remarks and trends for future research directions in this topic are presented in Sect. 5.

The HBM applied to hysteretic systems
Harmonic balance is an analytical method in the frequency domain that carries a simple proposal and exact physical meaning among several ways to approximate output signals in nonlinear systems. In opposite to what happens in linear systems, when a nonlinear system is excited by a mono-harmonic input signal, its output does not follow the input's mono-harmonic frequency but is distorted by the presence of high-order harmonics [22]. 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. [11], 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 y,ẏ andÿ represent the displacement, velocity and acceleration, respectively; m is the mass constant, c and k are the viscous damping and linear stiffness, respectively; Z(t) 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 is the input amplitude, and ω the frequency. For this input signal, it is assumed to be a steady-state harmonic output that generates a restoring force with these same features, allowing them to be approximated by an expansion in Fourier series: [a n cos(nωt) + b n sin(nωt)] , where κ is the number of harmonic terms considered in the approximation, and a n , b n , A n and B n 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 Eqs. 3 to 5 in 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 [23] 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.
Further details related to this general series representation can be seen in Krack and Gross [24], in which the authors provide a comprehensive textbook that covers from the theory of harmonic balance and Fourier series to practical aspects of their computational implementation.

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 chosen for further modeling applications in this work. A dynamical system coupled to a hysteretic dissipation element can be modeled as a Bouc-Wen oscillator, whose the nonlinear force term Z(y,ẏ) is obtained by solving the first-order differential equation ofŻ(y,ẏ) [25]: in which α, γ , δ and ν are called Bouc-Wen parameters. They are responsible for inducing and controlling the memory effects and elastoplastic behavior of the Bouc-Wen model. The parameter α is a proportional elastic term that controls the amplitude and the slope of the linear paths of the nonlinear restoring force (paths with a constant slope in the hysteresis loop; for illustrative purposes, in Fig. 1 they correspond to the paths B→A and C→D). On the other hand, δ and γ control the hysteretic relation between y and Z as well as the shape of the plastic regime. Therefore, together they are responsible for determining, for instance, if the system acts with a hardening or softening nonlinear effect. Lastly, ν determines the smoothness of the elastic-plastic transition [26,27]. As pointed out by Jalali [28], the term Z(y,ẏ) does not offer an explicit expansion in terms of the output y or even its derivativesẏ andÿ, restricting the applicability of the HBM along with other complications associated with the non-smoothness features of hysteresis models. To overcome such issues, a smoothing procedure was developed based on the integration of the differential equation that governs the hysteresis force [29]. 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 y. 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: (i) path:ẏ 0, Z 0 (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, y 0 was also defined as the displacement to Z = 0, as shown in Fig. 1. Besides that, the smooth equations obtained allow a finite expansion in Taylor series around the term y 0 , that is: Thus, Eqs. 8 to 11 become: With the hysteresis force wholly described in the polynomial form, its challenges 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. 17 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 [22]. Once the series has been calculated, obtaining and solving the system usually proceed as described in the general case.

Problem definition
This paper proposes 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. These 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 l ∈ D, as illustrated by Miguel et al. [13].
The construction of the stochastic Bouc-Wen model is then performed. The calibrated parameter values for the deterministic model are used to propose a prior uniform distributions π( ) ∼ U (1 − )θ * , (1 + )θ * , with ∈ R, and is the randomized version of the vector θ . These probability density distributions are then updated to posterior ones π( |D) by the MCMC algorithm based on the learning inferred from experimental observations.
As discussed by Ikhouane and Rodellar [30], 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 stepped-sine tests. Under these controlled periodic excitations, the system of interest behaves nonlinearly. Figure 2 shows the framework of the parameter estimation procedure proposed by this work. Note that the interest here lies in approximating the dynamics of the first vibrational mode of an experimental test-bench, such that a reduced-order model described by a single-degree-of-freedom (SDOF) system is suitable. For applications involving multi-degrees of freedom (MDOF), there are no theoretical impediments in applying the proposed method since it is possible to consider hysteretic components on a MDOF representation and, in addition, the HBM has no limitations in this regard either. However, it is worth noting that to numerically solve this analytical approach, this work considers the Newton-Raphson algorithm [31]. In an application involving MDOF systems, the number of unknowns would be substantially more significant, resulting in a problem of increased complexity and, consequently, increased computational effort. There are other methods available in the literature that serve as alternatives to the Newton-Raphson algorithm, such as the fixed-point algorithm, which does not require the calculation of the Jacobian matrix to solve the system of equations [32,33].

Deterministic identification procedure
The Bouc-Wen model from Eqs. (1) and (6) is rewritten in the mass normalized form, yielding: where ω n is the linear resonance frequency defined as ω n = α +k and ζ is the damping ratio; θ = α,γ ,δ is the set of parameters to be estimated 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 prior to obtain the posterior distribution of the parameters 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 (ω n , ζ ) can be identified through any linear modal analysis approach, such as the complex-exponential method applied in this work or other classical methods [34].
These parameters are then fixed as constraints on the nonlinear optimization problem, which aims to obtain a set of feasible parameters θ * constrained to the inter-val S = [θ min , θ max ] which minimizes an objective function θ ∈ S → R(θ ) [13]: where R(θ) is the total residue between features inside D l built of different vibration data and the predicted ones, yielding: HBM terms Assuming that the experimental response admits a first-order harmonic solution considering the complex representation: where a exp 1 (ω i ) and b exp 1 (ω i ) are the harmonic coefficients calculated directly from the Fourier transform Y exp and its complex conjugate Y exp measured at a frequency ω i . 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 harmonic amplitudes obtained experimentally and analytically assuming κ = 1 harmonic term and n = 0, 1, 2, 3 on Taylor series: in which a 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 t = t 0 when the nonlinear restoring force is zero coincides with the zero input value, one can conclude that the displacement term y 0 in Eqs. (8)-(16) results in y 0 = a exp 1 [13]. 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 (t) and its predictive counterpartẏ(θ ; t) numerically integrated by the fourthorder Runge-Kutta method [35] applied in Eqs. (19)- (20). As the HBM only considers steady-state behavior of the system response, this term is included in the calibration to also evaluate the contribution of the transient effects. It will be seen in Sect. 4 that these velocity terms are responses obtained through swept-sine tests around the frequency range of the first vibration mode. In Eqs.
To solve the nonlinear optimization problem addressed in Eq. (21), the CE method [36] is considered to find a feasible solution θ * . The method is based on evaluating sets of parameters sampled from a prior distribution and then updating it iteratively in order to narrow it down, such that it tends to a Dirac delta function around a deterministic result. Therefore, it consists of firstly selecting a probability distribution π( , v) defined by hyperparameters v = [μ, σ ], where μ is the mean and σ is the standard deviation of the random vector . In this work, the realizations θ k (k = 1, . . . , N s ) 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 (N s ) . At this step, an elite sample E i is composed of those {θ 1 , . . . , θ N e } for which the performances R (1) ≤ · · · ≤ R (N e ) present their N e 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 N e = ρ N s samples, in which 0 < ρ < 1 is the rarity parameter and it represents the ρ-quantile of performances [18].
Based on the subset E i 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 Eqs. (27)- (28) to sample θ k in the next iteration, a heuristic smoothing procedure is applied to weights a new update with the history of previous distributions. It prevents a fast decrease in the standard deviation, making the distribution stuck in a region far from the optimal solution. [16,18]: 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 in 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 [16]. Thus, an usual stopping criterion is to do the iterations while σ i ∞ > σ s , in which σ s is a small tolerance related to an ideal Dirac delta, in which σ s = 0, and • ∞ is the L ∞ -norm. Algorithm 1 summarizes the cross-entropy optimization procedure employed [16][17][18].

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 [14]. 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 the Bayesian paradigm, the amount of samples used for model learning is not exhaustive to measure, making this technique widely used to identify stochastic models. From the Bayes rule of conditional probability, the posterior PDF is given by [37]: where π( ) is an 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 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 n 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 mono-harmonic analytical responses along the frequency range of interest, with elements given by: Similarly, H(D l ) is the absolute value vector of experimental responses from a realization D l ∈ D along the same frequency range. Also, the summation indicates that R B ( ) is composed by the sum of the residues evaluated in all realizations of D. Thus, the likelihood function can be rewritten as: As can be seen in Eq. 36, as both residues reduce, the likelihood increases. In other words, minimizing the residue yields maximizing the likelihood, which measures how well a set of parameters fits the entire dataset. The most representative parameters θ = θ * are those that reach the maximum a posteriori probability (MAP): * = argmax 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 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 [19,21], as illustrated in Algorithm 2 [20].

Algorithm 2:
Markov Chain Monte Carlo/Metropolis-Hastings for parameter identification. 1 Set the counter i = 1, a first candidate X 0 = {0.5, 0.5, 0.5, 0.5} T , and define σ s and i s ; 2 Generate a candidate Y ∼ N (X i−1 , σ s ); 3 Interpolate the value Else, reject candidate X i = X i−1 . Set i = i + 1; 6 Do the steps 2 to 5 while i < i s . From Algorithm 2, X ∈ R i s ×n is the matrix that contains the i s sets of n parameters accepted by the algorithm to compose the posterior distributions. The operator L {θ } evaluates the likelihood function of Eq. 36 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 σ s adjusted to obtain ∼ 50 % acceptance rate of candidates. Finally, i s 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 i lines of the matrix H * ∈ R i×N ω , that stores the vector H( k ) to each accepted candidate k (k = 1, . . . , i), and which the dimension is increased at each iteration i up to i s . 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 not increase the computational cost too much. An example of this kind of methodology was carried out as a preliminary investigation of the linear parameters of this work. It can be seen at the beginning of Sect. 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 [38][39][40][41][42] 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 Sect. 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. The relative motion of micro-slip between the two surfaces in contact induces a global hysteresis behavior that the Bouc-Wen model can fit. Two M5 bolts and nuts connect the bolted joint with a tightening torque of 5 Nm assessed by a torque wrench after each experimental run. The clamped end, in turn, was made by two blocks (a fixed and a mobile one for adjustment purposes) that clamped the beam between them as two M8 bolts were tightened. The fixed block was attached in a robust structure by four M8 bolts and washers. The entire assembly was then directly fastened to the inertial table by two other M8 bolts. Before tightening the beam, the parallel alignment between the inertial table and the beam was checked by measuring the distance along its length. It is essential to mention that the torque applied to the clamping bolts was not measured, but the tightening was conducted up to the maximum torque allowed by the bolt body. Moreover, all the connecting components and the clamping structure were made of steel and cast iron, respectively. Figure 3 shows the complete experimental setup of the BERT beam. Figure 3b 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. Since the interest here lies in identifying a reduced-order model, it is considered only the first vibration mode. As its modal shape has more pronounced displacement amplitudes at the free end, the model is constructed based on the signal from the laser vibrometer. 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 Sect. 4.2. A Modal Shop 2400E electrodynamic shaker provided all the input signals with an integrated power amplifier coupled by threaded nylon rod stinger at 85 mm from the clamped end to minimize shaker/structure interactions. The shaker's instruction manual indicates that its resonance frequencies are above 9.5 kHz, without a specific value and depending on the interactions with the test structure. An LMS SCADAS acquisition board then acquires the signals, and as input signals were considered, the voltage supplied by the shaker amplifier converted in [N/kg] by a conversion constant previously identified A = 49.507 [NV −1 /kg]. Figure 4 depicts preliminary tests that justify using an equivalent hysteresis model to represent the BERT system and detect the nonlinear effects acting on the first bending mode. Figure 4a illustrates the frequency response curves of the BERT benchmark estimated from swept-sine tests with frequency linearly increased from 0 up to 40 Hz along 8 seconds considering different low input levels (0.05 V), medium (0.15 V), and high (0.25 V) amplitude and a sampling rate of 1024 Hz. 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. 4b 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 effect 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 Sect. 2.

Deterministic identification step
As previously stated in Sect. 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 [34]. To achieve this, a white noise test was carried out with a low amplitude level supplied in the shaker amplifier (0.05 V). This is a broadband excitation that, for low input amplitudes, generates responses close to the linear behavior. Figure 5    tified, which has a resonance frequency of ω n = 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 steppedsine test D l ∈ 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 a exp 1 (ω i ) and b exp 1 (ω i ) for each frequency ω i increment. Then, these values were compared to the respective HBM amplitudes a 1 (θ k ; ω i ) and b 1 (θ k ; ω i ) in function of a sampled candidate θ k to obtain the residues R c (θ ) and R s (θ) from Eqs. 24 and 25 . A swept-sine test was also carried out to evaluate the transient behavior through the residue R T (θ) from Eq. 26. Finally, having the residues, the objective function R(θ ) from Eq. 22 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 <γ ≤δ [10].
To illustrate the validity of the updated parameters, Fig. 6 compares the hysteresis loops on the displace- ment × nonlinear restoring force plane predicted by the deterministic Bouc-Wen model and the ones obtained from the experimental 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 loop's internal area is quite similar and increases according to the input level. As this area is proportional to the energy dissipation during each oscillation cycle, one can conclude that there is a good prediction of the hysteretic dissipative effects. These results are consistent with the ones observed in Fig. 4a, in which the peak amplitudes are presumably reduced 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, the coefficient of variation was preliminary analyzed, that is, a percentage that evaluates the dispersion of a random variable in a dimensionless way and is defined as: For the natural frequency, a coefficient of variation C V ω n = 0.06% was computed and for the damping ratio C V ζ = 9.37%. Since the resonance frequency does not vary considerably between the experimentally estimated samples, this parameter is fixed at its mean (c) High amplitude (0.20 V). value μ ω n = 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, a prior distribution π( ) is proposed. A natural consideration is to define this 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 posterior 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  Table 3 presents its maximum and minimum values by considering a = 30 % for each side of a central value.
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, 10 realizations were considered with a medium excitation level generated by simply adding a Gaussian noise directly on the experimental stepped-sine curve with standard deviation σ D = 1.0% of the RMS value of the curve.
Lastly, using σ s 2 = 8.6×10 −3 adjusted to ensure the acceptance rate ≈ 50%, the convergence of the method was 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 distribution covering all the assumed prior ranges, indicating difficulty in identifying the linear damping ratio. Despite this, the methodology was able to predict in which regions there are more suitable parameters. It is important to highlight that, by definition, the total area below the PDF represents the probability of occurrence of any of the values from the distribution. 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, Table 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 accepted 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 the parameters may differ slightly when it is intended to adjust the model to represent a broader set of experiments and not just a specific one. 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 their specific combinations. 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 amplitudes with the same test conditions as the responses shown previously in Fig. 4. For comparing the results, 10 swept-sine outputs were used for each excitation level, which 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 were 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. The responses obtained were computed as 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. Notwithstanding, the presence of relatively worse prediction in the model's output concerning the experimental response just after the resonance region was noticed, 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. Figure 10b, d and f reflects 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, Fig. 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 a 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, particularly the medium amplitude used to identify it. Thus, it can be stated that the model helps predict 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. [10], the set of experimental realizations used for each work is different. However, a comparative analysis between both identification frameworks can be drawn, as summarized in Table 5.
Considering the deterministic procedure, there are no considerable differences in the experimental tests needed once both methodologies use white noise, transient swept-sine, and sinusoidal signals in a frequency range; however, while the method based on the higherorder FRFs estimated through the Volterra series expansion considers the sequential quadratic programming as an optimization algorithm, this work uses the CEmethod. Regarding the stochastic identification, the framework proposed by Teloli et al. [10] identifies, firstly, the parametersα andδ considering sine inputs with excitation frequency around one-third of the resonance frequency (ω n /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 fre-quency. Still, the HBM requires fewer signal classes than the Volterra series formulation in stochastic identification, 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 [43]. In terms of time efficiency, the HBM is slightly hampered by the complexity of the nonlinear system of equations generated by the harmonic balance applied to the Bouc-Wen model. As it is numerically solved (in this case by the Newton-Raphson algorithm), this aspect may slow convergence. However, it is not prohibitive by itself; it just becomes costly due to being repeated over the entire frequency range in a context involving Monte Carlo simulations. Compared to the identification procedure based on the Volterra series, the HBM is slower when considered a single iteration of the identification scheme. Nevertheless, as the Volterra-based methodology is performed in two different steps (that also have two other experimental requirements), the time efficiency of both procedures becomes quite similar. Finally, besides the methodologies compared have some practical differences, it is important to highlight that both 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 the 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 [11]. As can be seen in Figs. 9 and 10 , the methodology was proved to be quite efficient in representing the nonlin-   ear effects of the hysteretic restoring force and reproduced the behavior of the structure with statistical confidence both in the time and frequency domains. These results show 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. [44].
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. [10], 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.
Finally, 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 of 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.