An investigation into model extrapolation and stability in the system identification of a nonlinear structure

Estimating a nonlinear model from experimental measurements of a vibrating structure remains a challenge, despite huge progress in recent years. A major issue is that the dynamical behaviour of a nonlinear structure strongly depends on the magnitude of the displacement response. Thus, the validity of an identified model is generally limited to a certain range of motion. Also, outside this range, the stability of the solutions predicted by the model are not guaranteed. This raises the question as to how a nonlinear model derived using data from relatively low amplitude excitation can be used to predict the dynamical behaviour for higher amplitude excitation. This paper focuses on this problem, investigating the extrapolation capabilities of data-driven nonlinear state-space models based on a subspace approach. The experimental vibrating structure consists of a cantilever beam in which magnets are used to generate strong geometric nonlinearity. The beam is driven by an electrodynamic shaker using several levels of broadband random noise. Acceleration data from the beam tip are used to derive nonlinear state-space models for the structure. It is shown that model predictions errors generally tend to increase when extrapolating towards higher excitation levels. Furthermore, the validity of the estimated nonlinear models become poor for very strong nonlinear behaviour. Linearised models are also estimated to have a complete view of the performance of each candidate model for each level of excitation.

experimental vibrating structure consists of a cantilever beam in which magnets are used to generate strong geometric nonlinearity. The beam is driven by an electrodynamic shaker using several levels of broadband random noise. Acceleration data from the beam tip are used to derive nonlinear state-space models for the structure. It is shown that model predictions errors generally tend to increase when extrapolating towards higher excitation levels. Furthermore, the validity of the estimated nonlinear models become poor for very strong nonlinear behaviour. Linearised models are also estimated to have a complete view of the performance of each candidate model for each level of excitation.

Introduction
System identification is an essential tool for understanding and modelling the behaviour of engineering structures, facilitating the mathematical description of real-life structures based on measured data. Historically, this process has been used to estimate linear models, taking advantage of linear theory based on the superposition principle [1]. However, many engineering structures are only linear to a first approximation and some do not behave linearly at all. The development of methods for system identification of such structures has gained more attention in recent years, driven by the need to capture inherent nonlinear effects [2,3]. They can be used to estimate a linearised model for a specific working condition, or to estimate a model to capture the full nonlinear characteristics of the structure. The first case generally aims to determine an equivalent linear model [4,5], sometimes referred to as the best linear approximation (BLA) [3,6]. This kind of model is generally reliable in the case of weak nonlinearity and only under the conditions of the training dataset, because the dynamical behaviour of a nonlinear structure significantly depends on the amplitude of vibration. A nonlinear model, however, has a wider region of validity, and can provide greater understanding of the system behaviour, including potential bifurcations in the response.
Nonlinear system identification is still a challenging task, even though several methods have been developed by the research community [7]. The identification process generally involves three steps: nonlinear detection, characterisation, and estimation. The first two steps can be addressed using ad-hoc methods or prior knowledge of the system [2], while the last step involves estimation of the model parameters from the experimental data. Except for the theoretical case when a perfect match exists between the actual structure and the identified model, the outcome of this step is related to the range of motion (vibration level) covered by the data itself. The reason is twofold: (i) even in the case of a successful estimation of a model structure, the true nonlinear phenomena might change when moving out of the training dataset; (ii) the estimated nonlinear model itself can behave differently when extrapolated, for instance in the case of a polynomial representations. This is why extrapolation from nonlinear models is generally not recommended [8]. Ideally, nonlinear system identification should be performed using training data that covers the operational range of motion of the structure under test. This is not always possible in real-life scenarios, and knowledge of the limitations of an identified model structure is crucially important information. This topic was first investigated in [9], and for extrapolated nonlinear models estimated using a polynomial nonlinear state-space representation [10].
This paper further investigates this subject by considering the extrapolation of data-driven statespace models estimated using a subspace approach [11,12]. State-space formulations are widely used in engineering applications, and the nonlinear problem addressed in this paper derives directly from the standard linear one. This makes the use of a subspace approach to estimate the system matrices quite convenient, since the essential definitions of classical subspace methods remain valid. A cantilever beam with geometric nonlinearity caused by the addition of magnets, is excited using an electrodynamic shaker with several levels of random excitation. This causes a response ranging from quasi-linear to strongly nonlinear. Linearised and nonlinear models are estimated from the measured data, and extrapolation (and interpolation) issues are discussed in both cases. In particular, linearised models show acceptable results only for weak nonlinear behaviours and around the selected working point, while they hugely fail in case of stronger nonlinearity or when trying to interpolate/ extrapolate. Nonlinear models generally outperform linear ones, but extrapolation remains an important limitation. To the authors' best knowledge, an experimental investigation on this topic comprising several excitation levels and a strong nonlinearity is not present in the scientific literature yet. The investigation described in this paper offers a useful perspective on experimental and methodological aspects. The intention is to provide comprehensive insight on the importance of accurately selecting an appropriate range of motion when designing the experimental setup. The results provide useful insight into the reliability and the limits of the estimated models both in the case of interpolation and extrapolation.
The paper is organised as follows. Section 2 presents the methods adopted in this paper and describes the experimental test rig. Section 3 presents the results obtained considering first linearised models and then nonlinear models. Finally, the conclusions of the present work are summarised in Sect. 4. Appendix A provides a numerical demonstration of the proposed methodology.

