BAYESIAN DESIGN OF EXPERIMENT FOR THE ESTIMATION OF VISCOELASTIC PARAMETERS PROJETO DE EXPERIMENTO BAYESIANO PARA A ESTIMAÇÃO DE PARÂMETROS VISCOELÁSTICOS

The present work has as its main objective the formulation and solution of inverse problems aiming at the estimation of viscoelastic parameters of a discrete mechanical system. Its application potential covers several areas in the most varied segments of engineering and industry. For the formulation and solution of the direct problem, a discrete mechanical system with viscoelastic damping is considered. The parameter estimation problem is then formulated as an inverse problem, whose objective is to estimate the viscoelastic properties of a discrete system using a Bayesian approach. First, a Bayesian design of experiment its carried out in order to identify optimal experimental conditions for solving the inverse problem, taking into account the positioning of the sensors and actuators. To solve the parameter estimation inverse problem, an Adaptive Monte Carlo Markov Chain Method is employed.


INTRODUCTION
Mechanical structures are subject to dynamic loading, which, under certain conditions, can generete excessive levels of vibration and may even cause structural failures, thus causing great damage to society. Therefore, it is essential to develop techniques for the construction of mechanical systems with adequate vibration control and isolation, thus allowing their applications in various engineering areas, for example, in the development of viscoelastic dampers used to control vibrations in machines, aeronautical structures and energy dissipation mechanisms in buildings, protecting them against earthquakes damage. These measures have a direct impact on reducing operating costs and maintaining adequate structural functioning.
In this context, the present work proposes the application of the Bayesian approach for the formulation and solution of the inverse problem of estimation of viscoelastic parameters of discrete mechanical systems. Due to their great applicability and efficiency, Bayesian methods are, increasingly, being used by the scientific community to solve problems in several areas of interest (MALAKOFF, 1999). In addition to allowing a priori information to be incorporated into the model, the ease of incorporating it into a formal decision context, the explicit treatment of problem uncertainties and the ability to assimilate new information in adaptive contexts are some of its advantages (COSTA, 2004). In this approach, the variables of problem are modeled as random variables and Probability Density Functions (PDF) are used to incorporate, in the estimation process, prior information about the parameters to be estimated. At the end of the process, an approximation of the a posteriori probability density function of the parameters of interest is obtained (TEIXEIRA et al., 2016).
Sampling techniques based on MCMC methods are usually employed to estimate the posteriori probability density of the parameters of interest. Although conventional MCMC methods are robust and effective in solving inverse problems, they can present some difficulties in properly exploring the search space of the parameters, thus making it necessary to use adaptive techniques, since such adaptive techniques are able to exploit in a way more efficient the search space for the parameters, based on careful adjustments of intrinsic parameters to the Metropolis-Hastings algorithm (TEIXEIRA, 2018). The inverse problem of viscoelastic parameters estimation is then based on changes in the system's impulsive response and its solution will be obtained through Bayesian Inference, in which it will be employed the Adaptive Markov Chain Monte Carlo method, proposed by Haario, Sakaman e Tamminen (2001) and implemented using the Metropolis-Hastings algorithm.

