## Abstract

The problem of estimating parameters of nonlinear dynamical systems based on incomplete noisy measurements is considered within the framework of Bayesian filtering using Monte Carlo simulations. The measurement noise and unmodelled dynamics are represented through additive and/or multiplicative Gaussian white noise processes. Truncated Ito–Taylor expansions are used to discretize these equations leading to discrete maps containing a set of multiple stochastic integrals. These integrals, in general, constitute a set of non-Gaussian random variables. The system parameters to be determined are declared as additional state variables. The parameter identification problem is solved through a new sequential importance sampling filter. This involves Ito–Taylor expansions of nonlinear terms in the measurement equation and the development of an ideal proposal density function while accounting for the non-Gaussian terms appearing in the governing equations. Numerical illustrations on parameter identification of a few nonlinear oscillators and a geometrically nonlinear Euler–Bernoulli beam reveal a remarkably improved performance of the proposed methods over one of the best known algorithms, i.e. the unscented particle filter.

## 1. Introduction

Dynamic state estimation of dynamical systems is of central importance in the areas of structural vibration control and system parameter identification. The pioneering work of Kalman (1960) provides the exact solution to the problem of state estimation in linear dynamical systems with Gaussian additive noises. The filter is formulated in the time domain, is recursive in nature, implemented in real time, and the filter estimate minimizes the mean square error. The filter is extensively studied and has found widespread applications in the field of signal processing. Nevertheless, its inability to treat nonlinearity and/or non-Gaussian noises impairs its applicability to most problems of practical importance. In the context of mechanical systems, such problems are often characterized by their large dimensionality (possibly owing to a finite-dimensional projected semi-discrete form of an infinite-dimensional mathematical model) and diverse forms of nonlinearity (possibly resulting from deformation geometry or material inelasticity). Further complications may arise through the possible emergence of non-Gaussian terms in the temporal discretization of the semi-discrete form, which mostly appears as a system of stochastic differential equations. There are essentially two alternative strategies in the existing literature to treat at least some of these difficulties. The first strategy employs some form of linearization to approximate the given nonlinear state-space model and subsequently apply the Kalman filter to the linearized equations (e.g. Saito & Hoyshiya 1984; Imai *et al*. 1989; Ghanem & Shinozuka 1995; Brown & Hwang 1997; Ghosh *et al*. 2007). The second strategy, on the other hand, aims at numerically treating the posterior non-Gaussian probability density functions (pdfs) of the system states. Interestingly, recursive formulae for the evolution of such posterior pdfs conditioned on available measurements are derivable in an exact form. These are valid, in general, for nonlinear systems and Gaussian or non-Gaussian noises that could be additive or multiplicative in nature (Ristic *et al*. 2004). These formulae contain a large number of multidimensional integrals that are often computationally intractable. Consequently, the formulae remain essentially formal in nature. On the other hand, such recursive relations may form the basis of numerical strategies to tackle the problem approximately. This, indeed, has been the focus of a wide range of numerical techniques styled variously as particle filters, Monte Carlo filters, population Monte Carlo algorithms and sequential Monte Carlo sampling methods. In these methods, the posterior pdf is represented by a set of sample realizations of the states (obtainable using Monte Carlo simulations) along with a set of associated weights, and the algorithm ensures that these weights evolve recursively as the measurements are assimilated sequentially (Gordon *et al*. 1993; Tanizaki 1996; Merwe *et al*. 2000; Doucet *et al*. 2001; Liu 2001; Ristic *et al*. 2004; Cappé *et al*. 2007). A particle filtering technique of specific relevance to the present study is the sequential importance sampling (SIS) filter (Doucet *et al*. 2000). This filter is based on importance sampling from a proposal density. The choice of this density function itself constitutes an important step in the implementation of the method. For the special class of filtering problems involving nonlinear process equations and linear measurement models (with only additive Gaussian noises), it has been shown that the ideal importance sampling pdf (ispdf), which minimizes the variance of weights conditioned on previous states and measurements up to the current time, has a Gaussian form.