Nonlinear state-space modelling using subspace identification
The nonlinear state-space approach is used in this paper to model the dynamic behaviour of a nonlinear structure. Different formulations of this have been proposed in the literature, such as [10][11][12][13]. In particular, the formulation in [11] is the cornerstone of the Nonlinear Subspace Identification (NSI) technique, and is adopted in this paper. NSI has been used previously to estimate nonlinear models from experimental measurements considering discrete [14][15][16] and distributed [17] nonlinearities, as well as doublewell systems [18,19]. NSI treats the nonlinear restoring force of the system as a feedback force on the underlying-linear system (ULS). Considering the discrete case of a mechanical system with N degreesof-freedom, yields the equation of motion where M, C v and K 2 R NÂN are the mass, viscous damping, and stiffness matrices respectively, while y t ð Þ and fðtÞ 2 R N are the generalised displacement and external force vectors respectively. The nonlinear part of the equation is described by the term f nl y; _ y ð Þ 2 R N , which generally depends on displacements and/or velocities. It is assumed that f nl can be decomposed into J distinct nonlinear contributions using a linear-in-the-parameters model, yielding Each term of the summation is defined by a coefficient l j , a nonlinear basis function g j y; _ y ð Þ and a location vector L j 2 R N . The elements of L j define the position of the j th nonlinearity, and can have values of À1, 1 or 0 and. Introducing the extended-input vector f e as and the state vector x ¼ y T ; _ y T Â Ã T , a discrete statespace formulation can be written down as where s is the sampled time and matrices A; B e ; C; D e are the state, extended input, output and extended direct feedthrough matrices, respectively. These matrices can be estimated by using a subspace approach and assuming knowledge of the nonlinear basis functions g j and the location vectors L j of the nonlinear terms. The reader is referred to [11,20] for extensive information on this process. Once the statespace matrices are estimated, the model described by Eq. (4) can be used to estimate both the underlyinglinear and nonlinear features of the system. The first step involves: • Estimating the modal parameters of the ULS by eigenvalue decomposition of hte matrix A • Estimating the FRFs of the ULS by first defining the extended-FRF G e x ð Þ where I 2N is the identity matrix of size 2N and i is the imaginary unit. The matrix G e x ð Þ has the same structure of the extended force vector f e , so that its first block G x ð Þ is the FRF matrix of the underlying-linear system: Note that this step provides full knowledge of the underlying-linear dynamics of the system without the need for a linear measurement.
The next step involves the determination of the coefficients of the nonlinear terms l j . They are frequency-dependent and complex-valued, and are determined from the remaining blocks of G e . However, the true coefficients should be real numbers with no dependence on frequency. This only occurs in complete absence of noise and nonlinear modelling errors. When real measurements are performed, there are generally non-zero imaginary parts which are much smaller than their corresponding real parts. The reader is referred to [17,19,20] for detailed information on this step.
Finally, the estimated state-space model of Eq. (4) can be used to predict the response of the system to a given input vector fðtÞ or initial state vector. This can be carried out by propagating Eq. (4) over time, and it must be performed recursively for each time step due to the nonlinear nature of the problem. If the input vector fðtÞ is from a validation set, the outputs generated by the method y sim should replicate the measurements y val . The deviation between the simulated and measured outputs can be generally attributed to noise and/or nonlinear modelling errors. The output simulation percentage error n can be defined for each output n as which represents the capability of the estimated model in predicting the system dynamics. This work focuses on the last point-to understand the limitations of this process due to the nonlinearity. A numerical example is proposed in Appendix A: numerical demonstration to validate the methodology in a case where both the true underlying-linear FRF and nonlinear restoring force are known. The numerical example consists of a two-degrees-of-freedom system with a non-smooth asymmetrical nonlinearity. Monte Carlo simulations have been performed to account for stochastic variations and model uncertainties.

