Time domain force localization and reconstruction based on hierarchical Bayesian method

Journal Pre-proof Time domain force localization and reconstruction based on hierarchical Bayesian method Wei Feng, Qiaofeng Li, Qiuhai Lu, Bo Wang, Chen Li PII:

S0022-460X(20)30053-5

DOI:

https://doi.org/10.1016/j.jsv.2020.115222

Reference:

YJSVI 115222

To appear in:

Journal of Sound and Vibration

Received Date: 13 August 2019 Revised Date:

22 January 2020

Accepted Date: 27 January 2020

Please cite this article as: W. Feng, Q. Li, Q. Lu, B. Wang, C. Li, Time domain force localization and reconstruction based on hierarchical Bayesian method, Journal of Sound and Vibration (2020), doi: https://doi.org/10.1016/j.jsv.2020.115222. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

Wei Feng: conceptualization, methodology, software, investigation, writing – original draft, visualization Qiaofeng Li: methodology, validation, resources, writing – review & editing, supervision, project administration Qiuhai Lu: conceptualization, methodology, resources, writing – review & editing, funding acquisition Bo Wang: investigation, resources Chen Li: investigation, resources

2

Time domain force localization and reconstruction based on hierarchical Bayesian method

3

Wei Fenga , Qiaofeng Lia,b,∗, Qiuhai Lua , Bo Wanga , Chen Lia

1

4 5

6

a Applied

Mechanics Laboratory, School of Aerospace Engineering, Tsinghua University, Beijing, 10084, China of Mechanical Engineering, Virginia Tech, Blacksburg, VA 24061, United States of America

b Department

Abstract In engineering practice, localization and reconstruction of external excitations is a tough problem in absence of prior knowledge. Traditional regularization approaches require assigning regularization parameters and shape parameters before implementing the reconstruction. However, a poor selection of these parameters generally leads to a poor reconstruction. In this paper a time domain hierarchical Bayesian method is proposed to locate and reconstruct the external excitations with automatic selection of parameters. The entire method is performed in two stages. First based on the property of posterior probability distribution function of shape parameters, a novel non-force criterion is proposed to determine the non-force locations quickly. Then based on the determined force locations, a reduced-dimension problem is constructed and the posterior distributions of parameters are sampled by a Metropolis-within-Gibbs sampler with nested blocking technique. Both numerical simulations and laboratory experiment of a cantilever beam under various load conditions are carried out to validate the proposed method.

8

Keywords: hierarchical Bayesian formulation, force localization and reconstruction, Markov chain Monte Carlo method, regularization, inverse problem

9

1. Introduction

7

10 11 12 13 14 15 16 17 18 19 20

Information of external excitations on structures is useful for predicting structure responses and establishing structural health monitoring formulations [1, 2]. For example, it is vital to quantify the external impacts induced by birds or other debris for structural health evaluation before the damage accumulates in aero-engine fan blades [1]. Wave loads can be identified in order to monitor the potential damages at remote or possibly subsea locations [2]. However, the excitations are in many conditions unknown or only partially known because direct measurements are implausible. In direct measurements, force transducers are generally considered to measure time histories of forces . The main drawbacks of the direct measurements are threefold. First, force transducers are sometimes difficult to install at target locations on engineering structures. Second, mounted force transducers may bring unexpected added stiffness or mass to the concerned structures. Last, when the force locations are unknown, the number and locations of force transducers can not be determined in advance. Indirect measurement is an alternative approach to solve these problems. ∗ Corresponding

author Email address: [email protected] (Qiaofeng Li)

Preprint submitted to Journal of Sound and Vibration

January 28, 2020

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

Generally, if information like mass, damping, and stiffness of the concerned structure and responses like displacements or accelerations induced by certain excitations can be obtained, the time histories of these excitations can be inversely reconstructed by some techniques. Unfortunately, the inverse problem is usually ill-posed which means an inevitable small amount of noise will cause a great variation of the reconstructed forces [3]. Thus a naive solution (e.g. least squares solution) will be unreliable. A widely investigated technique to solve this problem is regularization. It adds a penalty term to the original instinct least-squares problem to stabilize the solution. According to the types of penalty, different regularizations are investigated to handle different problems. The most well-known one is classic Tikhonov regularization [3–7], which has a regularization term of `2 −norm. Generally, when the force profile exhibits a certain continuity, the Tikhonov regularization has a great performance. The amount of penalty/regularization is controlled by the regularization parameter λ determined by selection criterions like L-curve and Generalized Cross Validation (GCV). Another attractive technique is the Lasso regularization [8–12], which considers a `1 −norm regularization term. This term introduces a sparse solution into this inverse problem, which means that only a few nonzero elements exist in the force history. However, in engineering practice, the external excitations exerted on a structure could be arbitrary. The sparsity or continuity assumption of the excitations is sometimes inapplicable in such conditions. A choice between `1 −norm and `2 −norm is demanded and is called `q −norm in some literatures. Li and Lu [13] proposed an automatic and adaptive method to determine the parameter q according to the force profile based on a Bayesian formulation. Aucejo [14] proposed a Generalized Iteratively Reweighted Least-Squares (GIRLS) algorithm to deal with lp −lq regularization problem and a great performance of this algorithm was shown in a simply supported plate validation. Liu et al. [15] proposed a novel technique for dynamic load identification based on shape function method of moving least square fitting. Feng et al. [16] proposed a Bayesian inference regularization in time domain to simultaneously identify bridge structural parameters and vehicle axle loads. Most above developed force reconstruction approaches assume that the locations of the excitations are known in advance, which is often not the case in engineering practice. A localization and reconstruction problem is generally more complicated than the reconstruction-only ones. The system matrix is not fully known and computation scale is large since the force locations are undetermined. Different methods have been presented in the literatures [17–24] to solve the force localization and reconstruction problems. Maes et al. [17] proposed a recursive algorithm to estimate the forces applied to a structure and the corresponding system states, with validation on a footbridge. Worden and Staszewski [18] proposed a Genetic Algorithm to train two artificial neural networks to locate and quantify the impact on a composite panel. Rezayat et al. [19] developed a new penalty function called mixed `1&2 −norm consisting of layered `1 and `2 to localize and reconstruct a force from seven potential locations on a cantilever beam. Maes et al. [20] proposed a joint input-state estimation algorithm identification several types of forces applied to the bridge deck. Li and Lu [22] proposed a two-step approach to localize and reconstruct the sources in time domain. The first step is to localize the force with a certain regularization parameter using a complex method. The second step is that repeat the localization process with different regularization parameters. However, most proposed methods need prior information on the shape parameters of forces. It is generally troublesome for researchers/engineers to make decisions on shape parameters before implementing the reconstruction process. A convenient and efficient method to deal with such problems is the Bayesian framework [25–32]. In the Bayesian framework, all concerned and unknown parameters are inferred as part 2

82

of the reconstruction problem. Li and Lu [26] adopted a hierarchical Bayesian method to reconstruct a single force in time domain. The shape parameter was regarded as a random variable to be inferred in the hierarchical structure. Markov chain Monte Carlo (MCMC) algorithm was proposed to draw samples from posterior distribution to make post-evaluation. Aucejo and De Smet [28] extended Bayesian formulation into a full Bayesian inference for several forces reconstruction in frequency domain. This method was validated on a steel parallelepiped box to reconstruct the amplitude spectrum of two kinds of forces. The credibility of the identified solution and posterior uncertainties of unknown parameters were discussed. In this paper, a novel hierarchical Bayesian formulation is proposed to localize and reconstruct multiple forces from multiple potential force locations in time domain. It considers all to-be-determined parameters as random variables to be inferred from the measurements of structure responses. Little prior information, including that on shape parameter, is required in this reconstruction process. So it enjoys automatic data-driven determination of both the force location and shape parameter. A MCMC algorithm, specifically a Metropolis-within-Gibbs sampler with nested blocking technique, is proposed to draw samples from the posterior distributions. To deal with the large dimension of the problem, the entire procedure is divided into two stages: first the force locations are determined based on a proposed novel non-force criterion, and then a full Bayesian inference is run on a reduced-dimension problem. This article is divided into five parts. The rest four parts are summarized as follows. The hierarchical Bayesian formulation is introduced in Section 2. The posterior distribution functions of to-be-determined parameters are presented. Section 3 introduces the Markov chain Monte Carlo sampling algorithm and the non-force criterion. Numerical and experimental validations are presented in Section 4 and 5. Section 6 gives the conclusion.

83

2. Hierarchical Bayesian formulation

84

2.1. Problem statement

62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81

85 86 87 88 89 90