Parameter identification of nonlinear structural dynamical systems has attracted extensive research attention. A wide variety of identification strategies have been formulated in time and frequency domains. Worden & Tomlinson (2001) and Kerschen *et al*. (2006) provide extensive reviews of these developments. However, these reviews do not cover the methods based on dynamic state estimation techniques. Indeed, the application of particle filters to problems of relevance in structural dynamics is not yet widely researched. A few studies in the recent years have been done by Ching *et al*. (2006), Manohar & Roy (2006), Namdeo & Manohar (2007) and Sajeeb *et al*. (2007, submitted). The application of the Bayesian inference approach for model selection in the context of nonlinear dynamical system identification has been considered by Kerschen *et al*. (2003).

We presently consider the problem of state estimation and parameter identification in nonlinear dynamical systems when the measurements are noisy and nonlinear functions of the system states. The unmodelled dynamics in the formulation of process equation is represented as additive and/or multiplicative white noise vector processes (i.e. formal derivatives of Weiner vector processes). The process and measurement equations are discretized into maps via Ito–Taylor expansions. The discretization involves a set of multiple stochastic integrals (MSIs) that are, in general, mutually dependent, zero-mean, non-Gaussian random variables (local martingales). The system parameters to be identified are treated as auxiliary state variables, thereby converting the parameter estimation problem into one of state estimation. Measurements are allowed to be sufficiently smooth nonlinear functions of the system states. Thus, for instance, the reaction transferred to the support in a nonlinear oscillator is a nonlinear function of displacement and velocity processes; similarly, strains measured on a beam or a plate, undergoing large dynamic deformations, would be nonlinear functions of the discretized displacement and velocity fields. The problem at hand is then characterized by the presence of non-Gaussian random variables in the discretized equations, and a nonlinear functional relationship between the measured quantities and the state variables. These features pose difficulties in developing sequential Monte Carlo schemes for the filtering problem, and the present study aims at addressing these difficulties. Specifically, we propose a strategy to select an ideal importance sampling density function for this class of problems. This function is given by a weighted sum of multidimensional Gaussian density functions and is shown to be asymptotically exact as the nonlinearity in the measurement model tends to vanish. Furthermore, subject to nonlinear terms being sufficiently smooth, the formal accuracy of the expansion improves as the number of terms in the truncated Ito–Taylor expansion increases. By way of illustration, we consider a few examples from low-dimensional dynamical systems and the problem of state and parameter estimation of an Euler–Bernoulli beam undergoing large-amplitude (geometrically nonlinear) oscillations, with measurements being made on bending strains at a set of points on the beam. The proposed method is shown to perform significantly better than the unscented particle method (Merwe *et al*. 2000; Ristic *et al*. 2004) in terms of both the numerical accuracy and the computational time.

## 2. Governing equations

We consider dynamical systems governed by a set of stochastic differential equation (SDE) of the generic form(2.1)where is the state vector; is the deterministic forcing function; is the unknown vector of system parameters to be identified; is the drift vector or the state transition function; is the diffusion coefficient matrix; is a vector of standard Brownian motion processes; and is the vector of initial conditions (possibly random, but assumed to be independent of d*B*_{I}(*t*)). *B*_{I}(*t*) is taken to account for any unmodelled dynamics (while being questionable, this is standard in the particle filters literature). To facilitate determining the vector ** ϕ**, we declare the elements of

**as additional states evolving according to the equation(2.2)where is a diagonal matrix of the noise intensities; is a vector of increments of standard Brownian processes independent of**

*ϕ*

*B*_{I}(

*t*); and is modelled as random variables independent of

*u*_{0}, d

*B*_{I}(

*t*) and d

*B*_{II}(

*t*). The time dependence of

**, as in equation (2.2), is an artefact to convert the problem of parameter estimation to one of state estimation. Now, we define the augmented state vector and combine equations (2.1) and (2.2) to obtain(2.3)where , with**

*ϕ***=**

*s**n*+

*p*is the augmented state vector; is the augmented drift function; , with , is the augmented vector of standard Brownian motion processes; is the augmented diffusion coefficient matrix; and is the augmented initial condition vector (assumed to be independent of

**(**

*B**t*)). The implementation of the state estimation algorithm generally requires that the governing state equation (2.3) be discretized in

*t*. We presently accomplish this by using the Ito–Taylor expansion (Kloeden & Platen 1992). Thus, we consider discretization times with ;

*j*=0, 1, 2, …. In order to stay focused on the basic elements of the method, we presently adopt a uniform step size . A 1.5 local-order strong Ito–Taylor expansion of the multidimensional SDE (2.3) has its