Experimental test rig
The test rig consists of a cantilever beam with total length of 400 mm and rectangular cross-section (width of 25 mm and thickness of 5 mm). There are three small neodymium magnets located at the tip as depicted in Fig. 1. Different nonlinear effects can be obtained by changing the polarization of the magnets such as hardening, softening, multi-stable or quasizero-stiffness [21][22][23]. A repelling configuration was considered in this work, generating a predominantly hardening nonlinearity. The system was excited using an electromagnetic shaker (Modal Shop -model K2007E01) that provides brown noise to the structure. This allows to have higher intensity at lower frequencies so that the analysis can be focused on the first mode. Six levels of excitation were considered from 0.08 N RMS to 1.16 N RMS following a geometric progression, as listed in Table 1. The excitation level was increased roughly by 70% from one level to the subsequent. An accelerometer A 1 (model PCB 352A24) was used to measure the response of the beam at the tip, the sampling frequency was 2000 Hz and each acquisition lasted 120 s. All the signals were recorded using a NI USB-4431 acquisition board from National Instrument. The last 20 s of the measured data are used as validation set in all the identification steps involved in the following sections, while the rest is adopted as a training set.
The measured accelerations in the time domain are depicted in Fig. 2a for each level of excitation. The corresponding displacements are obtained by integrating twice the accelerations and depicted in Fig. 2b. Note that the overdots denote differentiation with respect to time.
The experimental FRF obtained using the H1estimator and the corresponding coherence function are depicted in Fig. 3 for each level. A clear shift of the resonance peak to the higher frequencies can be observed for higher excitation levels, which is a sign of strong nonlinear hardening behaviour. Also note the diminishing level of coherence as the excitation level increases, which is an additional indication of the presence of a strong nonlinearity.