In this section, we establish the hierarchical Bayesian formulation for localization and reconstruction of the applied forces. Specifically, the forces realized by random process are not considered in this paper. The structure is supposed to be linear elastic and time invariant. Several forces are applied on this structure at unknown locations. First we consider the condition that a single force is exerted on the structure at a ¯ If the known location. The behavior of this structure can be completely described by the system matrix H. ¯ ¯ responses of the structure W are caused by the excitation F, the linear structure is described as follows: ¯ =H ¯F ¯ W

91

(1)

where ¯ = W

h

T

w (1)

T

w (2)

···

T

w (nt − 1)

92

¯= F 93 94 95

h

f (1) f (2) · · ·

T

w (nt )

f (nt − 1) f (nt )

iT

iT

¯ ∈ Rns nt ×1 is the computed acceleration vector. ns is the number of accelerometers. nt is the length W ¯ ∈ Rnt ×1 is the force vector. H ¯ can be constructed by state-space model [26] or of force history vector. F obtained experimentally by hammer test [13, 33]. 3

96 97 98

Since the structure is linear elastic, superposition principle can be applied. The structural responses of multiple forces are equivalent to the sum of the responses caused by each excitation individually. When r various forces are exerted on the structure at the same time, the responses of the structure will be:  W=

r X i=1

99 100 101 102 103 104 105

106

¯i= W

r X

¯ iF ¯i = H

h

¯1 H

¯2 H

108 109 110 111

i    

    = HF  

(2)

The responses of structure in this situation can also be represented as the product of assembled system matrix and the assembled force vector. Note here r force vectors are corresponding vectors at r different potential force locations. Generally we consider r to be larger than the number of forces actually exerted on the structure. So only a fraction of F are non-zero vectors, corresponding to real force locations. The rest zero vectors correspond to non-force locations. Practically, we add a zero-mean Gaussian white noise vector N ∈ Rns nt ×1 to simulate the noise and modeling errors. X = W + N = HF + N (3) The Signal-to-Noise Ratio (SNR) is defined as SNR = 20log10

107

···

i=1

¯r H

¯1 F ¯2 F .. . ¯ Fr

kWk2 kNk2

The force reconstruction is complicated since it is generally an ill-posed inverse problem and many parameters need to be carefully estimated. The Bayesian framework is an effective way to solve this problem. It regards the to-be-determined parameters of the inverse problem as random variables to infer. The estimation and uncertainty of each parameter can be quantified from its posterior probability distribution. First, we need to introduce the Bayes’ rule from mathematical standpoint. p (F |X ) ∝ p (X |F ) p (F)

(4)

120

where p (F |X ) is the posterior probability distribution function of the unknown force history, representing the probability that the force history is F if the structural response is X; p (X |F ) is the likelihood function, representing the probability that the structural response is X when the force history is F; p (F) is the prior probability distribution function of the unknown force vector, representing the prior information before implementing the force reconstruction. Generally, the choice of the likelihood function and the prior probability distribution function strongly influences the quality of the force reconstruction [28]. If we get the exact information of the likelihood function p (X |F ) and prior probability distribution function p (F), a reasonable estimate of the unknown forces can be obtained from the posterior distribution.

121

2.2. Likelihood function

112 113 114 115 116 117 118 119

122 123

As mentioned previously, the likelihood function represents the probability of measured responses when the forces are known. The relationship between the force vector and responses is influenced by the noise N. 4

124

According to Eq. (3), the likelihood function can be derived as follows. p (X |F ) = p (HF + N |F ) = p (N |F ) = p (N)

125 126 127 128 129 130

131

The last equality holds because the noise N is assumed irrelevant to the force history F. It means that the likelihood function is equivalent to the prior probability distribution function of the noise. Practically, we assign the prior distribution of the noise as a zero-mean Gaussian distribution [26, 28] with precision parameter τn . This choice can limit the number of hyperparameters of the Bayesian model and obtain satisfactory reconstructions. The components of N are assumed to be independent and identically distributed (i.i.d). r ns nt   τn 2 p (N |τn ) = exp −τn kNk2 (6) π where k·k2 is the `2 − norm. Substitute Eq. (3) into Eq. (6) r p (X |F, τn ) =

132

133 134 135 136 137

138 139 140 141 142 143 144

(5)

τn π

ns nt

  2 exp −τn kX − HFk2

(7)

2.3. Prior probability distribution function The prior probability distribution function reflects the prior knowledge of researcher on the to-be¯ i exerted on the structure, its distribution is assumed as a determined excitations. For a single force F multivariate generalized zero mean Gaussian distribution [28, 29] with precision parameter τfi and shape parameter qi . And its components are also assumed to be independent and identically distributed (i.i.d) [26]. nt  

qi 
qi ¯ τfi nt /qi exp −τfi ¯ Fi q (8) p Fi |τfi , qi = i 2Γ (1/qi ) R +∞ where Γ (x) = 0 tx−1 e−t dt is the Gamma function. The reconstruction performance strongly relies on the value qi . In general, if qi = 2, a continuous-profile force, such as cosine force, is more likely the solution of force reconstruction; if qi = 1, a sparse-profile force, such as an impact, is more likely the solution. If various excitations are simultaneously exerted at r possible locations of the structure, the prior probability distribution function of assembled force vector F can be written as the product of the distribution ¯ i since the latter is considered as independent and uncorrelated to one another functions of single force F [28].
r ¯ i |τfi , qi (9) p (F) = Π p F i=1

145 146

Since different forces have different locations and distributions, we need to assign specific parameters τfi and qi for each force. The final expression of prior probability distribution function is r

p (F |τfi , qi , (i = 1, . . . , r) ) = Π

i=1

147



qi 2Γ (1/qi )

nt τfi

nt /qi

 exp −

r X

!

qi τfi ¯ Fi q i

(10)

i=1

By substituting Eq. (7) and Eq. (10) into Eq. (4), we can obtain the Maximum a Posterior (MAP) estimate

5

148

of F as follows. 2

FMAP = arg max p (F |X ) = arg max log (p (F |X )) = arg min kX − HFk2 + F

F

F

r

qi 1 X τfi ¯ Fi q i τn i=1

(11)

150

From Eq. (11), we observe that in Bayesian framework the ratio between τfi and τn is equivalent to regularization parameters in point estimators, e.g., denoted as λ in Tikhonov regularization method.

151

2.4. Hyperprior probability distribution function

149

152 153 154 155 156 157 158 159 160 161 162 163 164

Since the parameters qi , τn and τfi have a major influence on the reconstructed F, the standard Bayesian formulation is limited in some aspects. Manual selection of these parameters may lead to a poor reconstruction. An automatic data-driven method is preferred in this paper. Therefore, the parameters qi , τfi , and τn are treated as random variables and their posterior distributions are also to-be-inferred. Since qi , τfi , and τn act as parameters of the prior distribution of F and N, their prior distributions are called hyperprior distributions. As we have mentioned, q is a random variable and inferred based on data. The prior probability function of q has been discussed in some literatures [28, 34, 35]. In [28, 29], the possible value of q is bounded and practically lies in the interval [0,2] for force reconstruction problems. However, the advantage of quasi-norm regularization (q < 1) over `1 regularization in time domain force reconstruction problems requires further investigation since the latter already provides good sparsity-inducing effect. What’s more, the convexity of the problem is not guaranteed when q < 1, while q ∈ [1, 2] can produce a unique solution for point estimators like Eq. (11). Consequently the the prior distribution for qi in this paper is chosen as p (qi ) ∝ U [1, 2] , i = (1, 2, . . . , r)

165 166 167

where U [·, ·] represents the uniform distribution. Similarly, τn and τfi are also random variables to infer. In this way, a manual selection of equivalent regularization parameter λi = τfi /τn is avoided. The prior distributions of τn and τfi are Gamma distributions, p (τn |αn , βn ) =

168

p (τfi |αfi , βfi ) = 169 170 171 172 173 174 175 176

(12)

βn αn αn −1 τ exp (−βn τn ) Γ (αn ) n

βfi αfi αfi −1 τ exp (−βfi τfi ) , i = (1, 2, . . . , r) Γ (αfi ) fi

(13)

(14)

where αn , βn , αfi and βfi are hyperparameters of the distributions. The reason for this selection is that Gamma distribution is the conjugate prior for a generalized Gaussian distribution [36]. A conjugate prior is an algebraic convenience, giving a closed-form expression for the posterior distribution which can simplify the inference. This property makes it efficient to draw samples by MCMC algorithms. What’s more, this choice can guarantee the positivity of τn and τfi . However the hyperparameters αn , βn , αfi , and βfi reflect the unknown prior information on the precision parameters. These values should be set small for an automatic, data-driven determination of the parameters. In this paper, αn = βn = αfi = βfi = 1 × 10−16 .

