Bayesian parameter estimation and model selection for strongly nonlinear dynamical systems

Abstract

The Bayesian model selection and parameter estimation framework developed by Khalil et al. (J Sound Vib 332(15):3670–3691, 2013, Bayesian inference for complex and large-scale engineering systems. Ph.D. thesis, Carleton University, 2013) and Sandhu et al. (J Comput Methods Appl Mech Eng 282:161–183, 2014, Bayesian model selection and parameter estimation for a nonlinear fluid–structure interaction problem. Master’s thesis, Carleton University, 2012, 54th AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference, 2013) is extended to handle strongly nonlinear systems. The evidence required to estimate the posterior probability of each proposed model is computed using the Chib–Jeliazkov method. The posterior parameter samples generated using Metropolis–Hastings (M–H) Markov Chain Monte Carlo (MCMC) simulations are required by this method. The M–H MCMC-based parameter estimation procedure is complemented by an efficient particle filter and an ensemble Kalman filter for the strongly non-Gaussian state estimation problem. A strongly nonlinear system having multiple fixed (equilibrium) points is analyzed to demonstrate the efficacy of the algorithm.

Appendices

Appendix 1: Bayesian formulation of joint state and parameter estimation

The joint state and parameter pdf conditioned on data for a given model $\mathcal {M}_i$ can be expressed by the multiplication of the conditional state pdf and the parameter posterior. The following equations explicitly depend on a given model $\mathcal {M}_i$, but the symbol $\mathcal {M}_i$ is removed from these equations for brevity.

(28)

The joint pdf can be expressed as [10, 12, 13, 16]

