Bayesian and Markov chain Monte Carlo methods for identifying nonlinear systems in the presence of uncertainty

1. Introduction

System Identification (SI) is a technique of considerable importance within the discipline of structural dynamics. In the absence of a complete physics-based description of a system or structure, SI can provide the missing pieces of information that allow the formulation of a descriptive or predictive model. When the structure of interest has linear dynamical behaviour, the problem of SI is well established, to the extent that authoritative text books and monographs exist [1,2]. In the case of linear dynamical systems, it is usually sufficient to consider sets of linear second-order differential equations (modal models) or first-order differential equations (state-space models) as the appropriate mathematical model structure. In that case, the SI problem is largely reduced to determining the correct number of equations and the numerical parameters in the model. Unfortunately, most structures will, in reality, display nonlinear characteristics to some extent, and the SI problem for nonlinear structures and systems is by no means solved. One of the main problems in nonlinear SI is the number and variety of possible model structures that arise once the variety of possible nonlinearities is taken into account [3,4].

It is not necessary here to provide a detailed classification of nonlinear SI models and approaches; however, it will prove useful to give a higher level breakdown of model structures based on their motivation. Predictive models can be divided into three classes: white, grey and black-box models.

• White-box models are taken here to be those whose equations of motion have been derived completely from the underlying physics of the problem of interest and in which the model parameters have direct physical meanings. Finite-element models constitute one sub-class of such models.

• Black-box models are, by contrast, usually formed by adopting a parametrized class of models with some universal approximation property and learning the parameters from measured data; in such a model, like a neural network, the parameters will not generally carry any physical meaning.

• Grey-box models, as the name suggests, are usually a hybrid of the first two types above. They are commonly formed by taking a basic core motivated by known physics and then adding a black-box component with approximation properties suited to the problem of interest. A good example of a grey-box model is the Bouc–Wen model of hysteresis. In the Bouc–Wen model, a mass-spring-damper core is supplemented by an extra state-space equation which allows versatile approximation of a class of hysteresis loops [5,6].

In all of these cases, measured data from the system or structure of interest can be used in order to determine any unknown aspects of the model, e.g. any necessary undetermined parameters can be estimated. The use of measured data often means that uncertainty is introduced into the problem. There are two main sources of uncertainty caused by consideration of measured data. The first source is measurement noise; in general, other sources (noise) will contribute to measurements of the variable of interest and the direct distinction between signal and noise will be impossible. The second problem is encountered when a measured variable is itself a random process. In this case, only specific finite realizations of the process of interest can be measured; variability between realizations leads to variability between parameter estimates and thus gives rise to uncertainty.

In the past, the SI practitioner would generally implement the classical algorithms (i.e. least-squares minimization) as an exercise in linear algebra and would usually treat the resulting set of crisp parameter estimates as determining ‘the model’. Even if a covariance matrix were extracted, the user would usually use this only to provide confidence intervals or ‘error bars’ on the parameters; predictions would still be made using the crisp parameters produced by the algorithm. Such approaches do not fully accommodate the fact that a given set of measured data, subject to the sources of uncertainty discussed above, may be consistent with a number of different parametric models. It is now becoming clear—largely as a result of the pioneering work of James Beck and colleagues and more recently from guidance from the machine learning community—that a more robust approach to parameter estimation, and also model selection, can be formulated on the basis of Bayesian principles for probability and statistics. Among the potential advantages offered by a Bayesian formulation are the estimation procedure will return parameter distributions rather than parameters; predictions can be made by integrating over all parameters consistent with the data, weighted by their probabilities; evidence for a given model structure can be computed, leading to a principled means of model selection.

Adoption of Bayesian methods first became widespread in the context of the identification of black-box models; the methods have recently begun to occupy a central position within the machine learning community [7,8]. Bayesian methods for training multi-layer perceptron neural networks are a good example of this trend [9]; the Gaussian process model is also achieving wide popularity [10]. Most machine learning algorithms, like the neural networks and Gaussian processes already mentioned, are used to learn static relationships between variables; however, they can easily be used to model dynamical processes by assuming a NARX or NARMAX form for the mapping of interest [11,12]. A recent example of nonlinear system identification using Gaussian process NARX models can be found in [13]; this study is of interest because it shows how physical insight might be gained from the black-box GP NARX models. There has also been a body of work concerned with Bayesian parameter estimation for polynomial NARMAX models, a recent contribution can be found in [14].