DIRECT PROBLEM FORMULATION
According to the theory of linear viscoelasticity (CHRISTENSEN, 1982), the constitutive model for a linear, isotropic viscoelastic material under constant temperature is given by where σ and ε are, respectively, the strain and normal deformation of the material and E rel represents the relaxation module. Different relaxation modules are proposed in the specialized literature (VASQUES; MOREIRA; RODRIGUES, 2010). In the present work, the relaxation module adopted is given by where E ∞ is the static relaxation module of the material, n i is the number of internal dissipation variables considered in the description of viscoelastic damping, Ω k is the inverse of the relaxation time of the material, under constant deformation and ∆ k represents the corresponding resistance to relaxation. For simplicity, a viscoelastic model with a single internal dissipation variable will be adopted here. Therefore, from now on, the subscript k will be omitted.
Considering Eqs. (1) and (2), and the use of a single internal variable, based on the Laplace transform, the viscoelastic constitutive model for the material can be described as where ξ represents the internal dissipation variable.
For the numerical analysis of the presented methodology, it is considered a problem of estimation of viscoelastic parameters of a simple discrete system of two degrees of freedom, illustrated in Figure (1). This system has masses m 1 and m 2 , elastic springs with stiffness constants k 1 and k 2 , viscous dampers with constants d 1 and d 2 and viscoelastic dampers with constants E ∞1 , ∆ 1 , Ω 1 and E ∞2 , ∆ 2 , Ω 2 . The behavior of the system is described by the coordinates y 1 and y 2 , and it is subjected to external loadings f 1 and f 2 . S1 is the sensor positioned at coordinate y 1 and S2 is the sensor positioned at coordinate y 2 , actuator A1 corresponds to the excitation of the structure at coordinate y 1 and actuator A2 corresponds to the excitation of the structure at coordinate y 2 .
Analogously to obtaining the viscoelastic differential constitutive model, given in Eq.
(3), the reaction force in a viscoelastic damper can be obtained as a function of the relative displacement of its ends. Thus, the equation of motion of the system shown in Figure ( where y is the vector with the displacements of system y 1 and y 2 , ξ is the vector containing the internal variables ξ 1 and ξ 2 , M, D and K are, respectively, the mass, damping and stiffness matrices of the purely elastic system and with viscous damping, K v is the stiffness matrix associated with viscoelastic dampers, the matrix B ζ is the input matrix of the internal variables, H 1 and H 2 are, respectively, the state and input matrices in the equation of evolution of the internal variables. The matrices K v , B ζ , H 1 and H 2 are given by The system solution given in Eq. (4) is obtained through the computational routines ss and impulse (available on the Matlab platform) that calculate the state model and the system's Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 226-01, 226-17, 2020 DOI: 10.21575/25254782rmetg2020vol5n21163 impulsive response, respectively.

Bayesian Design of Experiment
Previously to the application of inverse techniques to estimate the parameters of interest, it is recommended to carry out an experimental optimal design, in order to estimate the best possible conditions for conducting the numerical and experimental simulations. Previously to the application of inverse techniques to estimate the parameters of interest, it is recommended to carry out an experimental optimal design, in order to estimate the best possible conditions for conducting the numerical and experimental simulations. In this work, the Bayesian experimental optimal design is performed, which consists on determining an optimal experimental arrangement d * that allows to efficiently estimate the parameters of interest. For that purpose, metrics are used, being therefore possible to determine, for example, the quantity and the optimal positioning of sensors and actuators, among other sets of factors that lead to an optimal estimate of the parameters of interest (CAPELLARI; CHATZI; MARIANI, 2012; HUAN; MARZOUK, 2013; PAPADIMITRIOU; ARGYRIS, 2017; TSILIFIS; GHANEM; HAJALI, 2017).
We are interested in estimating the constants for viscoelastic dampers (unknown parameters) β = [E ∞1 , ∆ 1 , Ω 1 , E ∞2 , ∆ 2 , Ω 2 ] that characterize a mechanical system with viscoelastic damping. In this approach, the parameters are modeled as a vector of random variables. In a Bayesian scenario, this inference is made by updating the a priori distribution p(β) with a posteriori distribution, namely, p(β|Z E , d), given that some specific data (experimental data) Z E were observed for certain parameters of design d. Note that the design parameters must be determined during the experiment execution process, and have no influence on previous knowledge about β, which suggests that a priori p(β) is independent of d. The posterior distribution is updated according to Bayes' theorem where the multidimensional integral p(Z E |d) = p(Z E |β, d)p(β)dβ is referred to as the evidence.
Using a Shannon information approach, one can define the information gain after updating the distribution of β, as being the Kullback-Leibler (KL) divergence from the a priori and a posteriori distributions, given by Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 226-01, 226-17, 2020 DOI: 10.21575/25254782rmetg2020vol5n21163 where integration is taken in the space of the Θ parameters. We are interested in quantifying the information gain before the data is actually collected. Therefore, we define the expected information gain where the integration in relation to Z E is taken in the space of the Y outputs. This can be interpreted as follows: The contribution of any possible output Z E , which is used as a data set, to update the a posteriori distribution is given as the KL divergence of the two distributions. From a mathematical point of view, ] is a function of Z E , therefore, the information gain is a random quantity and its expectation in relation to Z E , which is U (d), gives a measure of the average information to be obtained. This is a function only of the design parameters d that admit values in a design space D and, of course, suggests that determining provides, on average, the most informative observations, and thus is the optimal experimental design. From Monte Carlo estimators, one can define the estimatê where {β i } and {β j }, i = 1, . . . , N , j = 1, . . . , M , are distributed samples independently and identically (iid) of p(β) and {Z i E } are iid samples of p(Z i E |β i , d), where M and N are the lengths of the samples. Thus, our goal is to find the optimal positioning of the sensor so that the expected information gainÛ (d) is maximized.