6

177

178 179

2.5. Posterior distribution By repeating Bayes’ rule, we can obtain the posterior distribution of each to-be-inferred parameter. The force vector F ! r X

qi 2 p (F |X, τn , τfi , qi , (i = 1, . . . , r) ) ∝ p (X |F, τn ) p (F |τfi , qi ) ∝ exp −τn kX − HFk − τfi ¯ Fi 2

qi

i=1

(15) 180

The precision parameters τn p (τn |F, X, τfi , qi , (i = 1, . . . , r) ) ∝ p (X |F, τn ) p (τn |αn , βn )       n n α + s t −1 2 2 ∝ τn n 2 exp − βn + kX − HFk2 τn = Γ αn + ns2nt , βn + kX − HFk2

181 182

183

2

ns nt 2

and kX − HFk2 . The

n n  s t 2 , kX − HFk2 2

(17)

When αn and βn are very small, they can be ignored compared to values posterior probability distribution function of τn can be rewritten as p (τn |F, X, τfi , qi , (i = 1, . . . , r) ) ≈ Γ The precision parameters τfi

!   n

qi   αf + q t −1 ∝ p (F |τfi , qi ) p (τfi ) ∝ τfi i i exp − βfi + ¯ Fi q τfi p τfi X, F, τn , τfj , qi i i=1:r j=1:r,j6 = i 

qi  nt = Γ αf + , βf + ¯ Fi i

184 185

186

187

qi

i

When αfi and βfi are very small, they can be ignored compared to values probability distribution function of τfi can be rewritten as

190 191

nt qi

qi and ¯ F q . The posterior i

!  

qi nt

¯

p τfi F, X, τn , τfj , qi ≈Γ , Fi qi qi j=1:r,j6=i i=1:r

(19)

!  nt 