In the context of white-box models, and in particular within the nonlinear SI community, the use of Bayesian methods has not been so widespread; however, their pedigree is as long. One can find references to Bayesian methods in a monograph on parameter estimation from 1974 [15], and dating from the same year is perhaps the first paper on Bayesian methods for structural dynamic SI [16]. To date, the most systematic and extensive development of Bayesian SI is the result of the work of James Beck and his various collaborators. Beck’s early work on statistical system identification is summarized in [17] and his transition to a Bayesian framework is given in [18]. This paper uses a Laplace approximation to remove the need to evaluate intractable high-dimensional integrals. Later, Beck & Au [19] introduce a Markov chain Monte Carlo (MCMC) method as a more general means of computing response quantities of interest represented by high-dimensional integrals. Bayesian methods of model selection are discussed in [20], and the paper also discusses the possibility of marginalizing over different model classes. A recent contribution [21] discusses identification and model selection for a type of hysteretic system model—the Masing model. Staying with hysteresis models, the paper [22] considers how MCMC can be used for Bayesian estimation of Bouc–Wen models and discusses a simple model selection statistic. Two recent developments which are of interest are the introduction of probability logic for Bayesian SI [23] and a method for potentially reducing computational expense for MCMC by selecting the most informative training data [24]. Bayesian methods for the system identification of differential equations have also been the subject of recent interest in the context of systems biology [25,26] and show considerable promise in the context of structural dynamics.

At this point, it is appropriate to define some notation. Here is used to represent a model structure. is then used to represent the vector of parameters within that model which requires estimation. Finally, is used to denote a set of observations which one has made about the system of interest, i.e. the measured data. As an example, one may consider the case study which is shown in §5, where one is attempting to create a white-box model of a dynamical system whose response is thought to be greatly influenced by friction effects. In this case, represents the hypothesized equation of motion of the system. Here θ represents the parameters within the equation of motion which require estimation—in the current example, this includes terms which modulate the level of viscous damping and friction in the system. The data, , consist of a time history of acceleration measurements which have been taken during a dynamic test. 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 θ as well as a model structure conditional on the data using

and respectively, where is a normalizing constant which ensures that integrates to unity. This is referred to here as the ‘marginal likelihood’ but can also be described as the ‘model evidence’ (because, as is shown in §2, it can provide evidence for candidate model structures). With equation (1.1), one converts an a priori probability density for the parameters θ into a posterior density having seen the data . If one desires a point estimate of the parameters, the usual course of action is to choose that which maximizes the posterior probability . Now, as the data is a constant of the identification problem, one is reduced to maximizing . It is often the case at this point, that an uninformative constant (and hence improper) prior is chosen, and this reduces the problem to that of maximizing , which is simply the likelihood of the data. The maximum a posteriori (MAP) estimate thus becomes maximum likelihood. If one were to further assume that the distribution of any measurement noise was Gaussian (with some extra conditions), the problem essentially becomes one of minimizing a least-squares cost/error function. The sequence of assumptions and approximations discussed above clearly loses much of the benefit of adopting a Bayesian approach in the first place. This paper will emphasize the benefits of a ‘full’ Bayesian methodology.

The layout of the paper is as follows. In §2, the fundamental principles behind a Bayesian approach to system identification are described and the benefits of using MCMC algorithms within a Bayesian framework are emphasized. Sections 3 and 4 are devoted to the description of various MCMC algorithms which can be used to address the issues of parameter estimation and model selection, respectively. These sections are not intended to be a thorough review but, instead, focus on those algorithms which have proved to be particularly useful and/or are based on unique concepts and methodologies.¹ Finally, in §5, a case study is used to demonstrate how MCMC can be used within a Bayesian framework to generate robust models of nonlinear dynamical systems.

2. Bayesian system identification