Best linear approximation (BLA) via subspace identification
A linear state-space model is estimated for each level of excitation using the linear subspace identification technique [24] with a model order of 2. To ensure that the estimated model is the best possible linear candidate, an a-posteriori optimisation is carried out using the validation set. This serves to estimate the BLA of the selected data, and is used in the extrapolation process of Sect. 3.4 as a comparative measure.
To this end, the system output is generated using the estimated linear state-space matrices (A, B, C, D) and the measured forcing input. A least-square minimisation problem is defined to reduce the residual between the measured output y val 1 and the simulated output y sim 1 by optimising the state-space matrices. The minimisation problem is defined by Þ and the vector operation vec Á ð Þ stacks the coloumn of a matrix on top of each other. The starting point of the optimisation is the set of matrices obtained from the subspace approach and the Levenberg-Marquardt algorithm is adopted with a step tolerance of 10 -8 .
The system FRFs for the different excitation levels of Fig. 3 are now depicted in Fig. 4 together with the corresponding BLA. The estimated natural frequencies, damping ratios and simulation errors are listed in Table 2. The output simulation errors are computed according to Eq. (7). It is clear that an increase in the excitation level is associated with a decrease in the BLA performance. The output simulation errors become extremely high from level 3 (light green line in the figure), denoting the inability of a linear model to capture the system dynamics in the case of severe nonlinear behaviour.
Since the model obtained using the BLA is a linearisation for a specific working range, it is generally not advisable to use this model away from  the training data set when the system behaves nonlinearly [3]. This is confirmed by the error matrix shown in Fig. 5, obtained from the output simulation error 1 when extrapolating and interpolating towards higher and lower levels, respectively. The values on the diagonal correspond to those reported in Table 2. The off-diagonal values dramatically increase when moving out of the BLA training set both in the case of interpolation (blue background) and extrapolation (red background). This result is expected for the latter case, but it is not so obvious in the case of interpolation. The reason can be seen by close examination of Table 2 and Fig. 4. The estimated natural frequencies and damping ratios of the linearised models change with the excitation level, and this occurs both when moving towards higher levels (i.e. extrapolating) or lower ones (i.e. interpolating). As an example, the error of the model BLA3 at excitation level 3 is 25.5% (row 3, column 3). The same model gives an error of 51.5% when reducing the excitation to level 2 (row 3, column 2), while it gives an error of 64.4% when increasing the excitation to level 4 (row 3, column 4). Consequently, the minimum errors for each identification level are those on the diagonal, corresponding to the excitation level at which the BLA is calculated. A surface plot of the error matrix is depicted in Fig. 6 to highlight this behaviour.

Qualitative nonlinear characterisation
The nonlinear behaviour of the system is first investigated using the Acceleration Surface Method (ASM) [25] to get a qualitative data-based visualisation of the system restoring force. Since the sensor is located close to the source of the nonlinearity (i.e., the  magnets), and is quite far from the excitation point, its restoring surface can give a rough visualisation of the system characteristic. In particular, a qualitative representation of the restoring force can be obtained by slicing the acceleration surface around the lowvelocity region. The points of the restoring force are therefore obtained as the ones that satisfy the condition where the tolerance e v is set to 10 À3 . The result is shown in Fig. 7. It can be seen that the system has a hardening restoring force, and there is consistency for the different excitation levels.
A polynomial representation of the nonlinear restoring force is therefore used to apply the NSI technique in the following section, considering the nonlinear basis functions in which maximum exponent P ¼ J þ 1 is selected independently for each excitation level to minimise the errors over the output residual given in Eq. (7).

Estimation of the nonlinear model parameters
A nonlinear state-space model of order 2 is estimated for each level of excitation using the nonlinear subspace identification technique. The following information is estimated for each level: 1. The modal parameters and FRFs of the ULS. These are compared with the lowest level of excitation (0.08 N RMS ), which is the closest to linear behaviour. The FRFs of the ULS are depicted Fig. 8a, and the natural frequencies and damping ratios are listed in Table 3. 2. The coefficients l j of the nonlinear basis functions g j and the nonlinear restoring force of the structure. The final shape of the estimated nonlinear restoring force is depicted in Fig. 8b for each level of excitation. 3. The errors over the output residual 1 as in Eq. (7).
To this end, the last 20 s of each test are used again as a validation set, and the estimated state-space model is adopted to generate the outputs given the measured forcing input.
The identification is repeated for each level considering different sets of candidate basis functions having the maximum exponent P Eq. (10) ranging  Table 3 for each test.
The results of the identification show significant consistency up to the fifth level (0.64 N RMS ) on both the underlying-linear parameters and the nonlinear estimation. Residuals and errors increase noticeably at level 6 (0.81 N RMS ), suggesting that the selected nonlinear basis functions are not appropriate when describing the nonlinear behaviour of the system at high-level responses. This is likely to be associated with the occurrence of new nonlinear phenomena such as geometric nonlinearity caused by large-amplitude displacements [17], which are not well captured by the polynomial functions used. The model estimated from level 6 is therefore considered as erratic and marked with the symbol Á Ã in Table 3.
Interestingly, the output simulation error of lowest excitation level (Level 1) is higher than the subsequent ones. A possible explanation is that vibrations are so low in this case that noise and/or boundary connections affect the result. A similar behaviour is observed and discussed in Sect. 3.4.