qi  qi τfi nt /qi exp −τfi ¯ p qi F, X, τn , τfi , qj Fi qi [[[1,2] (qi ) ∝ 2Γ (1/qi ) i=1:r j=1:r,j6=i

(20)

The shape parameters qi

where [[[1,2] (qi ) is the truncation function [[[1,2] (qi ) =

189

(18)

qi

(

188

(16)

1, qi ∈ [1, 2] 0, otherwise

3. Markov chain Monte Carlo algorithm In general, the posterior probability distribution functions (PDFs) determine the shapes of distributions of to-be-determined parameters. However, the analytical solutions of the posterior PDFs of forces are intractable due to the high dimension of the parameter space [28]. An efficient method to solve this problem 7

209

is using MCMC algorithm to draw samples from their posterior distributions. When enough samples are obtained, the estimated values and their uncertainties of to-be-determined parameters can be evaluated by the samples. This section will introduce the details of this algorithm, initialization of parameters, convergence criterion, non-force criterion, and nested blocking technique. This paper adopts Metropolis-within-Gibbs sampler to draw samples from the target posterior probability distributions. The main idea of the Metropolis-within-Gibbs sampler [37] is the same as the Gibbs sampler. A random-walking Metropolis (RWM) step is applied to draw samples from posterior distributions of F and q. A new candidate is generated by adding a Gaussian noise on the last sample. It is immediately accepted when the posterior probability of the new candidate is higher than that of the last sample. Otherwise, the new candidate is accepted with a probability. If the step size of the Gaussian noise is inappropriately small, a lot of samples will gather in a narrow region leading to poor depiction of the entire distribution. So the RWM step size is automatically adjusted by the instantaneous acceptance rate of samples of force during computation. The acceptance rate is calculated each time 1000 samples are drawn. If the acceptance rate is larger or smaller than the threshold values, the stepsize will be increased or decreased by multiplying a constant value (0.8 and 1.2). In this paper, the acceptance rate intervals for F and qi are [0.25,0.35] and [0.4,0.5] respectively [37]. Technically, this algorithm is divided into five main steps. The first four steps focus on drawing samples from their posterior distributions. Please refer to Algorithm. 1 for details.

210

3.1. Initialization

192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208

211 212 213 214 215 216

In general, all samples with different initial values should achieve the same stationary distribution. However, a poor initial value may lead to a long burn-in period because the chains are hard to mix well. Due to the limitation of desktop computer capability, a failure of force reconstruction may occur when too many invalid samples are generated. A good estimated initial value will accelerate the sampling process. In this paper, the initialization is performed with a classical Tikhonov regularization reconstruction. Here, all to-be-determined forces are treated as one assembled force to infer.
−1 T 2 2 F(0) = arg min kHF − Xk2 + λ kFk2 = HT H + λI H X

(21)

F

217 218

(0)

λ is the regularization parameter, which can be determined by L-curve method. The initial values for τn (0) and τfi are ns nt τn(0) = (22)

2 2 HF(0) − X 2

219

(0)

τfi = λ ∗ τn(0)

(23)

(0)

221

Eq. (22) can be derived from Eq. (17). τn is assumed as the mean of the posterior probability distribution when F is F(0) . The initial value for qi is randomly chosen from the interval [1,2].

222

3.2. Convergence criterion

220

223 224 225

The chains with different initial values are expected to eventually converge to the same stationary distribution, which is the target distribution. It is important to introduce a criterion to determine whether the chains have actually converged and properly mixed. Several techniques have been developed in scientific 8

244

literatures [38–40]. Practically, the researches on convergence diagnostics have been divided into two categories. The first one considers one long chain while the other one uses several parallel chains to determine the convergence. In this paper, Potential Scale Reduction Factor (PSRF) [40] is used as the convergence criterion. τn , τfi , qi , and F are monitored in this process. Since τn , τfi , and qi are one-dimensional variables, their PSRFs are easy to calculate. However, the PSRF of F is hard to calculate since it is a high-dimensional variable. We consider kHF − Xk2 as an alternative one to simplify the problem. If the chains of kHF − Xk2 are converged, the chains of F are also considered converged. The process of the convergence diagnostics consists of five main steps. 1. Generate 3 parallel chains with different initial values. 2. Discard the first-half samples of each chain. 3. Calculate the within-chain and between-chain variances of the second-half samples of each parameter. 4. Calculate the estimated variances which are the weighted sum of the within-chain and between-chain variances. 5. Calculate the Potential Scale Reduction Factor. In this paper, 3 parallel Markov ichains are generated with different initial values. The initial value of h (0) (0) τn for each chain are 0.1 1 10 ∗ τn respectively, τn is the result of Eq. (22). The qi is selected from [1,2] randomly. For simplicity, the initial values of τfi and F are same for 3 chains. Burn-in period is another important issue when implementing the MCMC algorithm. We discard the first-half samples of each chain which means that the first half samples are regarded as burn-in phase here.

245

3.3. Non-force criterion

226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243

246 247 248 249 250 251 252 253 254 255 256 257 258 259

In engineering practice, the locations of forces are generally unknown a priori, and are sparse in space. If the force at a potential location is zero, the corresponding τfi and qi could be any value theoretically. In such cases, undoubtedly the Markov chains will not converge easily, the algorithm will be time-consuming, and many invalid samples may be produced. It is important to propose a criterion to determine whether a real force is acting at a particular location before fully diving into force reconstruction stage. It will accelerate the MCMC algorithm since the number of forces and dimension of problem are reduced. The statistics of the shape parameter qi at each location are promising indicators for force/non-force determination. If the force at location i is zero, the shape parameter qi could be any value in the interval [1,2]. The samples of qi will fill up the interval, which looks like a uniform distribution. So a criterion can be set based on this observation. In statistics, a quantile - quantile (Q - Q) plot is used to compare two probability distributions by plotting their quantiles against each other. A point (x, y) on the plot corresponds to one of the quantiles of the second distribution (y−coordinate) plotted against the same quantile of the first distribution (x−coordinate). If the two distributions being compared are similar, the points in the Q - Q plot will approximately lie on the line y = x. An indicator Q is constructed to determine whether there is a force according to the Q - Q plot. Q= Q=

260 261

1 N 1 N

N P i=1 N P i=1

NqS (i) NqU (i) ,

when NqS (i) ≤ NqU (i)

NqU (i) NqS (i) ,

when NqS (i) > NqU (i)

(24)

NqS (i) is the ith value that partitions a finite set of values into N subsets of equal size in samples drawn from qi . NqU (i) is the ith value that partitions a finite set of values into N subsets of equal size in samples 9

262 263 264 265

of uniform distribution. The quotient of these two quantiles is set smaller than 1. The indicator Q indicates the similarity between to-be-determined samples and uniform distribution. If Q = 1, these samples follow a uniform distribution. In this paper, Q > 0.9 indicates a non-force location. The mean and variance of uniform distribution can serve as supplementary indicators. E (qi ) =

266

R2

qi dqi 1 2−1

(25)

= 1.5


2 V ar (qi ) = E qi 2 − [E (qi )] =

1 12

(26)

= 0.0833

275

If the force at location i is zero, the mean and variance of shape parameter qi will be close to 1.5 and 0.833 respectively. On the contrary, a specific force will have a unique mean value and the variance will be much smaller than the uniform distribution. These two indicators are used to enhance the reliability of our non-force criterion. In this paper, we draw 50000 samples from posterior distributions of qi to calculate the three indicators mentioned above in advance. If the Q of qi is smaller than the threshold Q = 0.9, there is a force at location i. Otherwise, there is no excitation. Then the non-force locations are eliminated from the problem. The system matrix H is rebuilt based on only the identified force locations. This force localization procedure is performed before force reconstruction stage.

276

3.4. Nested blocking technique

267 268 269 270 271 272 273 274

277 278 279 280 281 282 283 284 285 286 287 288 289

Markov chain Monte Carlo methods are often deemed too computationally intensive to be practical for big data applications. As mentioned above, the step size of MCMC algorithm strongly influences the computation time. A relatively large RWM step size is required to reduce the autocorrelation of the samples produced by Metropolis-Hastings (MH) algorithm. However this may increase the rejection of the proposal jump. In fact, either too large or small RWM step size will lead to a slow mixing Markov chain, in which a large number of samples are needed for a reasonable estimate of the target distribution. For a multidimensional variable, finding the suitable RWM step size can be difficult, as the RWM step size must be "just right" for all dimensions to avoid excessively slow mixing. Various approaches to scale up the MH algorithm in a Bayesian inference context have been proposed in literatures [41–43], but are not suitable for the update of F in our specific hierarchical structure. In this paper, we propose a nested blocking technique to solve this problem. The main idea of this technique is to divide the high-dimensional variable into several blocks and choose a new sample for each low-dimensional block by turn. With this technique, larger random walking step size can be considered without damaging the acceptance rate. The details are listed as follows. 1. Divide the vector into m blocks as in Fig. 1. F = éë f1 (1)

f1 ( 2 ) ⋯

f1 ( nt )

f 2 (1)

f2 ( 2) ⋯

f 2 ( nt ) ⋯

f r (1)

fr ( 2) ⋯

f r ( nt ) ùû

r ´ nt F = éë f (1)

f ( 2) ⋯

Block 1

f ( n1 )

f ( n1 + 1)

f ( n1 + 2 ) ⋯

Block 2

f ( n2 ) ⋯⋯⋯

⋯

f ( nm -1 + 1)

f ( nm -1 + 2 ) ⋯

f ( nm ) ùû

Block m

Figure 1: The illustration of nested blocking technique. The force vector F is divided into m blocks. 290

10

291 292 293 294 295 296 297 298 299 300 301

2. Add a Gaussian noise on block 1 and determine whether to accept the new sample. If accepted, the variable can be updated by new sample. Otherwise, keep the original one. 3. Add a Gaussian noise on block 2 and determine whether to accept the new sample. 4. Repeat the above steps until all blocks are updated. The Metropolis-within-Gibbs sampler with nested blocking technique is summarized in Algorithm. 1. It involves five main steps: draw samples from F, τn , τf , and q by turn, and check convergence of the Markov chains. The whole procedure of the proposed force localization and reconstruction method can be summarized in 3 steps. 1. Consider r possible forces in different locations and use Algorithm. 1 to draw samples of qi . 2. Only r0 real force locations are detected in this localization process according to the non-force criterion. 3. Reconstruct the r0 real forces with Algorithm. 1. Algorithm 1: Metropolis-within-Gibbs sampler Input: system matrix H, measured responses X, initial vectors h i T T T T (0) (0) (0) ¯ (0) ¯ (0) ¯ (0) F(0) = F , τn , τfi , qi , i ∈ (1, 2, . . . , r) F · · · F r 1 2 Output: Samples of F, τn , τfi , qi , i ∈ (1, 2, . . . , r) while do 1. draw F(k) ∝ p (F |X, τn , τfi , qi ) while j=1:m do
draw gj ∝ N 0, εfj I , FI = F(k−1) + g0 j i .  h  (k−1) (k−1) (k−1) (k−1) (k−1) (k−1) , τfi , qi (Eq. (15)) , τfi , qi p F(k−1) X, τn t = min 1, p FI X, τn if t1 < t then F(k) = FI else F(k) = F(k−1) end

302

end   (k−1) (k−1) ∝ p τn F(k) , X, τfi , qi (Eq. (16))   (k) (k) (k−1) (Eq. (18)) 3. For each location i, draw τfi ∝ p τf F(k) , X, τn , qi   (k) (k) 4. For each location i, draw qi (k) ∝ p qi F(k) , X, τn , τfi (Eq. (20)) 5. Monitor the convergence of Markov chains. end Note:   T x ∈ Rs×1 ∝ N (µ, Σ) = √ 1 s exp − 21 (x − µ) Σ−1 (x − µ) (k)

2. draw τn

(2π) |Σ|

I∈R 0

rnt m

×

rnt m

. "

rnt ×1

gj ∈R

0

0

is the corresponding vector of gj . g j =

|

0 ··· {z

block 1,··· ,j−1

For conciseness, i = 1, . . . , r are eliminated in above expressions.

11

0 }

gj |{z}

block j

0 |

0 ··· {z

#

0

block j+1,··· ,m

}

303

304 305 306 307 308 309 310 311 312 313 314 315 316 317

4. Numerical validation In this section, three numerical cases are carried out to validate the proposed method. A cantilever beam is simulated with 8 accelerometers for response measurements and 7 potential force locations. The structure is shown in Fig. 2. The beam has a length of 1 m and a rectangular cross section with a width of  3 cm and a height of 1 cm. The Young’s modulus is 206 GPa. And the density is 7800 kg m3 . The beam is discretized with 50 Euler-Bernoulli beam elements. We only consider the bending motion in the vertical plane. Rayleigh damping is used as the damping model. And the damping ratio is ζ=0.02. That means C = αM + βK, where α = 2ζω1 ω2 /(ω1 + ω2 ) and β = 2ζ/(ω1 + ω2 ). ω1 and ω2 are the angular frequencies of the first two modes. The sampling frequency is 2048 Hz. The SNR (Signal-to-Noise Ratio) is 25dB. We set nt = 64 in the following. The whole force history is divided into 8 blocks evenly in non-force detection process. In force reconstruction process (after non-force detection), the block number can be reduced to 4 for saving computation time. In all cases, two forces are exerted on two different locations simultaneously. The first case considers 2 sparse triangle forces and the second is 2 continuous sine forces. The last one is a sparse triangle force and a continuous sine force. In all cases of numerical validation, 7 potential forces with 8 accelerometers are considered for the force reconstruction. According to [44, 45], this problem is instantaneously invertible.

F2

a1

a2

F5 p1

a3

p2

a4

p3

a5

Unknown Force Locations

p4

p5

a6

p6

a7

p7

a8

Accelerometers

Z

3cm 1m

Z

Y

Z 1cm

X

X

Y

Figure 2: The schematic of the simulated cantilever beam. A total of 7 unknown force locations are considered. Real forces are applied at point 2 and 5. 318

319 320 321 322 323 324 325 326 327 328 329 330 331

4.1. Case 1: two sparse triangle forces Two sparse triangle forces are first applied on location 2 and 5 respectively. The non-force detection should be carried out at the very beginning of force reconstruction. 50000 samples of each variable are drawn from its posterior distribution. The chains of qi are shown in Fig. 3. The second-half samples are concatenated to calculate the values of indicators like Q, M ean (q), and V ar (q) listed in Table. 1. The Q Q plot is shown in Fig. 4. As we discussed in Section 3.3, the values of Q on non-force locations are larger than 0.9 while for locations 2 and 5 they are 0.6744 and 0.6936 respectively. Distinguishable difference is observed between the Q-Q plots of real force and non-force locations. The figures of non-force locations show that the blue lines are close to the black lines. However, for locations 2 and 5 the distributions of samples are quite different from uniform distribution. Other indicators are the means and variances of qi . Fig. 3 shows that qi can be any possible value within the interval [1,2] when no force is exerted at a location. The variances and means will tend to 0.0833 and 1.5 respectively. 12

2

2

2

2

1.8

1.8

1.8

1.8

1.6

1.6

1.6

1.6

1.4

1.4

1.4

1.4

1.2

1.2

1.2

1.2

1

1

1 1

2

3

4

5

1

2

3

4

1

5

10 4

2

3

2

2

1.8

1.8

1.8

1.6

1.6

1.6

1.4

1.4

1.4

1.2

1.2

1.2

1

1 2

3

4

1

5

1

5

1

2

3

10 4

2

1

4

10 4

2

3

4

5 10 4

1

5

1

2

3

4

5

10 4

10 4

4

10 4

Figure 3: The Markov chains of qi at all locations.

Table 1: Values of indicators for non-force detection of two sparse triangle force case.

P1 0.9800 1.5301 0.0839

P3 0.9895 1.5074 0.0903

P4 0.9567 1.5640 0.0750

2

1.8 1.6 1.4 1.2

1.8 1.6 1.4 1.2 1

1 1.5

1.6 1.4 1.2

1.2

1.4

1.6

1.8

2

1.2

1.4

1.6

Quantiles of Input Sample

1.6 1.4 1.2 1 1

1.2

1.4

1.6

1.8

Quantiles of Uniform Distribution

2

1.8

1.8 1.6 1.4 1.2

1

2

1.5

2

Quantiles of Uniform Distribution

Quantiles of Uniform Distribution

2

2

1.8

P7 0.9920 1.5112 0.0816

1 1

Quantiles of Uniform Distribution

2

Quantiles of Input Sample

1.8

1 1

2

Quantiles of Uniform Distribution

P6 0.9890 1.5003 0.0913

2

Quantiles of Input Sample

1

P5 0.6936 1.0049 0.00002

2

Quantiles of Input Sample

Quantiles of Input Sample

2

Quantiles of Input Sample

P2 0.6744 1.0061 0.00004

Quantiles of Input Sample

Loc Q M ean (q) V ar (q)

1.8 1.6 1.4 1.2

1.8 1.6 1.4 1.2 1

1 1

1.2

1.4

1.6

1.8

Quantiles of Uniform Distribution

2

1

1.2

1.4

1.6

1.8

2

Quantiles of Uniform Distribution

Figure 4: The Q - Q plot of qi at all locations for the two sparse triangle force case. A blue line is the quantile-quantile line of this data. A dashed reference black line is shown by y = x.

13

332 333 334 335 336 337 338 339 340 341

However, the samples of locations 2 and 5 are gathered in a stable area and have specific values close to 1 for both locations. This is reasonable since this value of qi is a sparsity-inducing norm value and suits the profile of the force well. According to the non-force criterion, locations 2 and 5 have real excitations. Till now, the computation scale for afterwards force reconstruction can be reduced into two force reconstruction. Two forces are reconstructed with the proposed MCMC sampler. The Markov chains of 4 kinds of monitored variables are shown Fig. 5. The parallel Markov chains are considered well mixed. 30000 samples for each parallel chain are generated when 4 PSRFs are all below the threshold value 1.1. The second-half samples are concatenated to form a single chain of 45000 samples to produce the post-evaluation. The histogram of each variable is shown in Fig. 6. The mean value, median value, and 95% credible interval of monitored parameters are listed in Table. 2. Table 2: Mean, median and 95% CI of kHF − Xk2 , τn , τf2 , τf5 , q2 , and q5 for the two sparse triangle force case.

M ean M edian 95%CI

q2 1.0125 1.0086 [1.0003,1.0466]

q5 1.0117 1.0080 [1.0003,1.0433]

τn 4.1777 4.1707 [3.6157,4.7807]

τf2 0.6774 0.6738 [0.5089,0.8665]

τf5 0.6983 0.6945 [0.5288,0.8888]

kHF − Xk2 7.8294 7.8230 [7.5949,8.1015]

Figure 5: The Markov chains of kHF − Xk2 , τn , τf2 , τf5 , q2 , and q5 for the two sparse triangle force case. 342 343 344 345 346

The results show that there are two sparse forces exerted on the cantilever beam. The shape parameter and precision parameters of each force can be estimated based on the samples. The posterior mean, `1 regularization reconstruction, `2 regularization reconstruction, and 95% credible interval of F are compared in Fig. 7. In practice, the `1 regularization holds a good estimate of the true forces since they are sparse. However the posterior mean is a better one compared with `1 and `2 regularization. The `1 and `2 regu14

Figure 6: The histograms of kHF − Xk2 , τn , τf2 , τf5 , q2 , and q5 for the two sparse triangle force case.

Table 3: The number of samples and relative errors in different blocks conditions.

Blocks Samples RE of f1 (%) RE of f2 (%)

1 130000 3.56 6.00

2 60000 3.40 5.73

4 30000 4.33 3.69

8 15000 3.52 4.88

16 7000 3.54 5.22

32 3000 3.56 5.59

64 1600 2.81 5.66

353

larization results have oscillatory components where the true force stays zero. In addition, the narrow 95% credible interval shows high confidence of the estimated force history. The statistical analysis of the number of blocks is also performed for the force reconstruction process. 7 conditions (1, 2, 4, 8, 16, 32, 64 blocks) are investigated for this force reconstruction. The number of samples and the reconstruction errors in each condition are listed in Table. 3. As shown in the table, over a wide range of block numbers, the method can converge and produce reasonable estimates. In practice, the number of blocks should be determined based on the researcher’s experience.

354

4.2. Case 2: two continuous sine forces

347 348 349 350 351 352

355 356 357 358 359

In this case, two continuous sine forces are exerted on the locations 2 and 5. The frequency of the sine forces is 94Hz. The non-force detection is operated in the same way as the first case. The chains of qi are shown in Fig. 8. The second-half samples are concatenated to calculate the values of indicators like Q, M ean (q), and V ar (q) listed in Table. 4. The Q - Q plots are shown in Fig. 9. According to the proposed criterion, only locations 2 and 5 have non-zero excitations.

15

25

25

20

20

15

15

10

10

5

5

0

0

-5 0

0.005

0.01

0.015

0.02

0.025

95% CI true force L1 regularization L2 regularization mean force

-5

0.03

0

0.005

0.01

0.015

0.02

0.025

25

25

25

25

25

0

0

0

0

0

0

0

0.03

0.03

0

0.03

0

0.03

0

0.03

0.03

Figure 7: The estimated force histories of `1 , `2 regularization and proposed method with 95% credible interval.

2

2 1.8

1.8

1.6

2

2

1.8

1.8

1.6

1.6

1.4

1.4

1.2

1.2

1.6 1.4 1.4

1.2

1.2

1 1

2

3

4

1

5

2

3

4

1

1

5

1

2

3

4

10 4

10 4

2

2

2

1.8

1.8

1.8

1.6

1.6

1.6

1.4

1.4

1.4

1.2

1.2

1.2

1

1 1

2

3

4

5

1

1

5

2

3

2

3

4

5

1

10 4

10 4

Figure 8: The Markov chains of qi at all locations.

16

1

4

5 10 4

10 4

2

3

4

5 10 4

2

1.6 1.4 1.2

1.8 1.6 1.4 1.2 1

1 1.5

1

2

1.5

1.4 1.2

2

1

1.5

Quantiles of Input Sample

1.8 1.6 1.4 1.2

1

1.5

2

Quantiles of Uniform Distribution

1.6 1.4 1.2

1

1.5

2

Quantiles of Uniform Distribution

2

1.8 1.6 1.4 1.2 1

1

1.8

1

2

Quantiles of Uniform Distribution

2

2

Quantiles of Input Sample

1.6

1

Quantiles of Uniform Distribution

Quantiles of Uniform Distribution

1.8

Quantiles of Input Sample

1

2

Quantiles of Input Sample

1.8

2

Quantiles of Input Sample

Quantiles of Input Sample

Quantiles of Input Sample

2

1.8 1.6 1.4 1.2 1

1

1.5

2

Quantiles of Uniform Distribution

1

1.5

2

Quantiles of Uniform Distribution

Figure 9: The Q - Q plot of qi at all locations for the two since force case. A blue line is the quantile-quantile line of this data. A dashed reference black line is shown by y = x.

Table 4: Values of indicators for non-force detection of two sine force case.

Loc Q M ean (q) V ar (q)

P1 0.9595 1.5597 0.0763

P2 0.7746 1.9321 0.0041

P3 0.9609 1.5567 0.0740

P4 0.9477 1.5797 0.0783

P5 0.7741 1.9335 0.0041

P6 0.9666 1.5486 0.0769

P7 0.9688 1.5474 0.0835

Figure 10: The Markov chains of kHF − Xk2 , τn , τf2 , τf5 , q2 and q5 for the two continuous sine force case.

17

360 361 362 363 364 365 366 367 368 369 370

MCMC algorithm is then applied to reconstruct these two forces. The Markov chains of 4 kinds of monitored variables are shown in Fig. 10. The parallel Markov chains are considered well mixed. 30000 samples for each parallel chain are generated when 4 PSRFs are all below the threshold value 1.1. The secondhalf samples are concatenated to form a single chain of 45000 samples to produce the post-evaluation. The histograms of variables are shown in Fig. 11. The mean value, median value, and 95% credible intervals of monitored parameters are listed in Table. 5. The results show that there are two continuous forces acting on the cantilever beam. The posterior mean, `1 regularization result, `2 regularization result, and 95% credible interval of F are shown in Fig. 12. The 95% credible intervals of F are different between the two locations. The uncertainty of F2 which is close to the fixed end is larger than that of F5 which is close to free end. Compared with `1 and `2 regularization results, the posterior mean has a better performance on the reconstruction of F. Table 5: Mean, median and 95% CI of kHF − Xk2 , τn , τf2 , τf5 , q2 and q5 for continuous sine force case.

M ean M edian 95%CI

371 372 373 374 375 376 377

q2 1.8722 1.9070 [1.5697,1.9963]

q5 1.8837 1.9145 [1.5947,1.9966]

τn 0.4188 0.4181 [0.3624,0.4788]

τf2 0.0042 0.0034 [0.0020,0.0108]

τf5 0.0039 0.0033 [0.0020,0.0098]

kHF − Xk2 24.7362 24.7112 [23.9519,25.6348]

In order to demonstrate the advantages of proposed non-force criterion and nested blocking technique, computation times of proposed method and full dimensional method (without non-force criterion and nested blocking technique) are compared in Table. 6. As we discussed in Section 3.3, the τfi and qi could be any value if force at corresponding location is zero. The Markov chains of these variables are hard to meet the convergence criterion. We reconstruct the forces with the same dataset by sampling 900000 samples using full dimensional method. However, some of PSRFs of τfi and qi can not be less than 1.1. On the contrary, the proposed method can deal with this problem efficiently. Table 6: Computation times of proposed method and full dimensional method. Note that the proposed method is converged with 80000 samples (50000 for localization, 30000 for reconstruction); the full dimensional is still not converged after 900000 samples.

method time/s

378

379 380 381 382

proposed 296.11

full dimensional 2509.52

4.3. Case 3: a sparse triangle force and a continuous sine force In this case, two forces with different types are applied on this cantilever beam on location 2 and 5 respectively. The forces are a sparse triangle force and a continuous sine force which are the same as the case 1 and case 2 respectively. After the process of non-force detection, the information of qi can be obtained as follows.

18

Figure 11: The histograms of kHF − Xk2 , τn , τf2 , τf5 , q2 , and q5 for the two continuous sine force case.

25

25

10

10

0

0

-10

-10

-25

-25 0

0.005

0.01

0.015

0.02

0.025

0

0.03

0.005

0.01

0.015

25

25

25

25

25

0

0

0

0

0

-25

-25

-25

0.02

0.025

0.03

95% CI true force L1 regularization L2 regularization

0

0.03

0

0.03

mean force

0

0.03

-25

-25 0

0.03

0

0.03

Figure 12: The estimated force histories of `1 , `2 regularization and proposed method with 95% credible interval.

Table 7: Values of indicators for non-force detection of a sparse triangle force and a sine force case.

Loc Q M ean (q) V ar (q)

P1 0.9761 1.5358 0.0829

P2 0.6948 1.0067 0.00005

P3 0.9567 1.5648 0.0779

19

P4 0.9688 1.5479 0.0847

P5 0.7747 1.9317 0.0043

P6 0.9453 1.5808 0.0693

P7 0.9650 1.5515 0.0776

2

2

2

2

1.8

1.8

1.8

1.8

1.6

1.6

1.6

1.6

1.4

1.4

1.4

1.4

1.2

1.2

1.2

1.2

1

1

1

1

2

3

4

1

5 10

2

3

4

5

4

10

1

2

3

4

4

10

2

2

2

1.8

1.8

1.8

1.6

1.6

1.6

1.4

1.4

1.4

1.2

1.2

1.2

1 2

3

4

5 10

1

2

3

4

5 10 4

1

1 1

1

5 4

1

2

3

4

4

1

5 10

2

3

4

5 10 4

4

Figure 13: The Markov chains of qi at all locations.

2

1.6 1.4 1.2

1.8 1.6 1.4 1.2

1.8 1.6 1.4 1.2

1

1 1.5

1.5

1

Quantiles of Uniform Distribution

Quantiles of Uniform Distribution

1.2

1.4

1.6

Quantiles of Input Sample

1.8 1.6 1.4 1.2

1

1.5

2

Quantiles of Uniform Distribution

1.4 1.2

1.8

1

2

1.5

2

Quantiles of Uniform Distribution

2

1.8 1.6 1.4 1.2 1

1

1.6

Quantiles of Uniform Distribution

2

2

Quantiles of Input Sample

2

1.8

1

1 1

2

Quantiles of Input Sample

1

Quantiles of Input Sample

1.8

2

2

Quantiles of Input Sample

Quantiles of Input Sample

Quantiles of Input Sample

2

1

1.2

1.4

1.6

1.8

Quantiles of Uniform Distribution

2

1.8 1.6 1.4 1.2 1 1

1.2

1.4

1.6

1.8

2

Quantiles of Uniform Distribution

Figure 14: The Q - Q plot of qi at all locations. A blue line is the quantile-quantile line of this data. A dashed reference black line is shown by y = x.

20

Figure 15: The Markov chains of kHF − Xk2 , τn , τf2 , τf5 , q2 , and q5 for the sparse triangle force and continuous sine force case.

Figure 16: The histograms of kHF − Xk2 , τn , τf2 , τf5 , q2 , and q5 for the sparse triangle force and continuous sine force case. Table 8: Mean, median and 95% CI of kHF − Xk2 , τn , τf2 , τf5 , q2 , and q5 for the sparse triangle force and continuous sine force.

M ean M edian 95%CI

q2 1.0143 1.0098 [1.0004,1.0536]

q5 1.8835 1.9167 [1.5847,1.9967]

τn 0.8450 0.8431 [0.7290,0.9687]

21

τf1 0.6488 0.6451 [0.4840,0.8348]

τf2 0.0039 0.0033 [0.0020,0.0100]

kHF − Xk2 17.4089 17.3915 [16.8360,18.0681]

383 384 385 386 387 388 389

The variances and means of qi s at non-force locations are close to 0.0833 and 1.5 respectively according to Table. 7. The more samples are drawn, the more accurate the criterion will be. According to the results, the locations of real forces are 2 and 5. The Markov chains of 4 kinds of monitored variables are shown in Fig. 15. The parallel Markov chains are considered well mixed. 30000 samples for each parallel chain are generated when 4 PSRFs are all below the threshold value 1.1. The second-half samples are concatenated to form a single chain of 45000 samples to produce the post-evaluation. The histograms of variables are shown in Fig. 16. The mean value, median value, and 95% credible interval of monitored parameters are listed in Table. 8. 25

25 95% CI true force L1 regularization L2 regularization mean force

20 15

10

10

0

5

-10

0 -5 0

0.005

0.01

0.015

0.02

0.025

-25

0.03

0

0.005

0.01

0.015

0.02

0.025

25

25

25

25

25

0

0

0

0

0

-25

-25 0

0.03

-25 0

0.03

0

0.03

0.03

-25

-25 0

0.03

0

0.03

Figure 17: The estimated force histories of `1 , `2 regularization and proposed method with 95 % credible interval. 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404

The posterior mean, `1 regularization result, `2 regularization result and 95% credible interval of F are shown in Fig. 17. The results show a great performance on the force reconstruction and the true forces are included in their 95% credible intervals. The shape parameters of these two forces are q2 = 1.0143 and q5 = 1.8835 respectively. Different types of forces are determined different shape parameters by the proposed algorithm. Compared with the `1 regularization and `2 regularization results, the proposed method enjoys the merit of automatic data-driven determination, and thus stronger robustness across different cases. Especially in situations where prior knowledge on the force nature is unavailable, the proposed method can adaptively choose the most suitable q, and achieve satisfactory results. In these three cases, 4 blocks are adopted in force reconstruction process (after non-force detection). And the iteration diagram of step size for these three cases is shown in Fig. 18. The step sizes have all become stabilized within the burn-in phase of MCMC. The reconstruction errors of the estimation of forces at different position in these three cases are summarized in Table. 9. It is observed that the proposed method does not always achieve the best results. But it has the benefit of adaptivity and robustness, i.e., the results are satisfactory for different types of cases.

22

0.08

0.1

0.1

0.07

0.09

0.09

0.08

0.08

0.07

0.07

0.06

0.06

0.06 0.05 0.04 0.03

0.05

0.02 0

1

2

0

3 10

1

2

0.05

3

4

10

0

1

2

4

3 104

Figure 18: The iteration of step sizes for three cases with 4 blocks.

Table 9: The reconstruction errors of three cases.

Case Case 1

Case 2

Case 3

405

406 407 408 409 410 411 412 413

method `1 `2 proposed `1 `2 proposed `1 `2 proposed

F1 1.6213 5.7863 0 13.8052 31.9876 0 6.0970 10.5291 0

F2 1.2606 9.3719 1.3986 18.6694 19.4095 8.1702 5.0548 12.5373 2.1297

F4 1.1206 7.6328 0 10.3437 19.1574 0 5.8798 14.4559 0

F5 1.1739 10.9413 1.0316 12.8158 18.8929 3.6842 9.9961 13.2722 3.4049

F6 1.1088 10.9959 0 8.0185 16.1334 0 6.0414 10.9240 0

F7 0.9054 3.0886 0 6.4539 11.5141 0 4.1594 10.3212 0

5. Experimental validation In this section, the proposed method is validated on a laboratory-scale cantilever beam. The experiment configuration is shown in Fig. 19. The steel cantilever beam has the dimension of 40 × 3 × 1cm3 . Four accelerometers are mounted on the midline of the beam evenly. A small mass is mounted on the free end. Three possible force locations are shown in the Fig. 19. Two hammers are used to excite the beam at the locations 1 and 2. The excitation force histories are measured by the force transducers. A Leuven Measurement & System (LMS) is used for data acquisition. The sampling frequency is 4096Hz. The system matrix is experimentally constructed according to [13, 33]. When a force of length nt excites the structure, the input-output system matrix can be express as     ˜ = H    

414

F3 0.8893 6.7078 0 12.0272 19.1637 0 5.9991 12.2778 0

h1 h2 h3 .. . hnt

 h1 h2 .. .

h1 .. .

..

···

h3

h2

       

0 .

(27)

h1

˜ We can reconstruct H ˜ by solving for the vector h = Note only nt independent variables exist in H. 23

Figure 19: The experiment configuration.

iT . Suppose that we excite the structure at the same location Ne times. We denote h1 h2 · · · hnt the response of a accelerometer as Xie and the excitation force as Fie (i = 1, 2, . . . , Ne ). The force vector can be rewritten as a matrix form. h

415 416 417

   i ¯ Fe =   

418

f1i f2i .. . fni t

0 f1i .. .

..

fni t −1

···

h

420 421 422 423 424 425 426 427 428 429 430 431 432 433 434

    

.

(28)

f1i

where fji is the jth element of vector Fie . And the vector h can be solved by h ∈ Rne ×1 = min

419



Ne

P

Xie − F ¯ ie h 2 2

(29)

i=1

In this paper, a classical Tikhonov regularization is considered for a stable solution to the Eq. (29). The measured responses and the calculated responses based on the identified system matrix are compared in Fig. 20. To validate the robustness of the proposed method, we add 30dB Gaussian noise into the responses. The assembled force history vector is divided into 20 blocks both in non-force detection and force reconstruction process. The non-force detection is carried out before force reconstruction. In this experiment, we draw 50000 samples of each variable to calculate the three indicators. The Markov chains of qi and Q - Q plot are shown in Figs. 21 and 22. The values of Q, means and variances of qi are listed in Table. 10. According to the non-force criterion, two forces act on the beam. The Markov chains of 4 kinds of monitored variables are shown in Fig. 23. The parallel Markov chains are considered well mixed. 50000 samples for each parallel chain are generated when 4 PSRFs are all below the threshold value 1.1. The second-half samples are concatenated to form a single chain of 75000 samples to produce the histogram of each variable shown in Fig. 24. The mean value, median value, and 95% credible interval of monitored parameters are listed in Table. 11. In this situation, the two impacts are exerted on the beam with different amplitudes and at time instants. The posterior mean, `1 regularization result, `2 regularization result, and 95% credible interval of F are shown in Fig. 25. The proposed method has a great reconstruction performance while traditional methods do not. 24

5

5 0

0 -5 -5 0

0.005

0.01

0.015

0.02

-10

0.025

0

5

0.005

0.01

0.015

0.02

0.025

0.015

0.02

0.025

10 5 0

0

-5 True Calculated

-10

-5 0

0.005

0.01

0.015

0.02

0.025

0

0.005

0.01

Figure 20: The measured responses and calculated responses based on the identified system matrix.

435 436

However, there are some misfits in the second force. This may be due to the imperfection of reconstructed system matrix. Neverless, the performance and robustness of the proposed method are validated in this experiment. 1.6 1.5

2

2

1.8

1.8

1.6

1.6

1.4

1.4

1.2

1.2

1.4 1.3 1.2 1.1 1

1 1

2

3

4

5

1

2

3

10 4

4

5

1 1

10 4

2

3

4

5 10 4

Figure 21: The Markov chains of qi at the three possible locations. 437

438

439 440 441 442 443 444

6. Conclusion In engineering practice, localization and reconstruction of external excitations is a tough problem due to the large problem dimension and lack of prior information on precision and shape parameters of forces. The selection of regularization parameters add additional complexity into the problem. In this paper, a hierarchical Bayesian formulation is proposed for both localization and reconstruction of forces in time domain. A Metropolis-within-Gibbs algorithm with nested blocking technique and a novel criterion for detecting non-force are also proposed to accelerate the identification procedure. One major merit of the 25

2

1.8

1.6

1.4

1.2

1

Quantiles of Input Sample

2

Quantiles of Input Sample

Quantiles of Input Sample

2

1.8

1.6

1.4

1.2

1 1

1.2

1.4

1.6

1.8

2

1.8

1.6

1.4

1.2

1 1

Quantiles of Uniform Distribution

1.2

1.4

1.6

1.8

2

Quantiles of Uniform Distribution

1

1.2

1.4

1.6

1.8

2

Quantiles of Uniform Distribution

Figure 22: The Q - Q plot of qi at all locations for the cantilever beam experiment. A blue line is the quantile-quantile line of this data. A dashed reference black line is shown by y = x.

Table 10: Values of indicators for non-force detection for the experiment case.

Loc Q M ean (q) V ar (q)

P1 0.6975 1.0115 0.0001

P2 0.6974 1.0121 0.0002

P3 0.9883 1.5600 0.0796

Figure 23: The Markov chains of kHF − Xk2 , τn , τf1 , τf2 , q1 , and q2 for the experiment.

26

Figure 24: The histograms of kHF − Xk2 , τn , τf1 , τf2 , q1 , and q2 for the experiment.

Table 11: Mean, median, and 95% CI of kHF − Xk2 , τn , τf1 , τf2 , q1 , and q2 for the sparse triangle force and continuous sine force.

M ean M edian 95%CI

q1 1.0098 1.0068 [1.0003,1.0369]

q2 1.0109 1.00745 [1.0003,1.0403]

5

20

0

15

-5

τn 33.9120 33.7973 [27.8438,40.6559]

τf1 0.6583 0.6560 [0.5208,0.8085]

τf2 0.6528 0.6498 [0.5157,0.8072]

kHF − Xk2 2.4321 2.4283 [2.2854,2.5997]

20

10

10 0

95% CI true force L1 regularization L2 regularization mean force

-10 -15

5 -10

0

-20 -5 0

0.005

0.01

0.015

0.02

0.025

0

0.005

0.01

0.015

0.02

0.025

-20 0

0.005

0.01

0.015

0.02

Figure 25: The estimated force histories of `1 , `2 regularization and proposed method with 95 % credible interval.

27

0.025

451

proposed method is that all to-be-determined parameters are determined in an automatic, date-driven manner, including force locations and shape parameters which have a significant impact on the reconstruction performance. By incorporating adaptivity into the identification method, its robustness is supposed to be greatly enhanced. Numerical and experimental validations on a cantilever beam under three different external excitation scenarios are carried out to validate the proposed method. Compared with the traditional regularization techniques, the results of the proposed method show better performance on both localization and force history reconstruction.

452

References

445 446 447 448 449 450

453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491

[1] B. Qiao, Z. Mao, J. Liu, Z. Zhao, X. Chen, Group sparse regularization for impact force identification in time domain, Journal of Sound and Vibration 445 (2019) 44–63. [2] M. Vigsø, J. C. Kristoffersen, R. Brincker, C. T. Georgakis, Indirect wave load estimates using operational modal analysis, Proceedings of the Twenty-eighth (2018) International Ocean and Polar Engineering Conference. [3] E. Jacquelin, A. Bennani, P. Hamelin, Force reconstruction: analysis and regularization of a deconvolution problem, Journal of Sound and Vibration 265 (1) (2003) 81–107. [4] A. N. Tikhonov, Regularization of incorrectly posed problems, in: Soviet Mathematics Doklady, Vol. 4, 1963, pp. 1624– 1627. [5] Z. Li, Z. Feng, F. Chu, A load identification method based on wavelet multi-resolution analysis, Journal of Sound and Vibration 333 (2) (2014) 381–391. [6] Y. Liu, W. S. Shepard Jr, Dynamic force identification based on enhanced least squares and total least-squares schemes in the frequency domain, Journal of Sound and Vibration 282 (1-2) (2005) 37–60. [7] D. Feng, M. Q. Feng, Identification of structural stiffness and excitation forces in time domain using noncontact visionbased displacement measurement, Journal of Sound and Vibration 406 (2017) 15–28. [8] R. Tibshirani, Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society: Series B (Methodological) 58 (1) (1996) 267–288. [9] B. Qiao, X. Zhang, C. Wang, H. Zhang, X. Chen, Sparse regularization for force identification using dictionaries, Journal of Sound and Vibration 368 (2016) 71–86. [10] S. Samagassi, A. Khamlichi, A. Driouach, E. Jacquelin, Reconstruction of multiple impact forces by wavelet relevance vector machine approach, Journal of Sound and Vibration 359 (2015) 56–67. [11] B. Qiao, J. Liu, J. Liu, Z. Yang, X. Chen, An enhanced sparse regularization method for impact force identification, Mechanical Systems and Signal Processing 126 (2019) 341–367. [12] B. Qiao, X. Zhang, J. Gao, X. Chen, Impact-force sparse reconstruction from highly incomplete and inaccurate measurements, Journal of Sound and Vibration 376 (2016) 72–94. [13] Q. Li, Q. Lu, Time domain force identification based on adaptive lq regularization, Journal of Vibration and Control 24 (23) (2018) 5610–5626. [14] M. Aucejo, Structural source identification using a generalized tikhonov regularization, Journal of Sound and Vibration 333 (22) (2014) 5693–5707. [15] J. Liu, X. Sun, X. Han, C. Jiang, D. Yu, A novel computational inverse technique for load identification using the shape function method of moving least square fitting, Computers & Structures 144 (2014) 127–137. [16] D. Feng, H. Sun, M. Q. Feng, Simultaneous identification of bridge structural parameters and vehicle loads, Computers & Structures 157 (2015) 76–88. [17] K. Maes, S. Gillijns, G. Lombaert, A smoothing algorithm for joint input-state estimation in structural dynamics, Mechanical Systems and Signal Processing 98 (2018) 292–309. [18] K. Worden, W. Staszewski, Impact location and quantification on a composite panel using neural networks and a genetic algorithm, Strain 36 (2) (2000) 61–68. [19] A. Rezayat, V. Nassiri, B. De Pauw, J. Ertveldt, S. Vanlanduit, P. Guillaume, Identification of dynamic forces using group-sparsity in frequency domain, Mechanical Systems and Signal Processing 70 (2016) 756–768. [20] K. Maes, K. Van Nimmen, E. Lourens, A. Rezayat, P. Guillaume, G. De Roeck, G. Lombaert, Verification of joint input-

28

492 493 494

[21]

495 496

[22]

497 498

[23]

499 500

[24]

501 502

[25]

503 504

[26]

505 506

[27]

507 508

[28]

509 510

[29]

511 512

[30]

513 514

[31]

515 516

[32]

517 518

[33]

519 520 521 522 523 524

[34] [35] [36] [37] [38]

525 526

[39]

527 528

[40]

529 530

[41]

531 532

[42]

533 534 535

[43]

536 537

[44]

538 539 540 541

[45]

state estimation for force identification by means of in situ measurements on a footbridge, Mechanical Systems and Signal Processing 75 (2016) 245–260. K. Maes, K. Van Nimmen, E. Lourens, A. Rezayat, P. Guillaume, G. De Roeck, G. Lombaert, Validation of joint input-state estimation for force identification and response estimation in structural dynamics. Q. Li, Q. Lu, Force localization and reconstruction using a two-step iterative approach, Journal of Vibration and Control 24 (17) (2018) 3830–3841. Q. Li, Q. Lu, Impact localization and identification under a constrained optimization scheme, Journal of Sound and Vibration 366 (2016) 133–148. V. K. Dertimanis, E. Chatzi, S. E. Azam, C. Papadimitriou, Input-state-parameter estimation of structural systems from limited output information, Mechanical Systems and Signal Processing 126 (2019) 711–746. E. Zhang, J. Antoni, P. Feissel, Bayesian force reconstruction with an uncertain model, Journal of Sound and Vibration 331 (4) (2012) 798–814. Q. Li, Q. Lu, A hierarchical bayesian method for vibration-based time domain force reconstruction problems, Journal of Sound and Vibration 421 (2018) 190–204. H. Sun, D. Feng, Y. Liu, M. Q. Feng, Statistical regularization for identification of structural parameters and external loadings using state space models, Computer-Aided Civil and Infrastructure Engineering 30 (11) (2015) 843–858. M. Aucejo, O. De Smet, On a full bayesian inference for force reconstruction problems, Mechanical Systems and Signal Processing 104 (2018) 36–59. M. Aucejo, O. De Smet, Bayesian source identification using local priors, Mechanical Systems and Signal Processing 66 (2016) 120–136. C. Faure, F. Ablitzer, J. Antoni, C. Pezerat, Empirical and fully bayesian approaches for the identification of vibration sources from transverse displacement measurements, Mechanical Systems and Signal Processing 94 (2017) 180–201. G. Yan, H. Sun, O. Büyüköztürk, Impact load identification for composite structures using bayesian regularization and unscented kalman filter, Structural Control and Health Monitoring 24 (5) (2017) e1910. Q. Li, Q. Lu, A revised time domain force identification method based on bayesian formulation, International Journal for Numerical Methods in Engineering 118 (7) (2019) 411–431. S. Atobe, S. Nonami, N. Hu, H. Fukunaga, Identification of impact force acting on composite laminated plates using the radiated sound measured with microphones, Journal of Sound and Vibration 405 (2017) 251–268. S. Boyd, L. Vandenberghe, Convex optimization, Cambridge University Press, 2004. M. Grasmair, Non-convex sparse regularisation, Journal of Mathematical Analysis and Applications 365 (1) (2010) 19–28. C. M. Bishop, Pattern recognition and machine learning, Springer, 2006. S. Brooks, A. Gelman, G. Jones, X.-L. Meng, Handbook of Markov chain Monte Carlo, CRC Press, 2011. J. Geweke, et al., Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments, Vol. 196, Federal Reserve Bank of Minneapolis, Research Department Minneapolis, MN, 1991. R. E. Kass, B. P. Carlin, A. Gelman, R. M. Neal, Markov chain monte carlo in practice: a roundtable discussion, The American Statistician 52 (2) (1998) 93–100. A. Gelman, D. B. Rubin, et al., Inference from iterative simulation using multiple sequences, Statistical Science 7 (4) (1992) 457–472. Z. Huang, A. Gelman, Sampling for bayesian computation with large datasets, Available at SSRN: https://ssrn.com/abstract=1010107. S. L. Scott, A. W. Blocker, F. V. Bonassi, H. A. Chipman, E. I. George, R. E. McCulloch, Bayes and big data: The consensus monte carlo algorithm, International Journal of Management Science and Engineering Management 11 (2) (2016) 78–88. C. Andrieu, A. Doucet, R. Holenstein, Particle markov chain monte carlo methods, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72 (3) (2010) 269–342. K. Maes, E. Lourens, G. De Roeck, G. Lombaert, General conditions for instantaneous system inversion in structural dynamics (2013) 17–18. K. Maes, E. Lourens, K. Van Nimmen, E. Reynders, G. De Roeck, G. Lombaert, Design of sensor networks for instantaneous inversion of modally reduced order models in structural dynamics, Mechanical Systems and Signal Processing 52 (2015) 628–644.

29

Highlights • •

•

•

A novel hierarchical Bayesian method for time domain force localization and reconstruction is proposed. The proposed method can adaptively determine the unknown parameters based on only the vibration responses. That is, the existence and sparsity/continuity of forces can all be determined in an automatic, date-driven manner. A Metropolis-within-Gibbs algorithm, companioned with nested blocking technique and a novel non-force detection criterion for sampling acceleration, is proposed for sampling the unknown parameters. Credible intervals for the identification are drawn, which provides a more comprehensive understanding of the identification results than traditional point estimation methods.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: