Fast Bayesian identification of a class of elastic weakly nonlinear systems using backbone curves

This paper introduces a method for the identi ﬁ cation of the parameters of nonlinear structures using a probabilistic Bayesian framework, employing a Markov chain Monte Carlo algorithm. This approach uses analytical models to describe the unforced, undamped dynamic responses of structures in the frequency – amplitude domain, known as the backbone curves . The analytical models describing these backbone curves are then ﬁ tted to measured responses, found using the resonant-decay method . To investigate the proposed identi ﬁ cation method, a nonlinear two-degree-of-freedom example structure is simulated numerically and analytical expressions describing the backbone curves are found. These expressions are then used, in conjunction with the backbone curve data found through simulated experiment, to estimate the system parameters. It is shown that the use of these computationally-cheap analytical expressions allows for an extremely ef ﬁ cient method for modelling the dynamic behaviour, providing an identi ﬁ cation procedure that is both fast and accurate. Furthermore, for the example structure, it is shown that the estimated parameters may be used to accurately predict the existence of dynamic behaviours that are well-away from the backbone curve data provided; speci ﬁ cally the existence of an isola is predicted. & 2015 The Authors


Introduction
As nonlinear behaviour in structural dynamics becomes increasingly important, so does the need to accurately describe its characteristics.Whilst a number of reliable and efficient approaches have been developed for the identification of the linear characteristics of nonlinear systems, for example [1,2], the identification of nonlinearities still poses a number of challenges.Existing approaches to the problem in multi-degree-of-freedom systems include the restoring force method [3], NARMAX methods [4] and the conditioned reverse path method [5].An extensive review of these, and other approaches, is given in [6].Methods that are of specific interest to this work include the nonlinear resonant decay method [7,8], and the use of Bayesian approaches [9][10][11][12].
Bayesian methods for parameter identification problems seek to determine the probability that a particular set of parameter values are "correct", given some experimental data.The resulting probability density function, referred to as the posterior, is difficult to find analytically for nonlinear systems.Alternatively a Markov chain Monte Carlo (MCMC) method may be used, alongside simulations of the system, to sample from the posterior.Due to the complexity of modelling nonlinear systems, these simulations often require numerical integration through time at great computational expense.This issue can be mediated, to an extent, through the use of surrogate models and algorithms which are suitable for parallelisation [13,14], or through using carefully selected subsets of the experimental data [12,15].However, computational expense still proves limiting for the application of these Bayesian methods to large nonlinear systems with many unknown parameters.
In this paper, a method is introduced that allows Bayesian system identification to be conducted in a computationally efficient manner without the need for model surrogates, or the need to use fewer experimental data.The method proposed here relies on inferring parameter estimates from experimentally measured backbone curves (i.e. the unforced, undamped response of nonlinear systems) and analytical expressions describing the curves.This approach results in significantly lower computational expense than those using time-domain numerical models.
Two nonlinear oscillators, each with two degrees-of-freedom, are considered in this work.These are introduced in Section 2, along with a procedure for the acquisition of data describing the backbone curves.This procedure is similar to the resonant decay method, as detailed in [8].A number of techniques may be used for developing analytical models for the backbone curves; for example, harmonic balancing, multiple scales or the first-order normal form technique [16][17][18].Here, however, we use the second-order normal form technique [19] as it is not only able to analytically describe modal interactions in backbone curves [20,21], but it is also suited for automation due to its matrix-based formulationan essential property for the identification of larger systems.A brief outline of the application of the second-order normal form technique to the example system is given in Section 3, and further details are provided in Appendix A.
Section 4 is dedicated to the discussion of the Bayesian parameter estimation approach taken here, with specific details regarding the formulation adopted to accommodate the analytical model.This discussion also examines the use of the residuals of the expressions describing the backbone curves, thus eliminating the need for explicit solutions which can be computationally demanding to calculate.The results of the parameter identification for the example systems are presented in Section 5 where it is shown that the predictions of the responses of the system, made using the estimated parameters, are accurate in response regions well-away from those described by the data.Specifically, backbone curves that are not described by the measured data can be obtained, using the estimated parametersdemonstrated here by the prediction of an isola.Finally, conclusions are drawn in Section 6 and discussion is given to the potential future applications of the approach outlined in this work.