Extrapolation, interpolation and model stability
The estimated state-space models have been used above to replicate the outputs of the system considering a validation set belonging to the same excitation level of the training set. It is, however, interesting to analyse the errors of the residuals when extrapolating and interpolating towards higher and lower levels, respectively. The results are illustrated in the error matrix of Fig. 9. The best linear approximation (BLA) is also reported for each level as a comparative measure, considering the diagonal values of Fig. 5. The linearised model shows similar performances of the nonlinear models 1-5 when considering the lowest excitation level (first column), confirming the quasi-linear behaviour in this range of motion. Interestingly, it outperforms nonlinear model 6 (marked with Á Ã ), which proved to be poor as indicated in Table 3. As for the nonlinear models, the errors tend to increase when extrapolating towards higher excitation levels (red background). This confirms the idea that extrapolation from nonlinear models is generally not advisable, especially in the case of polynomial expansions. Cells with a cross correspond to the cases  The Á* symbol indicates an erratic model when the simulation does not converge, making the extrapolation process not possible at all. Assuming that the underlying-linear state-space model is stable, this may occur for two reasons: (i) the extrapolated nonlinear restoring force becomes too large for convergence at a certain time step; (ii) the extrapolated nonlinear restoring force exhibits unstable behaviour, thus making the nonlinear model unstable. The latter case can be observed in Fig. 10 which shows the case when extrapolating from level 3 towards level 5. The origin of this behaviour is in the high polynomial terms of the selected nonlinear basis functions and must be carefully investigated when performing similar nonlinear data-based identifications.
It is worth highlighting that the estimated nonlinear models outperform the linear models in all the other cases, especially when no interpolation or extrapolation is involved (diagonal values). An interesting result can be observed by looking at the first column (lowest excitation level). The errors are generally higher than the subsequent columns for all the valuable nonlinear models (NSI1 7 NSI5). As mentioned in the previous section, this is possibly caused by a low signal to noise ratio, or to the occurrence of low-amplitude nonlinear effects related to the boundary connection of the beam, that are not included in the model.

Conclusions
The extrapolation and interpolation behaviour of datadriven state-space models from a nonlinear structure has been investigated in this paper. Analysis was carried out on a bench top experimental test rig consisting of a cantilever beam in which magnets were used to generate strong geometric nonlinearity. The beam was driven by an electrodynamic shaker using several levels of broadband random noise. Its behaviour ranged from quasi-linear to strongly nonlinear, and the data collected were used both as training and validations sets for subspace-based identification techniques (linear and nonlinear). The limitations of the identification process have been evaluated by comparing model predictions and measurements across the different levels of excitations. This was performed by defining first a measure of error based on the RMS deviation of the predicted and measured validation time histories. An error matrix was assembled which presents interpolation and extrapolation issues in one single view. This facilitates a clear picture of the limits of the estimated models when trying to simulate the behaviour of the structure in conditions that exceed the range of motion of the training data set. The results obtained with nonlinear models are also compared with linearised ones to have a complete view of the performance of each candidate model for each level of excitation. It is shown that nonlinear model predictions errors generally tend to increase when extrapolating towards higher excitation levels. The results are in accordance with the expectations, and provide useful insight on the importance of accurately selecting an appropriate range of motion when designing an experimental setup. Furthermore, the validity of the estimated nonlinear models becomes poor for very strong nonlinear behaviour. This is clear from the identification performed at the highest excitation level, which can be related to the occurrence of new nonlinear phenomena not included in the model. Linearised models on the other side can represent a valid modelling tool provided that the nonlinear behaviour is weak and the operational range of motion is covered by the validation set, making them inadequate to both interpolation and extrapolation processes.
Funding Open access funding provided by Politecnico di Torino within the CRUI-CARE Agreement.
Data availability The datasets generated during and/or analysed during the current study are not publicly available but are available from the corresponding author on reasonable request.

