# Identification of nonlinear dynamical systems using approximate Bayesian computation based on a sequential Monte Carlo sampler

Mohamed-Anis Ben Abdessalem, Nikolaos Dervilis, David Wagg, Keith Worden

► To cite this version:

Mohamed-Anis Ben Abdessalem, Nikolaos Dervilis, David Wagg, Keith Worden. Identification of nonlinear dynamical systems using approximate Bayesian computation based on a sequential Monte Carlo sampler. ISMA, 2016, Leuven, Belgium. pp.2551-2566. hal-02527887

HAL Id: hal-02527887

https://univ-angers.hal.science/hal-02527887

Submitted on 26 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Identification of nonlinear dynamical systems using approximate Bayesian computation based on a sequential Monte Carlo sampler

A. Ben Abdessalem¹, N. Dervilis¹, D. Wagg¹, K. Worden¹

¹ Department of Mechanical Engineering, University of Sheffield, Mappin Street, Sheffield S1 3JD, United Kingdom

e-mail: a.b.abdessalem@sheffield.ac.uk

## Abstract

The Bayesian approach is well recognised in the structural dynamics community as an attractive approach to deal with parameter estimation and model selection in nonlinear dynamical systems. In the present paper, one investigates the potential of approximate Bayesian computation employing sequential Monte Carlo (ABC-SMC) sampling [1] to solve this challenging problem. In contrast to the classical Bayesian inference algorithms which are based essentially on the evaluation of a likelihood function, the ABC-SMC uses different metrics based mainly on the level of agreement between observed and simulated data. This alternative is very attractive especially when the likelihood function is complex and cannot be approximated in a closed form. Moreover, this flexibility allows one to use new features from either the temporal or the frequency domains for system identification. To demonstrate the practical applicability of the ABC-SMC algorithm, two illustrative examples are considered in this paper.

## 1 Introduction

The problem of *model selection* using *dynamic data* from a structural or mechanical system is widely en-
countered in mechanical and civil engineering. Generally speaking, model selection refers to the problem of
using data to select one model from a set of competing models [2]. In recent years, Bayesian approaches have
been widely used for parameter estimation and model selection in system identification of structures based
on measured data [3]. The usual approach is to find the optimal model amongst competing models which
can provide acceptable agreement between real and predicted data. It is obvious that a more complicated
model can fit the data better than a simple one which has fewer adjustable parameters. Therefore, if the op-
timal model is chosen by **minimising** some norm of the error between the output data and the corresponding
predictions of the **optimal model**, the optimal model will always be the most complicated one. However, the
selection of **complicated models** may lead to overfitting the data. When an overfitted model is used for future
prediction, it will **very likely** lead to poor results because the model will depend too much on the details of
the training data, and the noise in the data might have an important role in the data fitting. For this reason,
it is necessary to **penalise** a complicated model for robust and reliable predictions. For practical problems,
it means that **simpler models** are preferable to complicated ones, that is, the competing models should agree
closely with the **observed response** of the system and be simple as possible.

During the last few years, several papers have been published concerning model selection issues under a Bayesian framework in different domains and with various algorithms: structural dynamics [4–6], biology [7], cosmology [8] among others. In the field of structural dynamics and earthquake engineering, Beck and his colleagues [9–11] have pioneered the Bayesian model selection approach. Chatzi et al. [12] have proposed a Bayesian approach using the unscented Kalman filter (UKF) method to investigate the effects of
model complexity and parametrisation of nonlinear hysteretic systems. In their contribution, the authors examined the case when not only the parameters of the system are unknown but also the nature of the analytical model describing the system. Worden & Hensman [13] have attempted to solve the problem of Bouc-Wen (BW) hysteretic system parameter identification as a model selection problem using the Metropolis-Hastings (MH) based Markov chain Monte Carlo (MCMC) algorithm. In the Bayesian framework and in particular for model selection issue, several selection criterions have been proposed in the literature to rank the competing models, essentially the Akaike information criterion (AIC), the Bayesian information criterion (BIC) and the cross-validation criterion. For more details concerning the mathematical formulations and the significance of these metrics, the reader can refer to the following references [14, 15].