The problem of SI is easily stated: given measured data from a structure or system, how does one infer the equations of motion which ‘generated’ the data. This problem is not at all easy to solve; it is essentially an inverse problem of the second kind and can be extremely ill-posed even if the underlying equations are assumed to be linear in the parameters of interest [3]. Furthermore, the ‘solution’ may not even be unique. If the equations of motion are not linear in the parameters of interest, the difficulties multiply. Another issue is concerned with confidence in derived parameter estimates. This issue is a result of the fact that measurements or data from a system will, in reality, almost always be contaminated by random noise. Given a set of data of sampled system inputs x_i and outputs y_i, if there is no measurement noise, an identification algorithm should yield a deterministic estimate of the system parameters θ,

where the function id represents the application of the identification algorithm to the data . Now, if noise ϵ(t) is present on the input or output data (or both), θ will become a random vector conditioned on the data. In this context, one no longer wishes to find an estimate of θ, but rather to specify one’s belief in its value. If it is assumed that the noise is Gaussian with (unknown) standard deviation σ_ϵ, then the parameter σ_ϵ can be subsumed into θ, and inferred along with the model parameters. In probabilistic terms, instead of equation (2.1) one now has where represents the choice of model.

The usual objective of system identification is to provide a predictive model, i.e. one which can estimate or predict system outputs if a different system input is provided. In the probabilistic context described above, the best that one could do is to determine a predictive distribution. Suppose a new input sequence x* were applied to the system, one would wish to determine the density for the predicted outputs

noting all the dependencies.² The mean of this distribution would give the ‘best’ estimates for the predictions and the covariance would allow one to establish confidence intervals for them. However, one notes the presence of the parameter vector θ. In practice, one might use the θ value corresponding to the mean or the mode of the posterior parameter distribution; however, a truly Bayesian viewpoint on the prediction would require one to marginalize over the parameter estimates, i.e. to derive where —the posterior parameter distribution—is given by equation (1.1).

This is a very powerful idea: allowing for a fixed model structure, one is making predictions using an entire set of parameters consistent with the training data, with each point in the space of parameters weighted according to its probability given the data. In practice, there are considerable problems in implementing the full Bayesian approach, i.e. performing the intractable integral (2.4). One of the main advantages of using MCMC algorithms is that they allow one to generate samples from the posterior parameter distribution, even when the geometry of is complex and its probability density is very concentrated relative to the prior. These samples can then be used as part of Monte Carlo simulations, allowing one to propagate one’s parameter uncertainties without evaluating equation (2.4).

Another potential advantage of a Bayesian approach is that it may be possible to assess the relative evidence for a number of competing model structures. Suppose one believes that the true model structure is one of a finite number (the discussion here will closely follow [25]). In principle, one could imagine computing the probability of observing the data , given the particular model structure. If this quantity were available then, by defining a prior probability on each model structure, one could use equation (1.2) to select the model with the highest probability. Even more in the spirit of Bayesian inference, one could marginalize over all possible model structures weighted according to their probability; in terms of prediction, one would have

Furthermore, if one appeals to Bayes theorem in the form of equation (1.2) and assumes equal priors on the models, one arrives at a comparison ratio or Bayes factor which weights the evidence for two models in terms of marginal likelihoods of the data given the models.

The Bayesian approach to model selection is particularly attractive as the marginal likelihood rewards models for being high fidelity while also penalizing them for being overly complex. By automatically embodying Ockham’s Razor with regard to model selection, it follows that the adoption of a Bayesian approach can help to prevent overfitting. An intuitive explanation of this property is provided by MacKay [7], where it is suggested that a complex model will be capable of replicating a larger range of predictions than a simple model with relatively few parameters. As the probability density function must always be normalized it follows that, in a region where both models are able to replicate the same data, the marginal likelihood will be larger for the simpler model (figure 1). It is in this way that can be used to provide evidence for candidate model structures.

• Download figure
• Open in new tab
• Download PowerPoint

For a more detailed explanation, an information theoretic analysis of the marginal likelihood was originally discussed by Beck & Yuen [29] before being generalized in [21]. Noting that , it follows that the logarithm of the marginal likelihood can be written as