The example systems
Fig. 1 shows a two-degree-of-freedom oscillator with a symmetric, linear structure and three nonlinear springs.The underlying linear structure is composed of two identical springs and viscous dampers that ground the masses, as well as a spring and a viscous damper connecting the masses.The grounding springs and connecting spring have stiffness constants k 1 and k 2 respectively, and the grounding dampers and connecting damper have damping constants c 1 and c 2 respectively.Additionally, two nonlinear cubic springs, with constants α 1 and α 3 , ground the masses and a nonlinear cubic spring, with constant α 2 , connects the masses.The two masses have displacements x 1 and x 2 , and both are of mass m.Sinusoidal forcing at frequency Ω f is applied to the masses at amplitudes P 1 and P 2 as shown.
The equation of motion for this system may be written as where M, C and K are f2 Â 2g mass, damping and stiffness matrices, respectively.The f2 Â 1g vectors x, Γ x and P x describe the physical displacements, nonlinear terms and the forcing amplitudes respectively, and are written as In this paper, two different parameter sets for the oscillator shown in Fig. 1 are considered.These parameter sets, labelled 1NL and 3NL, are given in Table 1.Parameter set 1NL describes a system that has previously been investigated in [22], and has only one nonlinear spring, i.e. α 2 ¼ α 3 ¼ 0. Parameter set 3NL describes a system with three nonlinear springs, i.e. all nonlinear parameters are non-zero.Aside from the nonlinear coefficient, the only parameter that differs between the two systems is the coefficient of the linear spring connecting the masses, k 2 , resulting in different linear natural frequencies between 1NL and 3NL.As the parameters describing the external forcing -P 1 , P 2 and Ω fare not characteristics of the physical system, they are not listed in Table 1.Furthermore, for the technique presented here, it is not necessary for the external forces to be measured.Instead, the forcing parameters are defined whenever a forced response is considered.