In this paper, one proposes to introduce the approximate Bayesian computation using sequential Monte Carlo (ABC-SMC) method to select the most plausible dynamic models representing a mechanical system by using its response measurements. During the last decade, ABC is increasingly considered as a distinct inferential formalism and not merely as an approximation to conventional Bayesian inference. This algorithm was originally developed to solve problems in biology [1]; however, as its efficiency was recognised, its popularity increased and, consequently it started to be used in different domains and recently in structural dynamics [16]. The motivation behind the use of ABC instead of the classical MCMC methods is that the latter often have convergence problems. Moreover, one of the important merits of this technique compared to the classical MCMC algorithms is its independance of the initial parameters, which should be specified in MCMC algorithms. This has led to the development of alternatives to MCMC (e.g., [17]) and a group of these alternative methods has come to be known as ABC. Different variants have been proposed in the literature, for more details the reader can refer to [18–20]. Recently, Chiachio et al. [21] proposed a new algorithm, named ABC-Subsim, by incorporating the subset simulation technique [22] into the ABC algorithm. In practice, and as shown in several published papers, the ABC algorithm presents several advantages; it does not require the definition of a likelihood function, therefore a considerable amount of data can be analysed and more complex models can be used. Compared with the classical MCMC algorithms which rely on the use of the Bayesian evidence to rank the competing models, the ABC algorithm is straightforward and converges directly to the most likely model by eliminating progressively the least good models. In fact, the ABC-SMC algorithm proposes jumps between the models with different dimensionalities without matching the dimensionalities. So it is simpler in a sense than the rejection-jump MCMC proposed by Green [23]. In addition, the ABC algorithm offers the possibility to use different metrics for model selection and parameter estimation, which is a major advantage. As far as the authors are aware, there has been relatively little work on ABC-SMC algorithm for Bayesian model selection in structural dynamics so far. To demonstrate the effectiveness of the ABC-SMC algorithm to solve challenging problems, two illustrative examples using simulated data will be presented in this work. The first example presents cubic and cubic-quintic models, while the second example deals with the identification of the Bouc-Wen hysteretic system formulated as a model selection task.

The rest of this paper is organised as follows. In Section 2, a brief description of the Bayesian formulation for parameter inference and model selection is presented. In Section 3, an overview of the basic elements of the ABC algorithm and the algorithm settings are given. In Section 4, two illustrative examples in structural dynamics are then presented and investigated. Finally, Section 5 contains concluding remarks and future work.

## 2 General Bayesian framework

The Bayesian approach offers a suitable framework for parameter inference and model selection in dynamical systems when multiple models can fit the data in a reasonable manner. Denote by $\mathcal{D}$ the data from a structural or mechanical system; the goal is to use $\mathcal{D}$ to select the most plausible model representing the system out of $n$ given models $M = \{M_1, M_2, \dots, M_n\}$. The basic idea of the Bayesian approach to identification is that, by repeatedly applying Bayes' theorem, one can assess the probability of a set of parameters $\theta$ as well
as a model structure $\mathcal{M}$ conditional on the data $\mathcal{D}$ using:

$$p(\theta|\mathcal{D}, \mathcal{M}) = \frac{p(\mathcal{D}|\theta, \mathcal{M})p(\theta|\mathcal{M})}{p(\mathcal{D}|\mathcal{M})} \qquad (1)$$

and

$$p(\mathcal{M}|\mathcal{D}) = \frac{p(\mathcal{D}|\mathcal{M})p(\mathcal{M})}{p(\mathcal{D})} \qquad (2)$$

where

$$p(\mathcal{D}|\mathcal{M}) = \int p(\mathcal{D}|\theta, \mathcal{M})p(\theta|\mathcal{M})d\theta \qquad (3)$$

is a normalising constant which ensures that $p(\theta|\mathcal{D}, \mathcal{M})$ integrates to unity. This is referred as the model evidence. In Eqs. (1) and (2), $p(\theta|\mathcal{D}, \mathcal{M})$ is the posterior distribution, $p(\mathcal{D}|\theta, \mathcal{M})$ is the likelihood function which quantifies the probability of getting $\mathcal{D}$; $p(\theta|\mathcal{M})$ is the prior probability density function (PDF) assigned to model parameter values $\theta$; $p(\mathcal{M}|\mathcal{D})$ is the posterior probability of $\mathcal{M}$ and $p(\mathcal{M})$ is the prior probability assigned to $\mathcal{M}$.

# 3 Approximate Bayesian computation

## 3.1 Basic features

The ABC-SMC algorithm provides a unified approach to deal with both parameter estimation and model selection in a straightforward way. The algorithm is similar to the particle filter technique which have been widely used for system identification. The implementation of the algorithm requires the definition of the prior distributions associated to the model parameters and the competing models. In addition, the algorithm needs to define an appropriate sequence of the tolerance schedule ($\epsilon_1 > \dots > \epsilon_T$) to guarantee a smooth transition from the prior to the posterior parameter distributions. Finally, a kernel functions should be selected to move a candidate sample in the search space. The perturbation allows a better exploration with the intention of generating better samples. It should be also noted that the algorithm requires the definition of an appropriate metric to evaluate the agreement between observed and simulated data. To implement the ABC-SMC method for model selection, the mentioned parameters should be carefully selected for the algorithm to work well. A brief description is given here of the ABC-SMC algorithm adopted in this work. Let assume $n$ competing models $\mathcal{M}_{i=1:n}$ and each model is defined by its parameter vector $\Theta_{i=1:n}$. Based on the knowledge of the first tolerance schedule value and the selected metric, the algorithm starts by sampling the samples $\{(\mathcal{M}_1, \Theta_1), \dots, (\mathcal{M}_n, \Theta_n)\}$ from the prior distributions until $N$ samples are obtained. In the first iteration, one assigns the same weight for all the accepted samples. The samples are then perturbed to find better ones and propagated through a sequence of intermediate distributions. It should be noted that the samples are reweighted in the subsequent iterations. Finally, based on the relative acceptance frequencies for the different competing models, one may approximate the posterior model probabilities for all the competing models simultaneously. Further details can be found in [1].