INVERSE PROBLEM FORMULATION
The final solution of the inverse problem can be affected by uncertainties, which are inherent to experimental errors and/or to the computational model, for example. Such uncertainties can be taken into account from a probabilistic approach, where the parameters of the problem are modeled as random variables. This type of problem characterization is adopted Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 226-01, 226-17, 2020 DOI: 10.21575/25254782rmetg2020vol5n21163 in the so-called Bayesian inference, whose objective is to generate a posteriori probability distribution of the parameters of interest, based on a priori information on their values and on the basis of available experimental data.

Bayesian Inference
After determining the optimal positioning of the sensors, that is, after determining d * , one can then solve the inverse problem of parameter estimation considering d * . As of now, the probability distributions no longer depend on d * . Therefore, from the Bayesian point of view, the solution of the inverse problem, given the experimental observations Z E , is a function of a posteriori probability density of β, denoted by p(β|Z E ), which can be written, according to Bayes' theorem, as where p(Z E |β) is the likelihood function, p(β) is the a priori probability density of β and p(Z E ) is the marginal density (SCHWAAB; PINTO, 2007).
When a posteriori distribution can be obtained through analytical expression, the solution to the problem follows directly from Bayes' theorem. Otherwise, its determination becomes too complicated or impossible. In these cases, sampling techniques are widely used, where the simulation of samples of the distribution of interest is sought p(β|Z E ), so that inferences for this probability distribution can be made, such as obtaining statistical moments of mean and standard deviation.

Adaptive Monte Carlo Markov Chain Method
Specific algorithms are used to obtain Markov chains, where in this work the Metropolis-Hastings algorithm will be used, which makes use of an auxiliary probability density function q, to obtain samples of the posterior probability distribution p(β|Z E ). The algorithm can be specified by the following steps (KAIPIO; SOMERSALO, 2004).
The initial states, generated until the distribution equilibrium is reached, are called burn-in samples (KAIPIO; SOMERSALO, 2004), whose length will be denoted by N b . Such samples are discarded and the remainder of the states of the chain are considered in the inference of the statistical properties of the parameters of interest.
Although the Metropolis-Hastings algorithm, used to implement the MCMC, is robust and presents, in general, good results, it is known that its efficiency can be improved through a careful adjustment of the auxiliary probability density q, from which samples In this work, an adaptation of the Metropolis-Hastings Algorithm, proposed by Haario, Sakaman e Tamminen (2001) was considered, where the standard deviation of the auxiliary transition distribution is adjusted based on the covariance of the chain data. The algorithm is able to tune the auxiliary distribution as the chain explores the target distribution, along the construction of the chain.
In this adaptation, an auxiliary Gaussian distribution with mean in the current state (β j−1 ) was considered with the covariance C j = C j (β 0 , β 1 , . . . , β j−1 ). A crucial point for this adaptation is how the covariance depends on the history of the chain. For this, an initial covariance C 0 was considered strictly positive and after an initial non-adaptive period t 0 < N mcmc of iterations the covariance is calculated as follows where S Np is a parameter that depends only on the dimension of the problem N p and > 0 is a constant that can be chosen very small in relation to the covariance of the data, ensuring the non-singularity of the matrix, and I Np is the identity matrix with dimension N p . After the burn-in stage, the adaptation is no longer performed, and then the covariance remains constant and equal to the covariance obtained at the end of burn-in stage.