An experimental method for measuring backbone curves
The data used for the identification process in this paper have been collected through simulated experiment.This employs a numerical integration routine, specifically MATLAB's integration algorithm implementing Runge-Kutta (4,5) [23], to simulate the procedure detailed below.Additionally, Gaussian white noise is added to the data before processing, to simulate noise generated by sensors.
This paper concerns the identification of nonlinear characteristics, and is based on the assumption that the underlying linear system is either known or can be found, for example using the methods proposed in [1].Here, as the linear parameters of the system are known, see Table 1, this step is omitted.
Fig. 2 demonstrates the procedure that is used to find the backbone curves of system 1NL.The panels ða 1 Þ and ða 2 Þ show the experimental step of this procedure for the first and second backbone curves, labelled S1 and S2 respectively.Here, the system is forced at a constant amplitude and the forcing frequency, Ω f , is varied such that the response may approach the backbone curves.As measurements of the forced responses are not required, the frequency, Ω f , may be varied continuously in time.This results in lower experimental cost in comparison to techniques that require the system to reach steady-state between discrete frequency steps.Further details of how this may be achieved can be found in [24].Once the response of the system is close to a backbone curve, the forcing is released, i.e. the forcing amplitude is set to zero, such that the transient response of the system then decays along, or close to, the backbone curvesee [24] for further details.
In panels ða 1 Þ and ða 2 Þ, a thin-black line and dotted-green line show the forced response branches and backbone curves respectively.The paths taken along the forced response branches to approach the backbone curves are shown by a thickblack line and the points of release are represented by blue crosses.The points on the backbone curve along which the system response then decays are shown by large green dots (although here these points are illustrative and would not be realised in the frequency-amplitude projection until after the data-processing step).
The next step involves transforming the time-domain data describing the displacements of the physical coordinates, x, into the displacements of the linear modal coordinates, q.This step is necessary as the second-order normal form technique, used to generate an analytical model of the system, uses a linear modal transform to decouple the linear terms in the equation of motion, see Section 3 for further details.The linear modal displacements are found using the inverse linear modal transform, written as q ¼ Φ À1 x where Φ is a modeshape matrix whose nth column describes the modeshape of the nth linear mode.Due to the symmetry of the underlying linear structure of the systems considered here, the modeshape matrix may be written as Panels ðb 1 Þ and ðb 2 Þ in Fig. 2 show the transformed time-domain decay data for the backbone curves S1 and S2 respectively.Both panels are in the projection of time, t, against the modal displacements q 1 and q 2 , where q 1 ðtÞ and q 2 ðtÞ are represented by solid-blue and dotted-red lines respectively.Embedded plots in these panels show a portion of the decay data in detail, and illustrate that for the decays along both backbone curves, the linear modes respond at the same frequency.Furthermore it can be seen in panel ðb 1 Þ, showing S1, that linear modal coordinates are in anti-phase, and panel ðb 2 Þ shows that q 1 and q 2 are in-phase along S2. 1 This frequency and phase information is used later for the application of the second-order normal form technique.The final step in the experimental acquisition of backbone curve involves the conversion of the time-displacement data into the frequency-amplitude domain.This is achieved using a moving-window discrete Fourier transform, with a window size of 5 periods, [25].For the approach taken here the harmonic responses of this system are neglected and only the fundamental responses are considered.The fundamental responses of q 1 and q 2 are denoted u 1 and u 2 respectively, and are assumed to be sinusoidal with frequency Ω and amplitudes U 1 and U 2 respectivelywhere Ω, the response frequency, is distinct from Ω f , the forcing frequency.As the estimated frequency is likely to be more accurate for a signal of larger amplitude, and as it is known that both linear modes share response frequencies, the dominant mode is used to obtain the response frequency for both modes.It can be seen from panels ðb 1 Þ and ðb 2 Þ in Fig. 2 that q 1 4 q 2 on S1 and q 2 4q 1 on S2.
Panels ðc 1 Þ and ðc 2 Þ in Fig. 2 show the frequency-amplitude data for backbone curves S1 and S2 respectively, and represent the data used for the identification of the nonlinear parameters in system 1NL.These are both shown in the projection of the common response frequency, Ω, against the fundamental response amplitudes, U 1 and U 2 .Due to the noise applied to the time-domain signal, the translation into the frequency domain is inaccurate for low-amplitude signals.Therefore, the decay data along each of these backbone curves have been truncated as the low-amplitude mode approaches zero.As can be seen in panels ðc 1 Þ and ðc 2 Þ, this results in the truncation of the dominant modes at relatively high amplitudes.The data are also truncated at the start of the decay in order to remove the transient effects as the response Fig. 2. A diagram of the procedure used to find the backbone curves of system 1NL through simulated experiment.The left-hand panelslabelled ða 1 Þ, ðb 1 Þ and ðc 1 Þand right-hand panelslabelled ða 2 Þ, ðb 2 Þ and ðc 2 Þdepict the process of data acquisition for backbone curves S1 and S2 respectively.converges towards the backbone curves.This transient behaviour is reduced if the system is released close to the backbone curve, allowing more high-amplitude data to be used for identification.
The systems considered here are discrete, however the same technique may be applied to continuous systems, provided that the linear modal displacements may be estimated.Additionally, as the technique presented here is concerned with identifying nonlinear parameters, only backbone curves that exhibit nonlinear behaviour, and their constituent linear modes, require experimental measurement [26].This approach may be used to reduce the experimental and computational costs.