therefore The first term in the above equation is the posterior mean of the log-likelihood which is a measure of the average data fit of model . It follows that achieving a good fit to the training data will provide evidence for a candidate model structure. The second term in equation (2.10) represents the relative entropy between the prior and posterior. The marginal likelihood therefore penalizes models which are ‘complex’ where, a complex model is defined as that which is able to extract large amounts of information about the parameters θ from the data . It is important to note that the marginal likelihood is a function of the prior —it is possible to alter the evidence for by altering the prior distribution while maintaining the same model structure.

Unfortunately, the marginal likelihoods themselves (equation (1.3)) are often analytically intractable and numerically challenging because of their high-dimensional nature [26]. While one can resort to less informative model selection indicators which are simpler to compute (for example, as used in [22], a Bayesian generalization of the Akaike Information Criterion (AIC) [30] known as the Deviance Information Criterion (DIC)), it will be shown in §4 that there now exist MCMC methods which can be used to estimate the marginal likelihoods of different models/generate samples directly from .

As a final comment on the issue of model selection, it should be noted that there are already examples of the use of Bayesian strategies for model selection in the structural dynamics literature, the ‘Ockham factor’ defined in [19] being one of these. In [31], the authors use a Bayesian model screening approach in order to determine the appropriate nonlinear terms to include in a system model. The book [32] discusses Bayesian model selection in some detail.

3. Markov chain Monte Carlo for the posterior parameter distribution

The first set of algorithms reviewed here are those which are designed to generate samples from the posterior parameter distribution (equation (1.1)), while circumventing the need to evaluate the marginal likelihood (equation (1.3)). These methods involve the creation of an ergodic Markov chain—a Markov chain whose probability distribution tends to a functional form which is independent of time—which evolves through the parameter space (see [33] for a comprehensive discussion on the convergence of Markov chains). Simply stated, MCMC involves ‘forcing’ the stationary distribution of the Markov chain to be equal to (or at least proportional to) some target distribution such that, having allowed the chain to become stationary, it is effectively generating samples from the target. In the context of this section, the target distribution is .

Throughout the following text, one’s target distribution is written as π(θ)=π*(θ)/Z where π* is used to denote the unnormalized distribution and Z is the corresponding normalizing constant.

4. Markov chain Monte Carlo for the posterior model distribution

The algorithms described in §3 are all designed to generate samples from the posterior parameter distribution while circumventing the need to evaluate the marginal likelihood. While these methods are undoubtedly powerful, they do not allow one to evaluate the posterior model distribution and so can only be used to evaluate what is commonly referred to as the ‘first level of Bayesian inference’. Here, three algorithms are described which, using quite different methods, can be used to address both levels of inference—parameter estimation and model selection.

5. Case study

The case study shown here was originally conducted as part of a collaborative project with the University of Southampton (full findings are published in [62]); it is included here as it clearly demonstrates how using MCMC methods within a Bayesian framework can be used to quantify and propagate the uncertainties involved in modelling nonlinear dynamical systems.

Figure 4a shows a rotational energy harvester—a device which, via a ball–screw mechanism, is designed to convert low-frequency translational motion into high-frequency rotational motion (which can then be transformed into electrical energy). The device is mounted on a electro-hydraulic shaker while accelerometers are attached to the shaker and the oscillating mass (for a more detailed description of the experiment, see [62,63]). With the measured inputs and outputs (x and y) being the acceleration of the base and mass, respectively, the aim was to use a set of experimentally obtained data to infer a robust model of the device.

• Download figure
• Open in new tab
• Download PowerPoint

Referring to the schematic in figure 4b, the equation of motion of the energy harvester is

where z=y−x is the relative acceleration between the mass and base, l is the ball screw lead, c_m is mechanical damping, k is spring stiffness, m is the oscillating mass and J represents the system’s moment of inertia. The function is a friction model which is to be identified. In this case, two models, and , were considered. With the first, it was assumed that all of the parasitic losses in the device could be modelled using a linear viscous damper (which is equivalent to setting ) while, with the second, a hyperbolic tangent friction model was hypothesized. Assuming that M, k and l were already estimated with sufficient accuracy and employing a Gaussian likelihood with standard deviation σ_ϵ, the identification of models 1 and 2 involved estimating the parameter vectors {c_m,σ_ϵ} and {c_m,F_c,β,σ_ϵ}, respectively. Aside from obtaining probabilistic estimates for the parameters in each model, the aim here was to establish whether it is worth including the additional complexity of the hyperbolic tangent friction model.