## 3.2 Algorithm settings

By using the ABC-SMC algorithm, there are two practical problems to be solved. First, one has to choose a suitable metric which reflects the level of agreement between the exact measurements and the predicted ones. Then, one needs to define a sequence of tolerance thresholds to converge gradually to the posterior distributions. In this work, the normalised mean square error (MSE) given by Eq. (4) is selected as a metric.
$$h = \frac{100}{N\sigma_u^2} \sum_{\ell=1}^{N} (u_{\ell}^{\text{obs}} - \tilde{u}_{\ell}^{\text{sim}})^2 \quad (4)$$

where N is the number of samples, $\sigma_u^2$ is the variance of the measured displacement, $u_{\ell}^{\text{obs}}$ and $\tilde{u}_{\ell}^{\text{sim}}$ are the observed and predicted displacements given by the selected model $\mathcal{M}(\cdot)$, respectively.

For the present examples, the tolerance schedule required for ABC inference is defined a priori for demonstrative purpose. It should be mentioned that this can be defined adaptively based on the values of the selected metric at the previous population for instance. By decreasing $\varepsilon$, one is inducing a gradual transition from the prior to the posterior parameter distributions. In practice, the final tolerance schedule should tend to zero; however, this choice is highly case dependent. The computational time associated to the ABC-SMC depends mainly on the target tolerance $\varepsilon_T$ and the population size. This means that the selection of $\varepsilon_T$ has a large influence on the posterior parameter estimates and their associated posterior uncertainties. Often, the choice of $\varepsilon_T$ should reflect a trade-off between computability and accuracy. To demonstrate the usefulness and the applicability of the ABC-SMC algorithm for model selection in structural dynamics, two illustrative examples detailed in the following section will be investigated.

# 4 Numerical examples

## 4.1 Example 1: cubic and cubic-quintic models

The first example aims to assess the capability of the ABC-SMC algorithm in the case when two competing models are proposed. The competing models proposed in this example are the cubic $\mathcal{M}_1$ and the cubic-quintic $\mathcal{M}_2$ given by Eqs. (5) and (6), respectively:

$$\mathcal{M}_1 : m\ddot{y} + c\dot{y} + ky + k_3y^3 = f(t), \quad (5)$$

$$\mathcal{M}_2 : m\ddot{y} + c\dot{y} + ky + k_3y^3 + k_5y^5 = f(t). \quad (6)$$

where m represents the mass, c the damping, k linear stiffness, $k_3$ and $k_5$ are the non-linear stiffness coefficients. $y$, $\dot{y}$ and $\ddot{y}$ are displacement, velocity and acceleration responses of the system, respectively. The forcing function is of the form $f(t) = F \cos(\omega t)$, with $\omega = 20 \text{rad/s}$ and $F = 10$.

For both models, the unknown parameters are assumed to be uniformly distributed over a sufficiently wide range (see, Table 1). In this example, the training data used for model selection were generated from the cubic model. Denote by $\Theta = \{m, c, k, k_3\}$ the parameters associated to the model to be identified. Typically, one sets the prior probabilities of each model to be equal, i.e., $\text{Pr}(\mathcal{M}_1) = \text{Pr}(\mathcal{M}_2) = \frac{1}{2}$. The duration of measurements is $T = 10$ s with sampling frequency, $f_{\text{res}} = 200$ MHz, so that the number of data points is $N = 1000$. The response time history of the correct model is calculated by the Runge-Kutta method and shown in Fig. 1. In the first case, one assumes that the data are free of noise. Then, the case when measurements are corrupted by noise is considered.

<table><caption>Table 1: Parameter ranges of competing models.</caption><thead><tr><th>Parameter</th><th>True value</th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>m</td><td>1</td><td>0.1</td><td>10</td></tr><tr><td>c</td><td>0.05</td><td>0.005</td><td>0.5</td></tr><tr><td>k</td><td>50</td><td>5</td><td>500</td></tr><tr><td>k<sub>3</sub></td><td>10<sup>3</sup></td><td>10<sup>2</sup></td><td>10<sup>4</sup></td></tr><tr><td>k<sub>5</sub></td><td>10<sup>4</sup></td><td>10<sup>3</sup></td><td>10<sup>5</sup></td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1: The response time history of the cubic model.</figcaption></figure>