Declarations
Conflict of interest The authors declare that no funds, grants, or other support were received during the preparation of this manuscript. The authors have no relevant financial or non-financial interests to disclose.

Open Access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A: Numerical demonstration
The system considered for the numerical simulations is depicted in Fig. 11. It has two degrees-of-freedom and a mechanical stop positioned 2 mm away from the second mass. This creates a strong, non-smooth asymmetric nonlinear behavior of contact type when the gap is filled.
The theoretical nonlinear restoring force is therefore given by the following relationship: The numerical example and the system parameters are taken from [16], except for the gap value d ¼ 2mm that is assumed to be known in this paper with a stochastic uncertainty, as described later. The system is excited at DOF 1 by a zero-mean Gaussian random input considering 10 equally spaced levels between 0.5 N RMS and 5 N RMS. The outputs are corrupted by 3% Gaussian noise and 100 Monte Carlo simulations are conducted for each level of excitation. For each simulation, the numerical integration of the equation of motion is performed using the Newmark method with a sampling frequency of 4096 Hz and considering a time span of 60 s. The first 40 s of the acquisition length are used to estimate the state-space model with NSI and the last 20 s are used as validation set. The uncertainty in the gap value is accounted for in the selection of the nonlinear basis functions of NSI that read: where n is sampled from a zero-mean Gaussian distribution having 3r ¼ 0:2mm. Therefore, the adopted gap value d Ã is different for each Monte Carlo simulation with a maximum variation of AE0:2mm from the nominal value (with a confidence interval of 99.7%). The errors on the estimated modal parameters and nonlinear coefficients are listed in Table 4 for each level of excitation. Results are averaged over the Monte Carlo simulations, with standard deviations written in brackets. The errors on the modal parameters of the ULS remain quite low among the excitation levels and are generally higher for the first mode of the system. The reason is that the nonlinearity mainly affects the first mode. This is clear from Fig. 12, showing the H1-estimator for several excitation levels. As for the nonlinear coefficients k 1 and k 3 , errors and standard deviations are very high for the lowest excitation levels, while they tend to stabilize for higher levels. The first level (0.5 N) does not provide any estimation of the nonlinear coefficients because the output displacement y 2 is not high enough to cover the gap value d, therefore, the nonlinearity is not activated. Moreover, it can be observed that errors on k 1 are generally higher than the errors on k 3 . This is to be addressed to the uncertainty on the gap estimation, that has a higher impact on the estimation of the stiffness k 1 , especially at low levels.
The estimated state-space nonlinear model of each Monte Carlo simulation is used to replicate the outputs of the system for all the excitation levels. The output   Fig. 13 Error matrix on the validation set of DOF 2 averaged over MC simulations. The values refer to the percentage errors defined in Eq. (7). The darkness of the cells is proportional to the error level. Blue background (lower-triangular matrix) indicates simulations with interpolation; red background (upper-triangular matrix) indicates simulations with extrapolation; bold type indicates same-level simulations higher than 1 for the considered system, only the error matrix of 2 is depicted in Fig. 13. The values on the diagonal correspond to the same-level identification and generally correspond to the lowest error for each model, as expected. Errors tend to increase when extrapolating (red background) or interpolating (blue background), with the exception of the first excitation level (first column). The reason is that it behaves linearly, therefore a good estimation of the ULS is enough to achieve acceptable errors. For the same reason, the model NSI1 generates huge errors when extrapolated towards every other excitation level.