Calculating the backbone curves
To find the backbone curves of the system described by Eq. ( 1) the associated unforced, undamped dynamics are considered, which may be written as The second-order normal form technique is now used to transform the equation of motion, Eq. ( 4), into a set of timeinvariant equations describing the backbone curves.An outline of the approach is given here, and a more detailed description is given in Appendix A. This approach is similar to that used in [20,21], and for a complete description of the second-order normal form technique, see [19].This technique is limited to weakly nonlinear systems, i.e. the nonlinear terms are small relative to the linear terms.The error resulting from this assumption may be reduced by computing the technique to a higher order of accuracy [27].
As described in Section 2.2, the first step in the second-order normal form technique is the linear modal transform x ¼ Φq where Φ is a modeshape matrix which, for the systems considered here, is given in Eq. ( 3).Applying this transform to Eq. ( 4) allows us to write where N q is a vector of nonlinear terms and Λ is a diagonal matrix whose nth diagonal term is the square of the nth linear natural frequency, ω nn 2 .For the systems considered here, one may write and for both systems one may write The next step is to apply the nonlinear near-identity transform q ¼ u þ h where u and h describe the fundamental and harmonic contents of q respectively.Here, the harmonics are neglected as only the fundamental responses of the decay along the backbone curve are measured.However in Appendix A the harmonics are retained at this stage.By removing the harmonics, one may write q ¼ u, which may be substituted into Eq.( 7), giving N q ðqÞ ¼ N q ðuÞ.As it is assumed that the fundamental response of the nth linear mode, u n , is sinusoidal, one may write where U n and ϕ n are the amplitude and phase of u n respectively.As, in the cases considered here, the linear modes respond at the same frequency (as shown in Fig. 2), both modes share the fundamental response frequency Ω.Note that the subscripts p and m correspond to the positive and negative (plus and minus) signs of the complex exponents respectively.The nonlinear near-identity transform results in the resonant equation of motion, written as where N u is a vector of resonant nonlinear termspopulated with the resonant terms from N q ðuÞ.The process used to determine which terms in N q ðuÞ are resonant is detailed in Appendix A, where it is found that N u ðuÞ may be written as where α p ¼ ðα 1 þ α 3 Þ=2m and α m ¼ ðα 1 À α 3 Þ=2m.Substituting Eq. ( 10) into Eq.( 9) forms the resonant equation of motion, in which all terms resonate at frequency Ω.As detailed in Appendix A, this can then be used to find the set of time-invariant equations describing the backbone curves which, for these systems, are given by i.e. p ¼ þ1 for backbone curves where the modes are in-phase, and p ¼ À1 for backbone curves where the modes are in anti-phase.Eqs. ( 11) and ( 12) are used in conjunction with the resonant decay data to identify the nonlinear parameters.The approach used for this has a Bayesian framework, as discussed in the following section.
Here, we find the parameters describing the equation of motion in the physical coordinates, Eq. ( 4); however, if one requires a modal model then the parameters describing the modal equation of motion, Eq. ( 5), may be identified.