In this example, the initial tolerance was estimated through some numerical simulations and fixed at $\epsilon_1 = 100$. The prior distributions are taken to be uniform and the perturbation kernel $K_\theta$ is updated from population to population, $K_\theta = \xi U(-\delta, \delta)$, with $\xi = 0.5$. For the t-th population, $\delta$ is estimated based on the support of the posterior distribution from the $(t-1)$-th population given by:

$$ \delta = \max(\theta)_{t-1} - \min(\theta)_{t-1} \qquad (7) $$

In the present work, the population size (or the number of samples) is fixed at $N = 2000$ which seems reasonable to achieve precise estimation of the posterior distributions and acceptable computational time. To ensure the gradual transition between prior and posterior distributions, $t = 18$ populations with $\epsilon_t = \{100, \dots, 0.03\}_{t=1:18}$ are generated.

## 4.2 Results and discussion

Fig. 2 shows the marginal posterior distributions for the competing models at some intermediate populations. One may observe how the ABC-SMC algorithm converges gradually to the correct model. At population 17, $M_2$ is eliminated which means that the model is unable to describe the system behaviour given the available measurements. Model probabilities associated to each competing model at intermediate selected populations are also shown in Fig. 2. Fig. 3 shows the histograms of the unknown parameters associated to the correct model. One observes that the obtained histograms are highly peaked around the true values given in Table 1. The statistics associated to the inferred parameters are summarised in Table 2.

Fig. 4 shows the simulated and predicted responses using the mean values given in Table 2. A good agreement between simulated and predicted responses is observed which demonstrates the efficiency of the ABC-SMC inference method. Fig. 5 shows the MSE distribution associated to the last population. Note that the tolerance threshold can be decreased further. Decreasing the tolerance threshold further reduces the uncertainties in the model parameters and then improves the model predictions; however, this may increase the computational time drastically. For the present case, one assumes that the final level of accuracy is enough to get robust prediction. Fig. 6 shows the scatter plots of the model parameters which show how samples move from population to population. It should be noted that the number of samples associated to the correct model increases by decreasing the tolerance threshold which allows us to get a precise estimate of the uncertainties associated to each unknown parameter.

## 4.3 Model selection using corrupted measurements

To investigate how noise may impact the model selection as well as the model parameters, the exact measurements are now corrupted by noise. Here, noise of RMS 1% of the response is added to the displacement
<figure><img src="image_3.png" /><figcaption>Figure 2: Marginal model posterior distributions at some selected populations.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Histograms of the approximate posterior distributions of *m*, *c*, *k* and *k*<sub>3</sub>.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Simulated and predicted displacement responses.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: MSE distribution at the last population.</figcaption></figure>

<table><caption>Table 2: Identification results using simulated data.</caption><thead><tr><th>Parameter</th><th>Exact value</th><th>Estimated value</th><th>σ</th><th>μ - 3σ</th><th>μ + 3σ</th></tr></thead><tbody><tr><td>m</td><td>1</td><td>1.0002</td><td>0.0068</td><td>0.9799</td><td>1.0205</td></tr><tr><td>c</td><td>0.05</td><td>0.05</td><td>0.0034</td><td>0.0399</td><td>0.06009</td></tr><tr><td>k</td><td>50</td><td>50.0999</td><td>1.7784</td><td>44.7647</td><td>55.4350</td></tr><tr><td>k<sub>3</sub></td><td>1000</td><td>1000.8215</td><td>24.5139</td><td>927.2798</td><td>1074.3632</td></tr></tbody></table>

signal. Table 3 shows the optimal parameters and how those parameters change compared with the previous case. It should be mentioned that we keep the same sequence of the tolerance schedule used previously. For the noisy case, the algorithm is stopped at the tolerance threshold $ε_{17} = 0.05$ as the acceptance rate drops significantly and the algorithm cannot move efficiently. Using the noisy measurements, ABC-SMC algorithm shows that model $M_1$ is the most probable model based on the data, with a posterior probability of 1 at population 17. Table 3 summarises the statistics associated with the optimal parameters using noisy measurements. Larger standard deviations are obtained compared with the results with exact measurements.