Using TMCMC to generate 1000 samples from the posterior parameter distribution of each model, figure 5 shows the histograms of the resulting samples (where row 1 is model 1 and row 2 is model 2). Using these samples as part of Monte Carlo simulations, figure 6a,b shows the ability of model 1 and model 2 to replicate the training data (the filled grey regions in figure 6 represent 3σ confidence bounds). These results seem to indicate that the additional complexity of model 2 has allowed it to form a more accurate representation of the training data. Using TMCMC to analyse the marginal likelihood, the finding that confirms that model 2 is preferable. Finally, figure 6c shows the ability of model 2 to replicate a set of ‘unseen’ acceleration time history (data which were not used to train the model).

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

It is interesting to note that, although model 2 provides a better fit to the data, the confidence bounds on the predictions made by each model are of a similar magnitude. This is essentially due to a poor ‘initial phrasing’ of the problem. Specifically, when the likelihood was defined, it was assumed that the probability of witnessing each data point was a Gaussian distribution, centred on the prediction made by the model and with variance . The choice of a Gaussian distribution is justified somewhat by the Principle of Maximum Entropy [64] from which one finds that, having assumed the first 2 moments of the likelihood, the Gaussian distribution is that which minimizes the amount of additional information that must be assumed. However, having completed the analysis using such a likelihood, it can be observed that model 2 actually appears better able to replicate the experiment when at low amplitudes. From this, one may conclude that the probability of witnessing a data point, conditional on the model, actually varies with amplitude. As a continuation to this study one could propose a more complex likelihood before repeating the analysis. Potentially, one could then adopt a Bayesian approach to the selection of different likelihoods, thus preventing the selection of overly complex ‘error-prediction’ models (e.g. [65]).

6. Future work

Ultimately, with each sample generated using MCMC requiring a model run, the applicability of MCMC to Bayesian system identification problems is limited by computational cost. This places several restrictions on the types of problems which can be addressed. For situations where one’s model is expensive, a current stream of research is aimed towards the development of MCMC algorithms which are suitable for large-scale parallelization [66], and those which are able to reduce computational cost via the exploitation of interpolation methods (see [67] for example, where kriging is integrated into TMCMC). Further interest has been directed towards the scenario where one is confronted with large datasets from which to infer models. The work [24] proposes a method which allows the selection of small, highly informative subsets of data while, in [47,68], MCMC methods are proposed which allow the tracking of one’s parameter estimates as more data are analysed (helping to establish when a sufficient amount of data has been used).

7. Conclusion

In this paper, the authors have presented arguments for the adoption of a Bayesian framework for the system identification of nonlinear dynamical systems in the presence of uncertainty. Specifically, it has been highlighted how a Bayesian approach allows one to realize probabilistic parameter estimates in the presence of measurement noise, select high fidelity models which are not overfitted and make predictions which are marginalized over one’s parameter estimates and, in some cases, over a set of candidate model structures. It is then shown how many of the potential difficulties with such an approach can be overcome through the use of Markov chain Monte Carlo (MCMC) algorithms. A brief tutorial/review of six different MCMC algorithms is then given, each of which has been chosen because it has either proved to be particularly useful and/or is based on unique concepts and methodologies. The paper finishes with a case study, where an MCMC algorithm is demonstrated within a Bayesian framework to realize a model of a nonlinear, rotational energy harvester.

Data accessibility

The training data are available in the electronic supplementary material.

Authors' contributions

P.L.G. contributed the material on Markov chain Monte Carlo methods and carried out the analysis shown in §5 (Case study). K.W. contributed the material giving a general overview of Bayesian system identification and aided the article’s revision. Both authors contributed to the drafting of the manuscript and gave final approval for publication.

Competing interests

We declare we have no competing interests.

Funding

The authors would like to acknowledge the EPSRC Programme Grant ‘Engineering Nonlinearity’ EP/K003836/1 which funded the work in this paper as well as the collaborative project described in §5.