Bayesian system identification
This section provides a brief overview of the Bayesian approach employed here, which uses the analytical expressions describing the backbone curves, Eq. ( 11), to form a physical-law based model M. The model contains a vector of parameters θAΘ & R N θ which require estimation (such that, in this case, θ ¼ fα 1 ; α 2 ; α 3 g).These estimates are probabilistic such that future predictions made using M are robust against parameter uncertainties.Using the experimental data set, D, consisting of the backbone curve data, the parameters, θ, are inferred.The probability of the parameters, θ, conditional on the model, M, and the experimental data set, D, can be obtained using Bayes' theorem: where PðθjD; MÞ is the posterior distribution and PðθjMÞthe prioris a probability distribution describing one's knowledge of θ before the data were known.The denominator of Eq. ( 13) can be interpreted as a normalising constant which ensures that the posterior integrates to unity, and is therefore defined as The term PðDjθ; MÞ is referred to as the likelihood and describes the probability of witnessing the data set D given that one believes that the behaviour of the real system can be replicated using model M with parameters θ.Evaluating Eq. ( 13) therefore requires the definition of a prediction-error model whose parameters can be included in θ.
In this case, the data set consists of K measurements of both U 1 and U 2 , such that D ¼ fU ð1Þ 1 ; …; U ðKÞ 1 ; U ð1Þ 2 ; …; U ðKÞ 2 g.Assuming that the probabilities of witnessing separate data are always mutually independent and employing a Gaussian prediction-error model, the likelihood is where Û 1 and Û 2 represent predictions made by the model and σ 1 and σ 2 are parameters describing the standard deviations of the likelihood; which, in this case, are also included in θ.In its current form, evaluation of the likelihood requires one to assemble an explicit solution for Û 1 and Û 2 in terms of θthis can be algebraically demanding and may prevent the method proposed here from being applied to more complex systems.As an alternative, one can use the residuals of Eqs.(11), written as ϵ 1 and ϵ 2 and defined as where α 1 , α 2 and α 3 are taken from elements of θ, and the data set is now defined as 1 ; …; U ðKÞ 1 ; U ð1Þ 2 ; …; U ðKÞ 2 ; Ω ð1Þ ; …; Ω ðKÞ g.Deriving the probability of witnessing the residuals, given θ, is a complex problem.In this case the authors have chosen to use the Principle of Maximum Entropy [28] which states that, having defined its first two moments, choosing a likelihood of the form assumes the least amount of additional information.The consequences of selecting such a prediction-error model are discussed further in Section 5.2.
This reduces the complexity of the problem, as the residual terms, ϵ 1 and ϵ 2 , can be found directly from the substitution of the data and parameters into Eqs.( 16), without the need for algebraic manipulation.It is this approach that is adopted here, i.e.Eq. ( 17) is used to evaluate the likelihood.Furthermore we define the prior as having a uniform distribution.With the likelihood and prior defined, the next task is to generate samples from the posterior distribution.In the current paper this is achieved through the use of Markov chain Monte Carlo (MCMC) methods2 algorithms which involve the evolution of an ergodic Markov chain whose stationary distribution is proportional to the posterior distribution (this is a particularly useful feature as it allows one to generate samples from PðθjD; MÞ whilst circumventing the need to evaluate the N θ dimensional integral given by Eq. ( 14)).In this case, the authors employed a variant of the well-known Metropolis algorithm [29] detailed in [30].
It is important to note that, as each sample generated using MCMC requires a model run, the method proposed here is very computationally efficient, as each model run simply requires the evaluation of Eqs.(11).However, if one were attempting to infer parameter estimates using time history data, for example, each model run would require the numerical time-domain integration of the equations of motion of the system.The speed of the method proposed here is such that all of the MCMC simulations shown subsequently were conducted in a matter of secondsthis is without any parallel processing or the use of model emulators.

Parameter identification for the example systems
Here, a uniform prior distribution is utilised.Table 2 shows that the prior limits have been set such that the nonlinear parameters are positive, as the backbone curves depicted by the experimental data show a hardening nonlinear behaviour in both systems.

Table 2
The parameters limits used for the identification of systems 1NL and 3NL.

Parameter
Lower limit Upper limit 5.1.Parameter identification for system 1NL Fig. 3 shows histograms of the results of the identification procedure for system 1NL.Additionally, Table 3 shows, for each of the parameters, the true and mean estimated values along with the relative error of the mean estimated values.It can be seen from Fig. 3 that the variances of the estimated values of the nonlinear parameters are low.This indicates that the effects of these parameters are well-represented by the backbone curve data, and hence there is little uncertainty regarding their value.Fig. 3 shows that the true value of α 1 (which is 0.5) is outside of the limits of the estimated values, whereas one would expect the true value to lie within these limits.The reason for this is that the second-order normal form technique, used to describe the backbone curves, is an approximate method.Therefore, the parameter estimates describe an optimisation of the fit between the true system and the approximate model, leading to the discrepancy between the true and estimated values.The relative error of α 1 is low, at less than 2 percent, suggesting a successful identification procedure.Additionally, the estimated values of α 2 and α 3 are very low suggesting that their identification was also accurate.For cases where there is a greater discrepancy between the true and the approximate model, the second-order normal form technique may be computed to a higher order of accuracy [27].
In order to investigate the significance of the error in the estimation of the nonlinear parameters, the backbone curves were calculated for both the true and estimated parameters using the numerical continuation software AUTO-07p [31].These are shown in Fig. 4, where solid-green and dashed-red lines represent the backbone curves described by the true and estimated parameters respectively.It can be seen that the two sets of backbone curves are very similar and, importantly, it can also be seen that they are close at amplitude and response frequencies well-beyond those described by the data, which are represented by blue dots.