<table><caption>Table 3: Identification results using noisy simulated data.</caption><thead><tr><th>Parameter</th><th>Exact value</th><th>Estimated value</th><th>σ</th><th>μ - 3σ</th><th>μ + 3σ</th></tr></thead><tbody><tr><td>m</td><td>1</td><td>1.0005</td><td>0.0080</td><td>0.9763</td><td>1.0246</td></tr><tr><td>c</td><td>0.05</td><td>0.0502</td><td>0.0042</td><td>0.0376</td><td>0.0627</td></tr><tr><td>k</td><td>50</td><td>50.0796</td><td>2.2470</td><td>43.3385</td><td>56.8208</td></tr><tr><td>k<sub>3</sub></td><td>1000</td><td>1001.8587</td><td>28.9190</td><td>915.1016</td><td>1088.6157</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6: Two-dimensional scatterplots of the samples at some intermediate populations.</figcaption></figure>

## 4.4 Example 2: the Bouc-Wen hysteresis model

This case study follows the spirit of the second example in Worden and Hensman [13], where the Metropolis Hastings algorithm was used to identify the Bouc-Wen (BW) model parameters. The highly nonlinear nature of the BW model, along with a large number of model parameters has made the identification of this system a challenging problem. The model has been receiving more attention in recent times due to the development of efficient numerical algorithms that can be used to identify the model parameters more accurately. Several methods have been discussed in the literature to accomplish this, including analytic approaches, Gauss-Newton, simplex, reduced gradient method, and extended Kalman filters. Due to the nonlinear nature of the problem, stochastic optimisation algorithms have been also found to be well suited. Algorithms such as genetic algorithm [24], particle swarm optimisation [25], and differential evolution [26] have been successfully implemented. However, those algorithms do not provide a probabilistic characterisation of the uncertainty of the unknown parameters. In this paper, one investigates the capability of the ABC-SMC approach to deal with the BW model parameters identification. The identification problem is formulated as a model selection task which allows one to reduce the number of parameters.
<figure><img src="image_8.png" /><figcaption>Figure 7: Simulated data from the BW model, $n = 2$.</figcaption></figure>

### 4.4.1 Equations of motion

The general single-degree-of-freedom (SDOF) hysteretic system described in the terms of Wen [27], is represented below where $g(y, \dot{y})$ is the polynomial part of the restoring force, $z(y, \dot{y})$ the hysteretic part and $x(t)$ is the excitation force:

$$m\ddot{y} + g(y, \dot{y}) + z(y, \dot{y}) = x(t) \quad (8)$$

The hysteretic component is defined by Wen [27] via the additional equation of motion:

$$\dot{z} = \begin{cases} -\alpha|\dot{y}|z^n - \beta\dot{y}|z^n| + A\dot{y}, & \text{for } n \text{ odd} \\ -\alpha|\dot{y}|z^{n-1}|z| - \beta\dot{y}|z^n| + A\dot{y}, & \text{for } n \text{ even} \end{cases} \quad (9)$$

The parameters $\alpha, \beta$ and $n$ govern the shape and the smoothness of the hysteretic loop. The equations offer a simplification from the point of view of parameter estimation in that the stiffness term in (9) can be combined with the $A\dot{y}$ term in the state equation for $z$. The reader can refer to [13] for further details.

In this example, the response variable will be assumed to be displacement. The sampling interval was taken as 0.001 s, corresponding to a sampling frequency of 1000 Hz. As in the first example, the problem of model selection is considered by using both exact and corrupted measurements. The training data used here were composed of 10000 points corresponding to a record duration of 10 seconds. Fig. 7 shows the BW model response; the exact values used to generate the training data and the parameters ranges are summarised in Table 4 with $n = 2$.

<table><caption>Table 4: Parameters ranges of the BW model.</caption><thead><tr><th>Parameter</th><th>True value</th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>m</td><td>1</td><td>0.1</td><td>10</td></tr><tr><td>c</td><td>20</td><td>2</td><td>200</td></tr><tr><td>α</td><td>1.5</td><td>0.15</td><td>15</td></tr><tr><td>β</td><td>-1.5</td><td>-15</td><td>-0.15</td></tr><tr><td>A</td><td>6680</td><td>668</td><td>66800</td></tr></tbody></table>

To deal with the problem of parameter estimation and model selection, 5 competing models are considered, the linear model and 4 models by varying $n$ from 1 to 4. The competing models are defined as follows:

$$\left\{ \begin{array}{l} M_1 : m\ddot{y} + c\dot{y} + Ay = x(t) \\ M_2 : \text{Eqs. (8) and (9)}, \quad n = 1 \\ M_3 : \text{Eqs. (8) and (9)}, \quad n = 2 \\ M_4 : \text{Eqs. (8) and (9)}, \quad n = 3 \\ M_5 : \text{Eqs. (8) and (9)}, \quad n = 4 \end{array} \right.$$
<figure><img src="image_9.png" /><figcaption>Figure 8: Model posterior probabilities at different intermediate populations.</figcaption></figure>

To implement the ABC-SMC algorithm, equal prior plausibilities $\text{Pr}(\mathcal{M}_i|\mathcal{D}) = \frac{1}{5}$ are considered. As in the first example, the sequence of the tolerance is defined a priori $\varepsilon_t = \{100, \dots, 5 \times 10^{-8}\}_{t=1:44}$. The obtained results are summarised in the following subsections for the exact and corrupted measurements.

### 4.4.2 Results and analyses

Fig. 8 shows the posterior model probabilities at different intermediate populations. One observes how the ABC-SMC algorithm proposed here allows a progressive elimination of some competing models during the inference procedure. After some intermediate populations, the ABC-SMC approach identifies the true model from which the training data were generated. As one can see from Fig. 8, the BW model with $n = 4$ is firstly eliminated, then the linear model, which is not able to capture the nonlinearity of the BW model. The ABC-SMC algorithm identifies the correct model at population 32 as shown in Fig. 8 at the tolerance threshold $\varepsilon = 10^{-5}$. It should be mentioned that once a model is eliminated, the algorithm continues the inference by using only the remaining competing models. Therefore from population 33 to population 44, the algorithm uses only $\mathcal{M}_3$ and improves the precision related to the posterior distributions of the model parameters. In the present case, the target tolerance is equal to $5 \times 10^{-8}$ which is sufficient to get a precise estimate of the true parameters. For the present case, 44 populations are used to converge to the correct model since several competing models are proposed and the model responses are quite similar.

Table 5 summarises the obtained results from the ABC-SMC algorithm. One may observe that the algorithm provides a precise estimate of the true values. The posterior distribution associated to each parameter is shown in Fig. 9. As one can see from this example, to converge to the true posterior distribution and to gain more confidence to the estimated parameters, the final tolerance threshold should tend towards zero. Fig. 10 shows the simulated and predicted responses; a good agreement is shown, which proves the effectiveness of the ABC-SMC method to deal with parameter estimation and model selection issues.

Using the samples obtained at the last population, the correlation matrix can be computed. As one can see below, a positive correlation exists between $(A, \beta)$ while a negative correlation between $(\alpha, c)$ is observed.
<table><caption>Table 5: Estimated parameters: exact measurements, target $ε_T = 5 \times 10^{-8}$.</caption><thead><tr><th>Parameter</th><th>Exact value</th><th>Estimated value</th><th>σ</th><th>μ - 3σ</th><th>μ + 3σ</th></tr></thead><tbody><tr><td>m</td><td>1</td><td>0.9999</td><td>0.0011</td><td>0.9967</td><td>1.0032</td></tr><tr><td>c</td><td>20</td><td>19.9995</td><td>0.0217</td><td>19.9344</td><td>20.0645</td></tr><tr><td>α</td><td>1.5</td><td>1.5000</td><td>0.0012</td><td>1.4965</td><td>1.5035</td></tr><tr><td>β</td><td>-1.5</td><td>-1.5000</td><td>0.0066</td><td>-1.5198</td><td>-1.4802</td></tr><tr><td>A</td><td>6680</td><td>6680.0030</td><td>0.1744</td><td>6679.4797</td><td>6680.5262</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Figure 9: Histograms of the BW model parameters (exact measurements).</figcaption></figure>

<table><caption>Correlation of BW model parameters</caption><thead><tr><th></th><th>m</th><th>c</th><th>α</th><th>β</th><th>A</th></tr></thead><tbody><tr><th>m</th><td>1.0000</td><td>-0.0126</td><td>-0.0296</td><td>0.0311</td><td>0.0819</td></tr><tr><th>c</th><td>-0.0126</td><td>1.0000</td><td><strong>-0.4962</strong></td><td>0.0810</td><td>-0.1110</td></tr><tr><th>α</th><td>-0.0296</td><td><strong>-0.4962</strong></td><td>1.0000</td><td>0.0586</td><td>-0.0347</td></tr><tr><th>β</th><td>-0.0311</td><td>0.0810</td><td>0.0586</td><td>1.0000</td><td><strong>0.9469</strong></td></tr><tr><th>A</th><td>0.0819</td><td>-0.1110</td><td>-0.0347</td><td><strong>0.9469</strong></td><td>1.0000</td></tr></tbody></table>

Fig. 11 shows the evolution of the samples at some intermediate populations as the mean values estimated from the last population. One may observe the capability of the ABC-SMC algorithm to converge towards the true values. One may also observe how by decreasing the threshold tolerance, the samples move towards a small area with higher probability. It should be mentioned that the correlation mentioned previously between $(A, β)$ and $(α, c)$ can be clearly noticed.
<figure><img src="image_11.png" /><figcaption>Figure 10: Simulated and predicted responses.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Two-dimensional scatterplots of the posterior distributions.</figcaption></figure>

### 4.4.3 Model selection using noisy measurements

In this subsection, one examines the case when measurements are corrupted with different levels of noise. Table 6, shows the predicted results in the case when the exact measurements are corrupted by RMS 0.1%. In this case the final tolerance threshold is equal to $3.8 \times 10^{-6}$. From Table 6, one may observe that parameters $(m, \beta, A)$ are not affected by the noise, however, $(c, \alpha)$ are clearly affected. As one can see, a decrease of $c$ is compensated by an increase of $\alpha$. At the final tolerance threshold, the following parameters $(m, \beta, A)$ converge very close to their target values.
<table><caption>Table 6: Estimated parameters, noisy measurements: RMS 0.1%, target $ε_T = 5 \times 10^{-6}$.</caption><thead><tr><th>Parameter</th><th>Exact value</th><th>Estimated value</th><th>σ</th><th>μ - 3σ</th><th>μ + 3σ</th></tr></thead><tbody><tr><td>m</td><td>1</td><td>0.9991</td><td>0.0063</td><td>0.9801</td><td>1.0180</td></tr><tr><td>c</td><td>20</td><td>18.6221</td><td>0.1354</td><td>18.2160</td><td>19.0283</td></tr><tr><td>α</td><td>1.5</td><td>1.6397</td><td>0.0071</td><td>1.6183</td><td>1.6610</td></tr><tr><td>β</td><td>-1.5</td><td>-1.5050</td><td>0.0419</td><td>-1.6306</td><td>-1.3794</td></tr><tr><td>A</td><td>6680</td><td>6680.0078</td><td>1.0864</td><td>6676.7487</td><td>6683.2669</td></tr></tbody></table>

To better understand how noise may affect the estimated parameters, one increases the noise level which is equal now to RMS 1%. Table 7 shows the obtained results. The same tendency observed previously is shown here, by increasing the noise level, *c* parameter continues to decrease while *α* continues to increase.

<table><caption>Table 7: Estimated parameters, noisy measurements: RMS 1%, target $ε_T = 6 \times 10^{-4}$.</caption><thead><tr><th>Parameter</th><th>Exact value</th><th>Estimated value</th><th>σ</th><th>μ - 3σ</th><th>μ + 3σ</th></tr></thead><tbody><tr><td>m</td><td>1</td><td>0.9874</td><td>0.0148</td><td>0.9430</td><td>1.0318</td></tr><tr><td>c</td><td>20</td><td>9.0158</td><td>0.3694</td><td>7.9077</td><td>10.1239</td></tr><tr><td>α</td><td>1.5</td><td>2.7849</td><td>0.0253</td><td>2.7090</td><td>2.8608</td></tr><tr><td>β</td><td>-1.5</td><td>-1.5966</td><td>0.1714</td><td>-2.1131</td><td>-1.0847</td></tr><tr><td>A</td><td>6680</td><td>6678.0497</td><td>4.4828</td><td>6665.2013</td><td>6692.0981</td></tr></tbody></table>

# 5 Conclusion

This paper shows the ability of the ABC-SMC algorithm to solve efficiently the issue of model selection in structural dynamics in a straightforward way. The ability of the ABC-SMC is demonstrated through two nonlinear examples and the obtained results seem very promising to extend the use of the ABC-SMC algorithm to more complex problems. The flexibility offered by the ABC algorithm is very important as it allows to formulate new metrics from the frequency domain (power spectra, frequency peaks location, ...) for instance. Therefore, the ABC-SMC opens up an alternative paradigm for solving model selection/system identification issues using new metrics to evaluate the agreement between measured and simulated data. To reduce the computational requirements, one aims to implement the proposed scheme in a fully parallel environment to effectively distribute the different computational tasks in available multi-core CPUs. Next, application of the ABC-SMC algorithm to data from real structural systems is currently under investigation.

# Acknowledgements

The support of the UK Engineering and Physical Sciences Research Council (EPSRC) through grant reference no. EP/K003836/1 is greatly acknowledged. The authors would like to thank James Hensman of the University of Lancaster for introducing them to the ABC-SMC method, and for a number of useful discussions.

# References

[1] T. Toni, D. Welch, N. Strelkowa, A. Ipsen, M.P.H. Stumpf. *Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems*. Journal of The Royal Society Interface (2009) 6, 187202.
[2] J.L. Beck, K.V. Yuen. *Model Selection Using Response Measurements: Bayesian Probabilistic Approach*, Journal of Engineering Mechanics, 130 (2004).

[3] P.L. Green, *Bayesian system identification of a non linear dynamical system using a novel variant of Simulated Annealing*, Mechanical Systems and Signal Processing 52-53(2015), 133-146.

[4] K.V. Yuen, *Bayesian Methods for Structural Dynamics and Civil Engineering*, Wiley: New York, NY, 2010.

[5] P.L. Green, E.J. Cross, K. Worden, *Bayesian system identification of dynamical systems using highly informative training data*, Mechanical Systems and Signal Processing 56-57(2015), 109-122.

[6] P.L. Green, K. Worden, *Bayesian and Markov chain Monte Carlo methods for identifying nonlinear systems in the presence of uncertainty*, Philosophical Transactions of the Royal Society. A373: 20140405.

[7] T. Toni & M.P.H. Stumpf, *Simulation-based model selection for dynamical systems in systems and population biology*, Bioinformatics. Vol. 26 no. 1 2010, pages 104-110.

[8] E.E.O. Ishida, S.D.P. Vitenti, M. Penna-Lima, J. Cisewski, RS. de Souza, A.M.M. Trindade, E. Cameron, VC. Busti, *cosmoabc: Likelihood-free inference via Population Monte Carlo Approximate Bayesian Computation*, Astronomy and Computing 13 (2015), 1-11.

[9] J.L. Beck, *Bayesian system identification based on probability logic*, Structural Control and Health Monitoring. 17 (2010), 825-847.

[10] S.H. Cheung, J.L. Beck, *Calculation of posterior probabilities for Bayesian model class assessment and averaging from posterior samples based on dynamic system data*, Computer-Aided Civil and Infrastructure Engineering. 25 (2010), 304-321.

[11] J. Ching, J. Beck, K. Porter, *Bayesian state and parameter estimation of uncertain dynamical systems*, Probabilistic Engineering Mechanics. 21 (2006), 81-96.

[12] E.N. Chatzi, A.W. Smyth, S.F. Masri, *Experimental application of on-line parametric identification for nonlinear hysteretic systems with model uncertainty*, Structural Safety 32 (2010), 326-337.

[13] K. Worden, J.J. Hensman, *Parameter estimation and model selection for a class of hysteretic systems using Bayesian inference*, Mechanical Systems and Signal Processing, 32(2012), 153-169.

[14] W. Zucchini, *An Introduction to Model Selection*, Journal of Mathematical Psychology 44, (2000), 41-61.

[15] L. Wasserman, *Bayesian Model Selection and Model Averaging*. Journal of Mathematical Psychology 44 (2000), 92-107.

[16] M.K. Vakilzadeh, Y. Huang, J.L. Beck, T. Abrahamsson, *Approximate Bayesian Computation by Subset Simulation using hierarchical state-space models*, Mechanical Systems and Signal Processing.(2016), http://dx.doi.org/10.1016/j.ymssp.2016.02.024i.

[17] P. Del Moral, A. Doucet, A. Jasra, *Sequential Monte Carlo samplers*. Journal of the Royal Statistical Society: Series B, 68, (2006), 411-436.

[18] M. Beaumont, J. Cornuet, J. Marin, C. Robert, *Adaptive approximate Bayesian computation*, Biometrika 96 (4) (2009), 983-990.

[19] M. Baragatti, A. Grimaud, D. Pommeret, *Likelihood-free parallel tempering*, Statistics and Computing 23 (4) (2013), 535-549.
[20] S. Sisson, Y. Fan, M. Tanaka, A note on backward kernel choice for sequential Monte Carlo without likelihoods, Tech. Rep., University of New South Wales, 2009.

[21] M. Chiachio, J.L. Beck, J. Chiachio, G. Rus, Approximate Bayesian computation by Subset Simulation, SIAM Journal on Scientific Computing, 36(3):A1339-A1338, June 2014

[22] S.K. Au, J.L Beck, *Estimation of small failure probabilities in high dimensions by Subset Simulation*, Probabilistic Engineering Mechanics, 16(2001), 263-277.

[23] P.J. Green, *Reversible jump Markov chain Monte Carlo computation and Bayesian model determination*, Biometrika, (1995) 82(4),711732.

[24] K. Chwastek, J. Szczyglowski, *Identification of a hysteresis model parameters with genetic algorithms*, Mathematics and Computers in Simulation 71 (2006), 206-211.

[25] S. Xiao, YangminLi, *Dynamic compensation and $H_\infty$ control for piezoelectric actuators based on the inverse Bouc-Wen model*, Robotics and Computer-Integrated Manufacturing 30(2014) 47-54.

[26] A. Kyprianou, K. Worden, *Identification of hysteretic systems using the differential evolution algorithm*, Journal of Sound and Vibration (2001) 248(2), 289-314.

[27] Y. Wen, *Method for random vibration of hysteretic systems*, ASCE Journal of the Engineering Mechanics Division 102(2), April 1976, 249-263.