$$\begin{aligned}&p(\mathbf {x}_1,\ldots ,\mathbf {x}_k | \varvec{d}_1,\ldots , \varvec{d}_J, \varvec{\phi })p(\varvec{\phi } |\varvec{d}_1,\ldots , \varvec{d}_J) \propto \nonumber \\&\quad p(\varvec{\phi }) \left[ \prod _{j'=1}^{d(1)-1} p(\mathbf {x}_{j'} | \mathbf {x}_{j'-1}, \varvec{\phi }) \right] p(\mathbf {x}_{d(1)} | \mathbf {x}_{d(1)-1}, \varvec{\phi })\nonumber \\&\quad \quad \quad \times \,p(\varvec{d}_1 | \mathbf {x}_{d(1)}, \varvec{\phi })\nonumber \\&\quad \quad \quad \quad \quad \vdots \nonumber \\&\quad \quad \left[ \prod _{j'=d(J-1)+1}^{d(J)-1} p(\mathbf {x}_{j'} | \mathbf {x}_{j'-1}, \varvec{\phi }) \right] p(\mathbf {x}_{d(J)} | \mathbf {x}_{d(J)-1}, \varvec{\phi })\nonumber \\&\quad \quad \quad \times \,p(\varvec{d}_J | \mathbf {x}_{d(J)}, \varvec{\phi })\nonumber \\&\quad \quad \quad \quad \quad \vdots \nonumber \\&\quad \quad \left[ \prod _{j'=d(J)+1}^k p(\mathbf {x}_{j'} | \mathbf {x}_{j'-1}, {\varvec{\phi }}) \right] \end{aligned}$$

(29)

Hence the parameter posterior can be expressed by [10, 12, 13, 16]

$$\begin{aligned} p(\varvec{\phi } |\varvec{d}_1,\ldots , \varvec{d}_J)&\propto p(\varvec{\phi })p(\varvec{d}_1,\ldots , \varvec{d}_J |\varvec{\phi })\nonumber \\&\propto p(\varvec{\phi }) \prod \limits _{j=1}^J \int \limits _{-\infty }^{\infty } p\left( \mathbf {x}_{d(j)} | \mathbf {x}_{d(j)-1}, {\varvec{\phi }}\right) \nonumber \\&\quad \times \,p\left( \varvec{d}_j | \mathbf {x}_{d(j)}, {\varvec{\phi }}\right) d \mathbf {x}_{d(j)}. \end{aligned}$$

(30)

Equations (28–30) highlight the relationship between the parameter posterior pdf and the conditional state pdf. For a non-Gaussian system, the accuracy of the model selection does not only depend on the nonlinear state estimator but also on a robust MCMC procedure to sample from the parameter posterior pdf.

Appendix 2: EKF

Although the details of the EKF are widely available (i.e., [29, 41–44]), a brief description is presented in this “Appendix” for completeness. EKF linearizes the nonlinear model and measurement operators in Eqs. (3) and (4). In the forecast step, the physics of the system is described through the model operator. The next state of the system is predicted using the current state and the proposed model. In the update stage, the predicted state is combined with an incoming observation. The Kalman gain matrix achieves the balance between the predicted state and the incoming observation. We choose that the measurement error is $\epsilon _j \sim \mathcal {N}(\varvec{0}, \varvec{\varGamma }_j)$ and the model error is $q_k \sim \mathcal {N}(\varvec{0}, \varvec{Q}_k)$. When an observation is available, the update step takes place as follows  (e.g., [43, 44])

$$\begin{aligned}&\varvec{C}_j = \left. \frac{\partial \varvec{h}_{d(j)}(\mathbf {x}_{d(j)}, \varvec{\varepsilon }_j)}{\partial \mathbf {x}_{d(j)}} \right| _{\mathbf {x}_{d(j)} = \mathbf {x}_{d(j)}^f, \varvec{\varepsilon }_{j} = \varvec{0}} \end{aligned}$$

(31)

$$\begin{aligned}&\varvec{D}_j = \left. \frac{\partial \varvec{h}_{d(j)}(\mathbf {x}_{d(j)}, \varvec{\varepsilon }_j)}{\partial \varvec{\varepsilon }_j} \right| _{\mathbf {x}_{d(j)} = \mathbf {x}_{d(j)}^f, \varvec{\varepsilon }_j = \varvec{0}} \end{aligned}$$

(32)

$$\begin{aligned}&\varvec{K}_{d(j)} = \varvec{P}_{d(j)}^f \varvec{C}_j^T \left[ \varvec{C}_j \varvec{P}_{d(j)}^f \varvec{C}_j^T + \varvec{D}_j \varvec{\varGamma }_j \varvec{D}_j^T \right] ^{-1} \end{aligned}$$

(33)

$$\begin{aligned}&\mathbf {x}_{d(j)}^a = \mathbf {x}_{d(j)}^f + \varvec{K}_{d(j)} ( \varvec{d}_j - \varvec{h}_{d(j)}( \mathbf {x}_{d(j)}^f , \varvec{0} ) ) \end{aligned}$$

(34)

$$\begin{aligned}&\varvec{P}_{d(j)}^a = (\varvec{I} - \varvec{K}_j \varvec{C}_j) \varvec{P}_{d(j)}^f. \end{aligned}$$

(35)

The forecast step is given by (e.g., [43, 44])

$$\begin{aligned}&\varvec{A}_k = \left. \frac{\partial \varvec{g}_k(\mathbf {x}_k, \varvec{f}_k, \varvec{q}_k)}{\partial \mathbf {x}_k} \right| _{\mathbf {x}_k = \mathbf {x}_k^a, \varvec{q}_k = \varvec{0}} \end{aligned}$$

(36)

$$\begin{aligned}&\varvec{B}_k = \left. \frac{\partial \varvec{g}_k(\mathbf {x}_k, \varvec{f}_k, \varvec{q}_k)}{\partial \varvec{q}_k} \right| _{\mathbf {x}_k = \mathbf {x}_k^a, \varvec{q}_k = \varvec{0}} \end{aligned}$$

(37)

$$\begin{aligned}&\mathbf {x}_{k+1}^f = \varvec{g}_k(\mathbf {x}_k^a, \varvec{f}_k, \varvec{0}) \end{aligned}$$

(38)

$$\begin{aligned}&\varvec{P}_{k+1}^f = \varvec{A}_k \varvec{P}_{k}^a \varvec{A}_k^T + \varvec{B}_k \varvec{Q}_k \varvec{B}_k^T. \end{aligned}$$

(39)

Appendix 3: EnKF

The ensemble Kalman filter addresses two major limitations of EKF. It avoids the approximate closure used in EKF and circumvents the computational cost of storing the error covariance matrix for high-dimensional state-space model [16]. The statistical closure stems from the fact that higher-order moments in the error covariance are discarded leading to an unbound variance growth in some cases. In the ensemble Kalman filter, the initial state ensemble matrix is generated by drawing samples from the prior pdf of the initial state vector. When a measurement is available, the update takes place as (e.g., [10, 16, 30])

$$\begin{aligned}&\varvec{d}_{j,i} = \varvec{h}_{d(j)}\left( \varvec{x}_{d(j),i}^f, \varvec{\epsilon }_{j,i}\right) \end{aligned}$$

(40)

$$\begin{aligned}&\bar{\varvec{d}}_j = \frac{1}{N} \sum \limits _{i=1}^N \varvec{d}_{j,i} \end{aligned}$$

(41)

$$\begin{aligned}&\bar{\mathbf {x}}_{d\left( j\right) }^f = \frac{1}{N} \sum \limits _{i=1}^N \mathbf {x}_{d\left( j\right) ,i}^f \end{aligned}$$

(42)

$$\begin{aligned}&\varvec{P}_{xd} = \frac{1}{N-1} \sum \limits _{i=1}^N \left( \mathbf {x}_{d\left( j\right) ,i}^f-\bar{\mathbf {x}}_{d\left( j\right) }^f\right) \left( \varvec{d}_{j,i}-\bar{\varvec{d}}_j\right) ^T \end{aligned}$$

(43)

$$\begin{aligned}&\varvec{P}_{dd} = \frac{1}{N-1} \sum \limits _{i=1}^N \left( \varvec{d}_{j,i}-\bar{\varvec{d}}_j\right) \left( \varvec{d}_{j,i}-\bar{\varvec{d}}_j\right) ^T \end{aligned}$$

(44)

$$\begin{aligned}&\varvec{K}_{d\left( j\right) } = \varvec{P}_{xd} \varvec{P}_{dd}^{-1} \end{aligned}$$

(45)

$$\begin{aligned}&\mathbf {x}_{d\left( j\right) ,i}^a = \mathbf {x}_{d\left( j\right) ,i}^f + \varvec{K}_{d\left( j\right) }\left( \varvec{d}_j - \varvec{d}_{j,i} \right) . \end{aligned}$$

(46)

Note $\mathbf {x}_{d(j),i}^a$ denotes the ith analysis or updated sample. Each updated sample is then propagated in time using the full nonlinear model in Eq. (3)

$$\begin{aligned} \mathbf {x}_{k+1,i}^f = \varvec{g}_k\left( \mathbf {x}_{k,i}^a, \varvec{f}_{k}, \varvec{q}_{k,i}\right) . \end{aligned}$$

(47)

The forecast mean of the state vector at time step $k+1$ is

$$\begin{aligned} \bar{x}^f_{k+1} = \frac{1}{N} \sum \limits _{i=1}^N \mathbf {x}_{k+1,i}^f . \end{aligned}$$

(48)

Appendix 4: PF-EnKF

The particle filters are fully non-Gaussian estimation algorithm based on the Bayesian framework (e.g., [26, 27]). While each particle in EnKF has equal weight of 1 / N, the particles in PF have varying weights which depend on the likelihood, the prior and the proposal density. While the prediction step is the same for the PF and EnKF, the difference lies in the update step whereby the weight of each particle is modified as (e.g., [26, 27])

$$\begin{aligned} w_{d(j),i} \propto w_{d(j)-1,i} \dfrac{p\left( \varvec{d}_j|\mathbf {x}_{d(j),i}\right) p(\mathbf {x}_{d(j),i} | \mathbf {x}_{d(j)-1,i})}{q\left( \mathbf {x}_{d(j),i} | \mathbf {x}_{d(j)-1,i}, \varvec{d}_j\right) } \end{aligned}$$

(49)

where $q(\mathbf {x}_{d(j),i} | \mathbf {x}_{d(j)-1,i}, \varvec{d}_j )$ is the proposal distribution, $p(\mathbf {x}_{d(j),i} | \mathbf {x}_{d(j)-1,i})$ is the prior and $p(\varvec{d}_j|\mathbf {x}_{d(j),i})$ is the likelihood function. A number of proposal densities are reported in the literature, for instance, refer to [10, 26, 27] for reviews. In this investigation, EnKF is used to construct the proposal density as it retains some non-Gaussian features in contrast to a Gaussian proposal obtained using EKF [10, 17, 18]. Using this proposal, the particles are nudged using the current observation in order to estimate the updated particle weights.

The update steps of the particle filter with EnKF proposal are as follows (e.g., [10, 17, 18])

$$\begin{aligned}&\varvec{d}_{j,i} = \varvec{h}_{d\left( j\right) }\left( \varvec{x}_{d\left( j\right) ,i}^f,\varvec{\epsilon }_{j,i} \right) \end{aligned}$$

(50)

$$\begin{aligned}&\bar{\varvec{d}}_j = \sum \limits _{i=1}^N w_{d\left( j\right) -1,i} \varvec{d}_{j,i} \end{aligned}$$

(51)

$$\begin{aligned}&\bar{\mathbf {x}}_{d\left( j\right) }^f = \sum \limits _{i=1}^N w_{d\left( j\right) -1,i} \mathbf {x}_{d\left( j\right) ,i}^f\end{aligned}$$

(52)

$$\begin{aligned}&\varvec{P}_{xd}= \dfrac{1}{1 - \sum \nolimits _{i=1}^N w_{d\left( j\right) -1,i}^2} \sum \limits _{i=1}^N w_{d\left( j\right) -1,i}\nonumber \\&\qquad \quad \times \,\left( \mathbf {x}_{d\left( j\right) ,i}^f-\bar{\mathbf {x}}_{d\left( j\right) }^f\right) \left( \varvec{d}_{j,i}-\bar{\varvec{d}}_j\right) ^T \end{aligned}$$

(53)

$$\begin{aligned}&\varvec{P}_{dd} = \dfrac{1}{1 - \sum \nolimits _{i=1}^N w_{d\left( j\right) -1,i}^2} \sum \limits _{i=1}^N w_{d\left( j\right) -1,i}\nonumber \\&\quad \qquad \quad \times \,\left( \varvec{d}_{j,i}-\bar{\varvec{d}}_j\right) \left( \varvec{d}_{j,i}-\bar{\varvec{d}}_j\right) ^T \end{aligned}$$

(54)

$$\begin{aligned}&\varvec{K}_{d\left( j\right) }= \varvec{P}_{xd} \varvec{P}_{dd}^{-1} \end{aligned}$$

(55)

$$\begin{aligned}&\mathbf {x}_{d\left( j\right) ,i}^a = \mathbf {x}_{d\left( j\right) ,i}^f + \varvec{K}_{d\left( j\right) }\left( \varvec{d}_j - \varvec{d}_{j,i} \right) \end{aligned}$$

(56)

$$\begin{aligned}&\tilde{w}_{d\left( j\right) ,i} = w_{d\left( j\right) -1,i} p\left( \varvec{d}_j | \varvec{x}_{d\left( j\right) ,i}^a\right) \nonumber \\&\quad \times \,\frac{\sum \nolimits _{l=1}^N w_{d\left( j\right) -1,l} \frac{1}{|\varvec{H}_f|}K\left( \varvec{H}_f^{-1} \left( \mathbf {x}_{d\left( j\right) ,l}^f-\mathbf {x}_{d\left( j\right) ,i}^a\right) \right) }{\sum \nolimits _{l=1}^N \frac{1}{N} \frac{1}{|\varvec{H}_a|}K\left( \varvec{H}_a^{-1} \left( \mathbf {x}_{d\left( j\right) ,l}^a-\mathbf {x}_{d\left( j\right) ,i}^a\right) \right) } \end{aligned}$$

(57)

$$\begin{aligned}&w_{d\left( j\right) ,i} = \frac{\tilde{w}_{d\left( j\right) ,i}}{\sum \nolimits _{i=1}^N \tilde{w}_{d\left( j\right) ,i}}. \end{aligned}$$

(58)

where $K(\cdot )$ is a Gaussian kernel and $\varvec{H}_a$ and $\varvec{H}_f$ are the bandwidth matrices for the proposal and the prior, respectively. These bandwidth matrices are defined as follows [10, 45]

$$\begin{aligned} \varvec{H}_a&= N^{\frac{-1}{n+4}} \left[ \dfrac{1}{N-1} \sum \limits _{i=1}^N \left( \mathbf {x}_{d(j),i}^a-\bar{\mathbf {x}}_{d(j)}^a\right) \right. \nonumber \\&\quad \times \,\left. \left( \mathbf {x}_{d(j),i}^a-\bar{\mathbf {x}}_{d(j)}^a\right) ^T \right] ^{\frac{1}{2}} \end{aligned}$$

(59)

$$\begin{aligned} {\varvec{H}}_{f}&=N_{eff}^{\frac{-1}{n+4}} \left[ \dfrac{1}{1 - \sum \limits _{i=1}^N w_{d(j)-1,i}^2} \sum \limits _{i=1}^N w_{d(j)-1,i} \right. \nonumber \\&\quad \times \,\left. \left( \mathbf {x}_{d(j),i}^f-\bar{\mathbf {x}}_{d(j)}^f\right) \left( \mathbf {x}_{d(j),i}^f-\bar{\mathbf {x}}_{d(j)}^f\right) ^T \right] ^{\frac{1}{2}} \end{aligned}$$

(60)

where N is the number of samples, n is the dimension of the state vector, and $N_{eff}$ is the effective ensemble size defined in Eq.(62). In the forecast step, the full nonlinear model in Eq. (3) is used to propagate each particle forward in time

$$\begin{aligned} \mathbf {x}_{k+1,i}^f = \varvec{g}_k\left( \mathbf {x}_{k,i}^a, \varvec{f}_k, \varvec{q}_{k,i}\right) . \end{aligned}$$

(61)

Degeneracy is a common problem with particle filters [26, 27], whereby only few particles with higher weights are clustered in the high-probability region. The particles in the low-probability region have negligible weights. A measure of this phenomenon is the effective ensemble size described as (e.g., [10, 26, 27, 46])

$$\begin{aligned} N_{eff} = \dfrac{1}{\sum \nolimits _{i=1}^N w_{d(j),i}^2}. \end{aligned}$$

(62)

One way to reduce the degeneracy is to increase the number of particles which in turn enhances the computational cost. Another option is to resample when the effective ensemble size is below a certain threshold. While several resampling algorithms are available [27, 47], the multinomial resampling [47, 48] has been successfully used in this work.