*i*th component given by(2.4)where the operators are given by(2.5a)(2.5b)The quantities and appearing in equation (2.4) are known as MSIs. They are given by(2.6)(2.7)For a given

*k*, the MSIs constitute a set of zero-mean random variables that are mutually dependent and non-Gaussian in general. The lower-order MSIs,

*viz*.

*I*

_{j},

*I*

_{(j,0)}and

*I*

_{(0,j)}, may be shown to be jointly Gaussian while the higher-order MSIs, and , are non-Gaussian. We refer to Kloeden & Platen (1992, pp. 351–356) for the details of the evaluation of MSIs. Thus, the discretized version of equation (2.3) may be cast in the form(2.8)where ; ; ; ; and . The explicit dependence on time has been suppressed in equation (2.8) for notational conciseness. While

*w*

_{k}contain lower-order MSIs that are strictly Gaussian (with ),

**ψ**_{k}is the set of higher-order non-Gaussian MSIs. It may be shown that elements of

**ψ**_{k}are expressible as nonlinear functions of a set of Gaussian random variables, denoted by the vector

*α*_{k}. The joint distribution of the random variables

**w**_{k}and

*α*_{k}is determinable, and we illustrate the relevant details through specific examples later in this paper.

The identification of ** ϕ** is based on a set of measurements made on functions of the system states. Such measurements are invariably via the sampling of the system response at a set of discrete time instants. We assume that the sampling instants coincide with the discretization times used in discretizing the process equation (2.3). Accordingly, we postulate a model for the measurement given by(2.9)where ; ; ; and . Moreover,

**m**_{k},

*k*=1, 2, …, is a sequence of independent vector-valued Gaussian random variables with zero mean and given covariance. These random variables represent the effect of measurement noise and unmodelled mechanics of the problem in relating measurement

**y**_{k}to the system state vector

*x*_{k}. The functions

**h**_{k}(

**x**_{k}) and

*q*

_{k}(

**x**_{k}) are, in general, nonlinear.

Equations (2.8) and (2.9), taken together, constitute the governing equations for the estimation of the states. We now introduce the vector sequences *x*_{0:k}={**x**_{0}, **x**_{1}, …, **x**_{k}} and **y**_{1:k}={**y**_{1}, **y**_{2}, …, **y**_{k}}. The problem is then to determine the multidimensional posterior pdf *p*(**x**_{0:k}|**y**_{1:k}) and the filtering density *p*(**x**_{k}|**y**_{1:k}). Using *p*(**x**_{k}|**y**_{1:k}), one may readily determine the moments of interest (e.g. the mean and the covariance of **x**_{k} conditioned on **y**_{1:k}). A formal solution to this problem is available (see Ristic *et al*. 2004). It consists of the following pair of equations corresponding to the prediction and updating steps:(2.10)The implementation of these equations is, however, not feasible in most cases as they involve a large number of multidimensional integrals with non-Gaussian pdfs in the integrand. Consequently, one resorts to either approximate analytical schemes or numerical solutions based on Monte Carlo simulations.

## 3. State estimation using the sequential importance sampling filter