Uncertainty propagation for system 1NL
An advantage of the Bayesian analysis adopted here is that, using the samples generated by MCMC, it is possible to propagate the uncertainties in one's parameter estimates into future predictions.This involves conducting a Monte Carlo simulation, where each model run uses parameters from the vector θ, which have been sampled from the posterior parameter distribution (using MCMC).Recalling that θ also includes parameters which were used to define the likelihood (σ 1 and σ 2 in this case), each model run is also "corrupted" with noise generated according to the prediction-error model Fig. 4. The backbone curves S1 and S2 for system 1NL.These have been calculated using AUTO-07p (a numerical continuation package) for both the true and estimated parameters, represented by solid-green and dashed-red lines respectively.Blue dots show the data used for the identification of these parameters.These results are shown in the projection of the common response frequency, Ω, against the amplitude of displacement of the second mass, X 2 .(For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

Table 3
The true and mean estimated parameter values, along with the relative error of the estimated values, for system 1NL.

Parameter
True value Mean estimated value Relative error with parameters σ 1 and σ 2 .Analysing the statistics of this ensemble of model predictions, Fig. 5 shows the 73σ confidence bounds that arose when predicting the backbone curves of system 1NL.It can be seen that, for the most part, the data are within the confidence bounds.Any discrepancies (where data are outside the confidence bounds) may be due to erroneous assumptions about the prediction-error model.Specifically, for future work, the authors aim to investigate the assumption that the standard deviations of the prediction-error model are independent of the amplitude response of the system.This could involve analysing a set of different prediction-error models, thus allowing the probability of these competing models to assessed (also using a Bayesian framework) [32].

Parameter identification for system 3NL
The data used for the identification of system 3NL were also found through simulated experiment, in a similar procedure to that applied to system 1NL, as described in Section 2. The results of this are shown in Fig. 6, where the blue and green histograms show the identified values for the nonlinear parameters and standard deviations of the likelihood respectively.Table 4 also gives the true, mean estimated, and relative error for each of the parameters for system 3NL.As with system 1NL, the mean estimated values are very close to the true values and the low variance suggests a high confidence (i.e. the effects of the parameters are well-represented by the data).Also, the approximate nature of the model again results in true parameter values that lie outside of the limits of the estimated values.
Table 4 shows that the maximum relative error in the estimation of the nonlinear parameters is less than 2.5 percent.To test the influence of these errors, Fig. 7 shows the backbone curves for system 3NL, calculated using the numerical continuation software AUTO-07p, for both the true and estimated parameters.These are represented by solid-green and Fig. 5. Uncertainty propagation for the backbone curve decay data for system 1NL.This shows the backbone curve decay data in the projection of the fundamental response frequency, Ω, against the fundamental response amplitudes of the linear modal coordinates, U 1 and U 2 .These are represented by blue and red dots respectively.The black lines show the confidence bounds at 7 3σ.(For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)Fig. 6.Histograms, for system 3NL, showing the estimated values of the nonlinear parameters and the standard deviations of the residuals of the expressions describing the backbone curves.The nonlinear parameters, α 1 , α 2 and α 3 , are shown in blue and the standard deviations, σ 1 and σ 2 , are shown in green.(For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)dashed-red lines respectively and blue dots show the data used for the identification of the parameters.One significant feature of these results is the detached backbone curve, denoted S i 1.This is referred to as an isolated backbone curve, and it can be viewed as an imperfect bifurcation of S2, caused by the asymmetry of this system.For further discussion of bifurcations in backbone curves see [20], and for discussion of isolated backbone curves see [33].It can be seen that the system identification procedure used no data describing the backbone curve S i 1; despite this, Fig. 7 shows that S i 1 using the real and estimated parameters are close.This demonstrates that the estimated parameters lead to accurate predictions of the backbone curves, even for responses that are very distinct from those represented by the data.