NUMERICAL RESULTS
In this work, as shown in Fig. 1, a simple mechanical system with two degrees of freedom was considered, in order to estimate the parameters that are related to viscoelastic damping, To obtain the experimental data used in the parameter estimation process, the impulsive response of the mechanical system, measured in terms of displacement, was considered. Markov chains with N mcmc = 20, 000 states were used, and a burn-in stage of N b = 12, 000 states was used to infer the statistical properties. The noise, added to the experimental data, is indirectly defined by the signal-to-noise ratio, given by where P s and P n correspond to the signal strength and the power of the noise, respectively, where in this work only a noise level of 30dB was used for the simulations. As an auxiliary distribution, it was considered a normal distribution with constant and equal standard deviation for all parameters σ k = 0.04, k = 1, . . . , N p , where N p = 6 is the number of estimated parameters. The cases considered in the simulations are described in Table 1. In this work, it was considered that the variable d, which describes the experimental design, that is, it describes the positioning of sensors and actuators, is discrete and is given by d = [d l,m ], where d l,m corresponds to the experimental design considering the actuator in the y l position and the sensor in the y m position, as they are two actuators and two sensors, the experimental design adopted is given by d = [d 1,1 , d 1,2 , d 2,1 , d 2,2 ]. Table 2 shows the expected information gainÛ (d) that was estimated for each case studied. Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 226-01, 226-17, 2020 DOI: 10.21575/25254782rmetg2020vol5n21163 Tables 3 and 4 show the results obtained from the inverse problem of estimation of viscoelastic parameters for the cases described in Table 1. The exact values (β exact ), the estimated averages (µ est,ij ), and the relative errors (E r,ij ) between the estimated value and the exact value are presented. The relative error is defined as  In Figure 2 are shown the Markov chains obtained in the process of estimating the viscoelastic parameters for the cases presented in Table 1 As can be seen in Table 2, the Bayesian experiment project identified the position d 2,2 of the experimental design d as being the one with the greatest gain of information, that is, placing the sensor and the actuator in the coordinate y 2 , the gain of information is greater in relation to the other positions. Tables 2 and 3 and Fig. 2 corroborate with this result, where it is observed that the values obtained in the parameter estimation process presented values closer to the exact values and, consequently, lower values of relative error, as well as converged Markov chains after approximately 10, 000 states. On the other hand, when positioning the actuator and the sensor at the y 1 coordinate, a smaller gain of information is observed in relation to the other positions, indicating that d 1,1 is the least indicated position to place the sensor and the actuator, thus leading to the estimation of values more distant from the exact values, higher values of relative error and non-converged Markov chains. The positions d 1,2 and d 2,1 presented similar values, beingÛ (d 1,2 ) >Û (d 2,1 ), thus showing that, for the problem addressed in this work, these two configurations present an intermediate gain of information. Figure 3 shows the posterior PDFs estimated considering the worst and the best experimental design, that is, Case 1 (d 1,1 ) and Case 4 (d 2,2 ), where it can be seen that, for the worst case, the estimated PDFs have low quality and multimodal behavior, showing the non-Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 226-01, 226-17, 2020 DOI: 10.21575/25254782rmetg2020vol5n21163 convergence of Markov chains, while, for the best case, the estimated PDFs showed behaviors close to a normal distribution, showing that the obtained chains are converged and that this design allows a better sampling of PDFs.

CONCLUSIONS
The present work had as its main objective the application of the Bayesian experimental project in the determination of an optimal experimental design/arrangement for the estimation of viscoelastic parameters of a mechanical system of two degrees of freedom. The formulations for the direct and inverse problems were presented. The solution of the inverse problem was obtained by the adaptive Monte Carlo Markov chain method proposed by Haario, Sakaman e Tamminen (2001). The experimental data used in the formulation of the inverse parameter estimation problem were corrupted with a noise of 30dB. It was observed that the Bayesian experiment design indicated as the optimal configuration, that is, the design that obtained the greatest gain of information, the one that considered the simultaneous positioning of the actuator and sensor in the y 2 coordinate. Furthermore, it indicated as the worst design the one that considered the simultaneous positioning of the actuator and sensor at the y 1 coordinate. Among the other designs, there was no significant difference in information gain. As suggestions for future work, it is proposed the application the Bayesian experiment design to more complex structures, i.e. systems with more degrees of freedom, such as beams and plates.