The SIS filter is a widely used method for dynamic state estimation via Monte Carlo simulations. We first present the salient features of this method and indicate its limitations, which the present study offers to overcome. Consider the problem of evaluation of the expectation(3.1)using Monte Carlo simulations. Here, *E*_{p}[.] denotes the expectation with respect to the pdf *p*(**x**_{0:k}|**y**_{1:k}). As there is no straightforward way of drawing samples from *p*(**x**_{0:k}|**y**_{1:k}), we rewrite equation (3.1) as(3.2)Here, *π*(**x**_{0:k}|**y**_{1:k}) is a valid pdf to be selected, so that (i.e. *π* is absolutely continuous with respect to *p*). While *p|π* is the Radon–Nikodym derivative, *π* is the ispdf or the proposal density function. It follows that:(3.3)and an estimate of this integral is given by(3.4)where **x**_{i},_{0:k} denotes the *i*th sample drawn from the ispdf *π*(**x**_{0:k}|**y**_{1:k}). Equation (3.4) can be recast as(3.5)In addition, by choosing the ispdf in the form(3.6)the weights *W*(**x**_{i,0:k}) may be shown to be recursively computable as(3.7)The implementation of the filter begins by assigning the value of 1/*N* to all the weights at *k*=0. However, given that weights are generated through a multiplicative scheme (3.7) and that they must be bounded, most weights, excepting one, go to zero as *k* increases. There is thus a gradual depletion in the number of samples with distinct values, and this compromises the ability to assimilate further measurements. Although using a very large *N* could address this issue in principle, it is impracticable. An alternative strategy involving a resampling step is implemented in which, at every *k*, an effective sample size, given by(3.8)is evaluated. Note that *N*_{eff}=*N* if all weights are equal and *N*_{eff}=1 if all weights excepting one are zero. Thus, when the effective sample at any time step falls below a threshold *N*_{thres}, a step involving resampling is implemented with a view to arresting the depletion of samples. The performance of the SIS filter crucially depends upon the choice of the ispdf *π*(**x**_{k}|**x**_{i,k−1}, **y**_{k}). Doucet *et al*. (2000) argue that the ispdf given by minimizes the variance of the weights conditioned on **x**_{0:k−1} and **y**_{1:k}. This optimal ispdf is, in general, not determinable. However, it can be shown that, when the process equation has additive Gaussian noise terms and the measured quantities are linear functions of the system states with the measurement noise being additive Gaussian, the optimal ispdf is also Gaussian. Thus, when the process and measurement equations are given, respectively, by(3.9)the ideal ispdf is Gaussian, with mean **M**_{k} and covariance matrix *Σ* given by(3.10)Owing to the non-availability of an analytical form of the ispdf in most cases, several approximate methods have been proposed (Doucet *et al*. 2000). However, each of them has specific disadvantages, e.g. a fixed pdf does not account for either process or measurement dynamics, the prior pdf does not account for the current observation, local linearization of the observation function relies on the accuracy of linearization, extended Kalman filter–particle filter and unscented particle filters (UPFs) are computationally expensive and depend on the accuracy of the extended and unscented Kalman filters, respectively.

In the present study, while the process equation (2.8) is allowed to have multiplicative Gaussian noise terms as well as non-Gaussian random noise terms, the measurement equation (2.9) may have terms that are nonlinear in the system states and multiplicative Gaussian noise terms. Clearly, these equations do not conform to the form given in equation (3.9) and, therefore, the problem of determining the optimal ispdf for this case requires the development of newer strategies.

## 4. An optimal ispdf based on the Ito–Taylor expansion

In order to implement the SIS filtering with the ideal ispdf for the problem at hand, we need to accomplish two tasks: (i) to determine the functional form of the ispdf and (ii) to develop a strategy to draw samples from this ispdf. In order to determine the ideal ispdf, we note that(4.1)For the process equation (2.8) and the measurement equation (2.9), the pdfs *p*(**x**_{k}, **y**_{k}|**x**_{i,k−1}) and *p*(**y**_{k}|**x**_{i,k−1}) are non-Gaussian and their functional form is not readily determinable. Consequently, drawing samples from the ideal ispdf would not be straightforward. We propose herein a strategy to overcome these difficulties in an approximate manner. We consider the nonlinear functions and in the measurement equation over and propose that these terms be Ito–Taylor expanded around the state . Accordingly, we obtain(4.2a)(4.2b)where are a vector and a matrix of MSIs, respectively; and *ρ*_{1} and *ρ*_{2} are the associated vector-valued and matrix-valued remainder terms, respectively. In order to simplify notations, we drop the time indices and, based on the above approximations, the measurement equation can be approximated as(4.3)where . The vectors and represent the resulting Gaussian and non-Gaussian terms as a result of the approximation, and .

From equations (2.8) and (4.3), it follows that the process equation and the approximated measurement equation have non-Gaussian noise terms, and we introduce the notation to collectively denote all these non-Gaussian random variables. It follows from equation (4.3) that is Gaussian in nature and is given by(4.4)where ; ; and is the covariance matrix of the resultant Gaussian noise vector in equation (4.3). Considering equations (2.8) and (4.3), it follows that **x**_{k} and **y**_{k} are jointly Gaussian given , i.e.(4.5)with , . If ** Θ** denotes the vector of Gaussian random variables consisting of the Gaussian vectors and denotes the corresponding coefficient matrix to be obtained from the joint distribution of

**x**_{k}and

**y**_{k}given , thenFrom the discrete map (2.8), it follows that(4.6)where ; ; and

*Q*is the covariance matrix of .

*Q*has the same functional form as

*C*

_{1}. Based on the theory of vector Gaussian random variables, it may be shown that(4.7)with and . Here, denotes the cross terms of the covariance matrix

*C*

_{3}.

As already noted, elements of the vector are non-Gaussian random variables that are described by nonlinear functions of Gaussian random variables. This means that, even though it is difficult to determine the pdf , drawing samples from this pdf is a simpler exercise. We use this in obtaining an approximation to the ispdf. Note that(4.8)Furthermore, we adopt a weighted Monte Carlo representation of the pdf , where are the sampled values and are the associated weights. It then follows from equation (4.8) that(4.9)This is the required optimal ispdf. This function, thus, turns out to be a weighted mixture of Gaussian pdfs. The strategy for drawing samples from such pdfs is well established (Liu 2001).

We also need to determine the terms and appearing in the equation for the weights (equation (3.7)). The determination of is straightforward since this density function may be shown to be given by , with (where is the covariance matrix of the measurement noise in equation (2.9)). On the other hand, owing to the presence of non-Gaussian terms in equation (2.8), determining is not so simple. To circumvent this difficulty, we follow the same strategy as has been used in deriving equation (4.9). Thus, we write(4.10)with being the weighted Monte Carlo representation for . The above pdf is again a weighted mixture of Gaussian densities, and it can easily be computed in order to derive the weights in a recursive manner via equation (3.7). A pseudo code of the proposed algorithm is given in appendix B in the electronic supplementary material.

In the case where the process and measurement equations conform to the format given in equation (3.9), the optimal non-Gaussian ispdf (equation (4.9)) may be shown to reduce to the ideal Gaussian ispdf valid for the system in equation (3.9). Details are available in appendix D in the electronic supplementary material.

Non-Gaussian noise terms in the present study essentially arise due to the application of Ito–Taylor expansion while discretizing the process equation or in the treatment of the nonlinear terms appearing in the measurement equation. The extent to which these terms appear in the formulation depends upon the nature of the process and measurement equations and, of course, the number of terms retained in the truncated Ito–Taylor expansion.

In the cases where the discretization of the process SDE or the treatment of the measurement equation does not lead to non-Gaussian MSIs, it is verified that the optimal ispdf in equation (4.9) reduces to a Gaussian pdf. In that case, drawing samples and calculating the corresponding weights become simpler. For further details, refer to appendix C in the electronic supplementary material.

It is well known that the Ito–Taylor expansion is a convergent series for sufficiently smooth drift and diffusion functions. Hence, for general nonlinear cases, the formal accuracy of the estimated states by the proposed method increases as the number of terms retained in the expansion increases (or as the time-step size decreases for a given number of terms in the truncated expansions).

In the context of the treatment of the non-Gaussian MSIs, it may be useful to recall the representation of continuous local martingales, vanishing at

*t*=0, as time-changed Brownian motions (Revuz & Yor 1991). In particular, considering any such non-Gaussian MSI (*α*being an appropriate multi-index), there is a time-changed Brownian representation for via the identity , where ( denotes the quadratic variation operator). The use of such representation should then enable one to use Gaussian characterizations of the MSIs, thereby simplifying the filter algorithm to a great extent. This point of view will however be explored elsewhere.

## 5. Numerical illustrations

We illustrate the formulations presented through three examples. First, we consider a nonlinear dynamical system governed by a scalar first-order SDE with parametric and external excitations. Second, we consider the problem of state and parameter estimation in an externally driven single degree-of-freedom Duffing oscillator, wherein measurements are made on the reaction transferred to the support. Finally, we consider state and parameter estimation in an Euler–Bernoulli beam undergoing geometrically nonlinear oscillations under external harmonic excitations and measurements are made on bending strains at a set of points on the beam surface. In all the examples, measurement time histories are synthetically generated using numerical simulations. Since the ‘true’ states are known in such cases, a metric of performance of the filter can be formulated. The implementation of the filter requires the sample size to be chosen and the estimates obtained are clearly subject to sampling fluctuations. To gain insights into the performance of the filter for a given sample size, we repeat the filtering process *N*_{f} times. We employ a measure of root mean square error (r.m.s.e.) given by(5.1)where is the true state of *j*th state in the *i*th performance of the experiment at and is the corresponding estimate obtained by filtering. The three examples taken together bring out the various features of the formulation, including the occurrence of non-Gaussian MSIs in process equations and the treatment of nonlinear measurement models. Limited comparisons with the results via the UPF are also made to highlight the conceptual and computational superiority of the proposed method.

### (a) Example 5.1: a nonlinear first-order system

We consider the process equation given by a first-order scalar SDE,(5.2)The measurement model is given by(5.3)with being mutually independent. The Ito–Taylor expansion for the process equation with an accuracy of order 1 (Milstein 1995) leads to the discrete map(5.4)From equations (5.3) and (5.4), we get(5.5)The vector of non-Gaussian random variables is given by(5.6)In order to treat the MSIs appearing in equation (5.6), we introduce a set of mutually independent Gaussian random variables with zero mean and unit standard deviation asFollowing Kloeden & Platen (1992), we obtain(5.7a)(5.7b)(5.7c)(5.7d)where denotes the number of terms retained in the series expansion of , given byProcess and the measurement equations are written in terms of ** θ** as(5.8)(5.9)with . Equations (5.8) and (5.9) are in the format as in equations (2.8) and (4.3), and, consequently, the formulation outlined in §4 becomes applicable to solve the state estimation problem. The measurement data are generated with

*K*=3 N m

^{−1},

*α*=9 N m

^{−3}, =15 N and

*λ*=2

*π*rad s

^{−1}. The process noise intensity is fixed at 1% of the r.m.s. value of the force,

*σ*=0.11 N. The simulation is carried out for 3 s with a step size of =0.01 s. The additive measurement noise level is taken to be 5% of the noise-free displacement amplitude (

*σ*

_{m}=0.056 m). Multiplicative noise intensities in the process and observation equations are taken as and , respectively. The sample size is fixed at 50. Number of samples from is taken to be 10. The threshold sample size below which resampling is needed is fixed at two-thirds of the original sample size. =3 is used while evaluating MSIs using equation (5.7

*a*)–(5.7

*d*).

Figure 1*a* compares the conditional expected value of with the reference trajectory along with the noisy observation . The estimated trajectory appears to closely follow the reference trajectory. The optimal ispdf in this case is non-Gaussian and is obtained as weighted-sum Gaussian pdfs. Figure 1*b* provides the details of these pdfs. The corresponding plots of histograms of the non-Gaussian random variables *θ*_{1} and *θ*_{2} are shown in figure 1*c*.

### (b) Example 5.2: a hardening Duffing oscillator

We consider a Duffing oscillator driven by harmonic excitations. The process equation in this case is(5.10)Here, (*c*, *K*, *α*) are the system parameters to be identified. We augment the above equations with additional state equationsWe have . The discretized equations based on a 1.5 local-order Ito–Taylor expansion are obtained as(5.11)The measurements are presently on the reaction transferred to the support and, accordingly, the measurement equation is obtainable as(5.12)The nonlinear terms appearing in equation (5.12) are expanded using Ito–Taylor's expansion and, retaining terms up to 1.5 local order, we get(5.13a)(5.13b)(5.13c)where . Note that the approximation of the measurement equation presently results in two non-Gaussian MSIs (*I*_{31} and *I*_{13}) in the second term and the filtering problem is handled within the framework outlined in §4. In the numerical work to synthetically generate the measurements, we take *m*=1 kg, *K*=17 N m^{−1}, *α*=120 N m^{−3}, *c*=2.47 N s m^{−1}, *n*_{f}=3, *p*_{i}=15 N (*i*=1, 2, 3), *λ*_{1}=*π* rad s^{−1}, *λ*_{2}=2*π* rad s^{−1} and *λ*_{3}=4*π* rad s^{−1}. The process noise intensity is fixed at 1% of the r.m.s. value of the force, *σ*=0.11 N. The simulation is carried out for 5 s at an interval of 0.01 s. The measurement noise level is taken to be 5% of the r.m.s. value of the true (noise-free) reaction (*σ*_{m}=4.21 N). The filtering problem is solved with 2000 samples and *N*_{thresh}=1330. The system parameters at *k*=0 are taken as distributed uniformly: ; ; . Figure 2*a* compares the estimate of the reaction transferred to the support along with the corresponding reference trajectory. Also shown in figure 2*a* is the contribution from only the linear terms to the reference trajectory. The influence of nonlinearity on the reaction is amply evident. The conditionally expected values of the parameters are shown, respectively, in figure 2*b*–*d*. The estimates are observed to approach their reference values as measurements are increasingly assimilated during the filtering. The evolution of the pdf of one of the system parameters, *viz*. *α*, is illustrated in figure 2*e*.

The state estimation problems being studied in this work may also be tackled using suboptimal filtering strategies currently available in the existing literature. Among these, the UPF is recognized to be one of the successful filtering methods. Thus, it is of interest to examine the performance of the present method vis-à-vis the existing UPF. For this purpose, we restrict ourselves to the problem of state estimation of Duffing's oscillator when all the system parameters are known. Accordingly, we derive the reduced process equation involving only *x*^{1} and *x*^{2}, with the measurement still being made on the reaction transferred to the support. For the purpose of illustration, we take *m*=1 kg, *K*=17 N m^{−1}, *α*=120 N m^{−3}, *c*=6.5 N s m^{−1}, *n*_{f}=1, *p*_{1}=*p*=15 N and *λ*=2*π* rad s^{−1}. The process noise intensity is at 1% of the r.m.s. value of the force, *σ*=0.12 N. The state estimation is carried out for 3 s using a step size of 0.01 s. To gain an insight into the influence of sampling fluctuations and robustness of state estimates as experiments are repeated, we consider the following two cases: (i) case I corresponds to repeated runs of the filter using different sets of particles with the same measurement time history and a fixed sample size, and (ii) case II corresponds to repeated runs of the filter, again with a fixed sample size but with different sets of particles and different measurement time histories. In each case, we consider two different measurement noise levels, namely 1% (sub-cases I(a) and II(a)) and 5% (sub-cases I(b) and II(b)) of noise-free reaction magnitude. While the Ito–Taylor expansion for the measurement equation is presently of (local) order 2.5, the process equation is expanded with 1.5 local order. Standard deviations of noise processes are fixed as *σ*_{m}=0.18 and 0.85 N. The sample sizes in these two cases are fixed as 90 (high measurement noise) and 300 (low measurement noise), respectively. In both the cases, the threshold sample size (below which resampling is needed) is two-thirds of the sample size. The number of repetitions of filtering is taken to be 20. In implementing the UPF, we use the procedure as outlined by Merwe *et al*. (2000) and take the scaling parameters =1, *β*=2 and *κ*=0. Results on the performance of the proposed filter and the UPF are summarized in table 1. It emerges that the proposed strategy is computationally more efficient (with the UPF consuming approx. five times more CPU time) and more accurate (as discerned from the values of r.m.s.e. for the two methods). From the sampling variance of the estimates of r.m.s.e., it is also observed that the present filter shows lesser sampling fluctuations than the UPF. Further studies are clearly needed to understand the nature of sampling fluctuations and convergence characteristics of the proposed method.

### (c) Example 5.3: state and parameter estimation in a nonlinear Euler–Bernoulli beam

We consider the dynamics of a simply supported, one-span, Euler–Bernoulli beam undergoing large-amplitude oscillations under the action of harmonic excitations. One of the simplified models for this system, incorporating nonlinear strain–displacement relations, is given by (Nayfeh & Mook 1979)(5.14)where is the independent spatial variable; is the transverse displacement field; *EI*(*x*) is the flexural rigidity; *EA*(*x*) is the axial rigidity; *m*(*x*) is the mass per unit length; *l* is the beam span; and *f*(*x*, *t*) is the externally applied load. Measurements are presently made on the strain component evaluated at a set of points on the beam surface. This strain component may be approximately related to the transverse displacement and axial displacement through the equation(5.15)To proceed further, we take the displacement field to be expressible by the series , with being the *j*th linear beam mode shape and being the *j*th generalized coordinate. Applying a Galerkin projection, the governing equations for the generalized coordinates (after including viscous damping) are(5.16)with and . We truncate the series for at the third term and account for the effect of unmodelled dynamics by a set of Gaussian white noise processes. Now, we define the augmented set of states as . Thus, Young's modulus *E* is declared as an additional state. This leads to the process equations(5.17)Assume that two (noisy) strain measurements are available at two different beam locations and that the channel noises are independent. The measurement equation is then given by(5.18)withwhere and are mutually independent Gaussian random variables with zero mean and unit standard deviation. The quadratic terms appearing in the measurement equation are further expanded in the Ito–Taylor series with 1.5 local order, and thus we obtain(5.19)Here, denotes the *j*th drift term in the governing SDE (5.17). Based on this, the measurement equation is obtained in the form of equation (4.3). It may be noted that all the MSIs that appear in the discretized process and approximate measurement equations are Gaussian in nature (see Saha & Roy (2007) for further details). This simplifies the filtering strategy considerably.

In the numerical work, we take the beam to have a rectangular cross section (0.025 m×0.003 m), *l*=0.75 m, mass density=7850 kg m^{−3}, 0.1 N s m^{−1} (*i*=1, 2, 3) and *E*=2.05×10^{11} N m^{−2}. The external load is given by , for ; *i*=1, 2, …, 5 with *f*_{0}=70 N m^{−1}, *λ*=10 rad s^{−1} and ={5.97, 1.45, 3.81, 3.05 and 5.60}^{T}. The process noise intensity is fixed at 1% of the r.m.s. value, with *σ*_{i}=0.1 N, for *i*=1, 2, 3. The simulation has been carried out for 1 s at an interval of 0.0001 s. The measurement noise level is 6% of the r.m.s. value of the maximum strain time history (*σ*_{m}=9.28×10^{−5}). The sample size of 500 with *N*_{thresh}=167 is used. As the initial guess, we assume *E* to be uniformly distributed, i.e. N m^{−2}. Axial strains on the bottom fibre of the beam at *x*=*l*/4 and *l*/2 are measured.