Table 4
The true and mean estimated parameter values, along with the relative error of the estimated values, for system 3NL.

Parameter
True -5.848 Â 10 À 4 -Fig.7. The backbone curves S1 and S2, along with the isolated backbone curve S i 1, for system 3NL.These have been calculated using AUTO-07p (a numerical continuation package) for both the true and estimated parameters, represented by solid-green and dashed-red lines respectively.Blue dots show the data used for the identification of these parameters.These results are shown in the projection of the common response frequency, Ω, against the amplitude of displacement of the second mass, X 2 .(For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)Fig. 8.The response of system 3NL when subjected to forcing at amplitude P T x ¼ ½0:002; 0:002.The solid-green line shows the response branch when the true nonlinear parameters are used, and the dashed-red line shows the response branch when the mean estimated parameters are used.These results are shown in the projection of forcing frequency, Ω f , against the amplitude of displacement of the second mass, X 2 .(For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

Forced responses for system 3NL
Although backbone curves are of great value for understanding the fundamental properties of a nonlinear system, the forced responses are generally of greatest importance in engineering.Therefore, to further test the significance of the errors of the estimated parameters, the responses of the system when subjected to forcing are considered.Fig. 8 shows the response of system 3NL when subjected to forcing at amplitude P T x ¼ ½0:002; 0:002, such that forcing is applied directly to the first linear mode.Responses using both the true and mean estimated parameters are represented by a solid-green and a dashed-red line respectively.These results have been calculated using AUTO-07p [31].As the forcing is only in the first linear mode, the response envelops the backbone curve S1 (in which the first mode is dominant) leading to a Duffing-like forced response when viewed in this projection.It can be seen in Fig. 8 that the forced response curve for the estimated parameters is very close to that which uses the true parameters, suggesting that the errors have little influence for the forcing considered here.
One significant property of an isolated backbone curve, such as S i 1 in Fig. 7, is that it may correspond to isolated forced solutions, or isola.These can be difficult to detect through experiment, or through modelling techniques that rely on continuation without prior knowledge of their existence, and hence isolas present significant challenges and potential risks in engineering [33,34].Fig. 9 shows the forced responses of system 3NL using the estimated and true nonlinear parameters, represented by dashed-red and solid-green lines respectively.These are subjected to forcing at amplitude P T x ¼ ½0:01; À 0:01, such that the forcing is applied directly to the second linear mode, and it can be seen that an isola exists, as indicated by the existence of S i 1. Fig. 9 also shows that the forced response branches are predicted accurately by the estimated parameters, even for those describing the isola, which exists well-away from the data used for estimation.
As, in this projection, the isola describes a high-amplitude response, predicting that its existence is important.Furthermore, reaching the isola experimentally would prove challenging, and hence, even if its existence is known, experimental data describing an isola are often not available for identification.Therefore, it is an important feature of identification methods to predict parameters with sufficient accuracy that they may, in turn, be used to predict the existence of such features, well-away from the available data.

Conclusions and future work
In this work, an approach for the Bayesian identification of the nonlinear parameters of dynamic structures has been introduced.This approach utilises the backbone curves of a system, which may be obtained experimentally, and modelled using analytical expressions derived using the second-order normal form technique.These analytical expressions are computationally efficient to evaluate, hence the identification procedure is carried out extremely rapidly.As such, the probabilistic estimation of the nonlinear parameters of two example cases demonstrated herewhere each system has three unknown nonlinear parametersis completed in a matter of seconds.
As discussed, one disadvantage of this approach is that the analytical expressions provide an approximate description of the backbone curves, thus introducing error into the model of the system.This is demonstrated in the results presented here, as the true values of the parameters lie outside of the range of estimated parameters.It is also demonstrated, however, that the resulting errors are small and their influence on the dynamic responses are negligible.This is shown through a comparison of the backbone curves of the systems using both the true and estimated parameter values, and shows a strong Fig. 9.The response of system 3NL when subjected to at amplitude T x ¼ ½0:01; À 0:01.The solid-green line shows the response when the true nonlinear parameters are used, and the dashed-red line shows the response branches when the mean estimated parameters are used.These results are shown in the projection of forcing frequency, Ω f , against the amplitude of displacement of the second mass, X 2 .(For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)agreement at responses well-beyond those represented by the experimental data.Specifically, an isolated backbone curve is predicted for one case, and it is shown that this corresponds to an isola in the forced response of the system.
Not only has it been shown that this approach is extremely fast and accurate, it also lends itself to application to much larger systems.This is due to its use of analytical models, which remove the need for expensive numerical solving routines, and the use of the residuals of the expressions, which overcome the need to find explicit solutions.Additionally, the secondorder normal form technique, used to find the analytical expressions for the backbone curves, has a matrix-based formulation that is well-suited to computer automation, hence allowing this approach to be extended to larger, more complex systems.
The matrix β is now introduced in order to determine which of these terms are resonant (i.e.jω ℓ j ¼ Ω).The fn; ℓg th element of β is written as It can be seen that any element in β with a value of zero must correspond to an element in ½n q that describes the coefficient of a resonant term.Hence, one may define the matrix ½n u which is equal in size to ½n q and whose elements are populated by the coefficients of the resonant terms, such that ½n u nℓ ¼ ½n q nℓ if β nℓ ¼ 0; 0 i fβ nℓ a0: From this, the vector of nonlinear termssee Eq. (A.3)is defined using where α p ¼ ðα 1 þα 3 Þ=2m, α m ¼ ðα 1 À α 3 Þ=2m and α ¼ α p þ 8α 2 =m.Now, from Eqs. (A.8) and (A.9), and using N u ¼ ½n u u Ã , N u may be written as Substituting Eq. (A.10) into the resonant equation of motion, Eq. (A.3), it can be seen that all terms are sinusoidal and resonating at frequency Ω. Therefore it follows that 11) where I is the identity matrix, and u p , u m , N þ u and N À u form two sets of complex conjugate vector pairs, given by U 2 e þ jϕ 2 ! ; and (A.12) From Eq. (A.11), it can be seen that where from Eq. (A.10) where ϕ d ¼ ϕ 1 Àϕ 2 .Substituting Eqs.(A.12) and (A.14) into Eq.(A.13) gives Taking the imaginary parts of Eqs.(A.15) shows that sin ϕ d À Á ¼ 0. Therefore e 7 j2ϕ d ¼ 1 and e 7 jϕ d ¼ 71 ¼ p, such that when p ¼ þ1 the two modes are in-phase, and when p ¼ À1 the two modes are in anti-phase, as seen in Section 2.2.This allows Eqs.(A.15) to be written as where α p ¼ ðα 1 þ α 3 Þ=2m and α m ¼ ðα 1 À α 3 Þ=2m have been used.

Fig. 1 .
Fig.1.A schematic diagram of an in-line, two-degree-of-freedom oscillator with a symmetric linear structure and three nonlinear cubic springs.

Fig. 3 .
Fig. 3. Histograms, for system 1NL, showing the estimated values of the nonlinear parameters and the standard deviations of the residuals of the expressions describing the backbone curves.The nonlinear parameters, α 1 , α 2 and α 3 , are shown in blue and the standard deviations, σ 1 and σ 2 , are shown in green.(For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

Table 1
The parameter values used in the two example systems, 1NL and 3NL.