In the numerical work, it is verified that nonlinearity in the strain–displacement equation indeed makes significant contributions to the strain. Figure 3*a*,*b* shows the estimates and pdf of Young's modulus of the beam. While the initial guess of the pdf of *E* has a wide support ( N m^{−2}), the conditional expectation of *E* begins to converge to the reference value, *E*=2.05×10^{11} N m^{−2}, within a short time of assimilating the measured strains. Estimates of displacement and velocity at *x*=*l*/4 are observed to compare well with the corresponding reference trajectories in figure 3*c*.

## 6. Conclusions

The present work deals with state and parameter estimation in dynamical systems with nonlinear and non-Gaussian state-space models. The emphasis is on dynamical systems of interest in structural mechanics. Nonlinearity in these systems typically arises in relating strains with displacements and/or stresses with strains. Moreover, in problems of parameter identification, the parameters to be identified are declared as additional state variables and the governing equations are thus nonlinear even though the mechanical system is strictly linear for a given set of parameters. In addition, measurements made on quantities, such as forces transferred to the supports and strains, are often nonlinear functions of system displacements and velocities. In the present study, we have treated the noise terms appearing in the process and the measurement equations as white noise processes and integrations of the resulting stochastic differential equations are performed within Ito's theory. The study reveals that, even as noise terms in the governing SDEs are Gaussian in nature, discretizations of the SDEs through Ito–Taylor series lead to maps that contain both Gaussian and non-Gaussian sequences of random variables. The study addresses the complications in state and parameter estimation problems arising out of the coexistence of nonlinear and non-Gaussian terms in the system equations. Specifically, we develop a scheme to formulate the optimal ispdf applicable to this general class of problems, and this ispdf is embedded into an SIS filtering strategy. We show that the optimal ispdf is non-Gaussian and admits a representation as a weighted sum of Gaussian pdfs, and that we can effectively sample from this representation. While selecting the optimal ispdf, the proposed method incorporates knowledge of both the process dynamics and the current observation. A crucial aspect in the formulation is the use of Ito–Taylor expansions for the nonlinear terms in the measurement equation. Such expansions permit (at least formally) an improvement of the accuracy of representation by retaining higher-order terms. Limited comparisons of the performance of the proposed method with the UPF indicate that the present method offers substantial computational and numerical advantages.

## Acknowledgments

C.S.M. thanks the Aeronautical R&D Board, Government of India, and D.R. thanks the Naval Science and Technological Laboratory, Government of India, for the financial support provided to conduct this work.

## Footnotes

Electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2007.0075 or via http://www.journals.royalsoc.ac.uk.

- Received June 13, 2007.
- Accepted September 13, 2007.

- © 2007 The Royal Society