In the present work, two issues that can complicate a damage detection process are considered: the uncertainties and the intrinsically nonlinear behavior . To deal with these issues, a stochastic version of Volterra series is proposed as a baseline model, and novelty detection is applied to distinguish the condition of the structure between a reference baseline state (presumed “healthy”) and damaged. The studied system exhibits nonlinear behavior even in the reference condition, and it is exposed to a type of damage that causes the structure to display a nonlinear behavior with a different nature than the initial one. Additionally, the uncertainties associated with data variation are taken into account in the application of the methodology. The results conﬁrm that the monitoring of nonlinear coefﬁcients and nonlinear components of the system response enables the method to detect the presence of the damage earlier than the application of some linear-based metrics. Besides that, the stochastic treatment enables the speciﬁcation of a probabilistic interval of conﬁdence for the system response in an uncertain ambient, thus providing more robust and reliable forecasts. Keywords Uncertainties, damage detection, stochastic Volterra model, nonlinear behavior Introduction Structural Health Monitoring (SHM) techniques aim to reduce the maintenance cost and increase the reliability and security of aerospace, civil, or mechanical engineering structures1. Moreover, within the hierarchy of complexity that SHM methodologies may achieve, damage detection is the ﬁrst step, and its performance is fundamental to the success of the subsequent application of higher forms of SHM in the hierarchy2. In this sense, many authors have studied and developed various damage detection techniques to be implemented for different structures and applications2–7. Of course, there is no general approach that can be used to detect damage in all real systems, since the ubiquitous challenge when the intrinsically nonlinear behavior of many systems8and the data variation related to uncertainties9,10is counted in the analysis. Many otherwise-linear structures can exhibit nonlinear phenomena induced by the presence of damage, and in this situation, any form of nonlinearity detector is akin to detecting damage11. Such damage that produces nonlinear behavior, for example, delamination12,13, rubbing and unbal- ance in rotor systems14–16, and opening cracks17–19may be detected through the observation of nonlinear behavior in the measured responses. However, as mentioned before, many structures fundamentally present nonlinear behavior even in the reference condition20, causing confusion in the nonlinearity (as a proxy for damage) detection process8. Additionally, systems’ measured output can exhibit data variability from sources such as environmental or input load variation, aleatoric noise, changes in boundary condi- tions, variations in the fabrication processes (i.e., unit-to-unit variability), and others21–23. These variations all confoundthe damage detection process, suggesting the use of prob- abilistic tools24–26, regression models5,27, machine learning algorithms28, probabilistic model selection approaches29–32, outlier analysis33,34and novelty detection methods5. Con- sidering the data variation problem is important in reducing the number of false alarms10,35, although there are situations where variability can mask positive detections as well. Villani et al. (2019)36introduced an extension of the deterministic V olterra series approach used by Bruce et al. (2017)37to detect damage in intrinsically nonlinear systems based on input/output measurements, recognizing data variation. The authors examined the identiﬁcation of the nonlinear data-driven model several times, using Monte Carlo simulations, to create a stochastic reference model capable of predicting the nonlinear performance and the ﬂuctuations in the response at the same time. Two main strategies were suggested to detect the presence of a crack in an initially nonlinear beam, considering simulated data, based on the random kernel contributions and random kernel coefﬁcients. The results presented an adequate performance in distinguishing the intrinsically nonlinear 1UNESP - Universidade Estadual Paulista, Faculdade de Engenharia de Ilha Solteira, Departamento de Engenharia Mec ˆanica, SP , Brasil 2UERJ - Universidade do Estado do Rio de Janeiro, NUMERICO – Nucleus of Modeling and Experimentation with Computers, RJ, Brasil 3UCSD - University of California San Diego, Department of Structural Engineering, La Jolla, CA, USA Corresponding author: Luis G. G. Villani, UNESP - Universidade Estadual Paulista, Faculdade de Engenharia de Ilha Solteira, Departamento de Engenharia Mec ˆanica, Av. Brasil, 56, Ilha Solteira, 15385-000, SP , Brasil. Email: luis.villani@unesp.br Prepared using sagej.cls [Version: 2017/01/17 v1.20] 2 Journal Title XX(X) behavior and the data variation from the nonlinearity caused by the damage. Although the simulated results showed a promising performance of the method, experimental practice suggested the use of a stochastic V olterra expansion to explore the model space, particularly in the more challenging present problem of inherent nonlinearity, rather than induced nonlinearity. The authors also hypothesized that a formulation using the kernels’ coefﬁcients and contributions concomitantly in the same index could be attractive to enhance the robustness of the method, leaving the method to detect damage with diverse characteristics without loss of performance. Additionally, in the simulations performed, only the variation of linear parameters was considered in the data variation scenario, a simpliﬁcation that might not reﬂect the behavior of structures that operate in the nonlinear regime of motion under the presence of uncertainties. This simpliﬁcation caused a low variation of the high-order kernels coefﬁcients, suggesting that a more realistic application encompassing the variation of the nonlinear components could complicate the application of the approach. Then, in Villani et al. (2019)38the authors showed the experimental application of the approach based on the V olterra kernels’ contributions to detect damage in an intrinsically nonlinear beam, considering data variation related to the reassembly of the experimental setup. However, the damage simulated was associated with the loss of mass that reﬂected in the variation of the natural frequencies of the equivalent linear system. In this situation, the damage did not have inﬂuence in the estimation of the kernels’ coefﬁcients. These results pointed to the demand for developing a hybrid approach capable of making use of the kernels’ contributions and coefﬁcients together in a robust index. Moreover, an experimental application considering damage with direct inﬂuence on the nonlinear components of the response could improve the performance of the approach with regard to the differentiation between the intrinsically nonlinear behavior and damage-induced nonlinear behavior. Hence, this paper aims to cover issues that were not yet considered in both previous works: (a) an experimental application of the stochastic V olterra series methodology with the presence of damage that produces nonlinear behavior to the system - a breathing crack emulation; (b) the observance of variation in the intrinsic nonlinearity of the structure including the data variability, thus not only considering its realization in the linear components; (c) the generalization of the approaches studied before with the development of a new hybrid method that considers both the kernels’ coefﬁcients and contributions simultaneously in the damage index; and, ﬁnally, (d) the construction of a theoretical distribution of the damage index calculated in the reference condition to reduce the number of experimental realizations needed to estimate the threshold value based on the Kernel Density Method used before (a practical problem when we consider real-world experiments). To the best knowledge of the authors, this is the ﬁrst paper that assumes nonlinear changes associated with an experimental mechanism of damage with the assumption of reference already nonlinear, but unlike Born et al. (2010)8, here considering the inherent uncertainties in the experimental setup to perform a rigorous stochastic SHM method.In this context, an intrinsically nonlinear beam is analyzed, with natural data variation in the full experimental context, to investigate the performance of the proposed methodology. The variation simulated in the data reﬂects merged changes both in linear and nonlinear components of the system response. Furthermore, the novelty detection is reshaped to hold the random kernels’ contributions and coefﬁcients in the same damage index, using Principal Component Analysis (PCA) and Mahalanobis distance metrics. A formal hypothesis test is presented to create a more robust damage detection methodology, based on a theoretical distribution for the Mahalanobis distance determined in the reference condition of the structure. The results exposed in this work have demonstrated the beneﬁcial performance of a nonlinear metric to detect damage in this situation and the capability of the stochastic V olterra series to predict the data variation in a probabilistic framework, improving the statistic conﬁdence of the method. To be familiar with the V olterra series model application in damage detection problems, ﬁnd the motivation to the study of the nonlinear phenomena in the procedure, and have more details about the V olterra series reformulation to predict the nonlinear responses considering data variation, the interested reader is referred to seeing Bruce et al. (2017)37, Villani et al. (2019)36and Villani et al. (2019)38. The present paper will be centered on the primary differentiation between the approach used here, and the one analyzed previously and in the experimental results obtained with the application of the method contemplating an initially uncertain nonlinear system subjected to the presence of damage that provokes nonlinear phenomena to the structure behavior, that was not done before. The present paper is organized as follows. Section “The damage detection methodology based on stochastic V olterra series” describes the stochastic mathematical model used to describe the nonlinear systems response and the approach proposed to detect the damage considering the data variation related to uncertainties. Section “Experimental setup” shows the nonlinear structure considered in this work, the damage simulated and the main characteristics of its behavior. Section “Application of the proposed methodology” shows the application of the methodology proposed and the main results obtained. Finally, the section “Final remarks” presents the conclusions of the work. The damage detection methodology based on stochastic Volterra series This section outlines the methodology proposed to be practiced in the damage detection problem in initially nonlinear systems, taking into account the data variation related to uncertainties. The development of the stochastic model is brieﬂy described, with more detail in Villani et al. (2018)36. Additionally, the reader can obtain more information about the deterministic V olterra series expanded using Kautz functions in Bruce et al. (2017)37. The stochastic version of the Volterra series In order to take into account the uncertainties in the model formulation, a parametric probabilistic approach is assumed. Thus, the model parameters are assumed as random Prepared using sagej.cls Villani et al. 3 parameters and the model response as a random process21–23. Therefore, a probability space (;;P)is considered, where represents the sample space, is a-algebra over , and Pis a probability measure36. In the discrete-time domain, assuming the presence of uncertainties, a single system output can be interpreted as a random process realization that is a consequence of a single deterministic input. The relationship between the deterministic input and the random output can be described, using the convolution notion39, through the stochastic version of the V olterra series y(;k) =1X =1N11X n1=0::: :::N1X n=0H(;n 1;:::;n)Y i=1u(kni);(1) where (;k)2Z+7!y(;k)represents the single ran- dom output that is consequence of the single deter- ministic input k2Z+7!u(k),(;n 1;::;n)2Z7! H(;n 1;:::;n)represents the random version of the - order V olterra kernel and Z+represents the set of integer positive numbers. The principal beneﬁt of the V olterra series model is the capability to reproduce the system output as a sum of linear and nonlinear contributions y(;k) =1X =1y(;k) = =y1(;k)|{z} linear+y2(;k) +y3(;k) +|{z} nonlinear;(2) where (;k)2Z+7!y1(;k)is the random output obtained using the ﬁrst random kernel, (;k)2Z+7! y2(;k)is the random output obtained using the second random kernel, and so on. In this work, the series will be truncated in the third order kernel because of the cubic characteristic of the nonlinear system response investigated, and the capability to classify linear and nonlinear contributions in the total response will be used in the damage detection procedure as feature sensitive to the presence of damage. On the other hand, as broadly addressed in previous works36,37,40, the central disadvantage of the approach is the challenge in achieving the convergence when a high number of terms is used. To solve this problem, the V olterra series can be extended utilizing the Kautz functions41,42, and the system random output can be represented as y(;k)1X =1J1X i1=1::: :::JX i=1B(;i1;:::;i)Y j=1l;ij(;k);(3) whereJ1;:::;Jrepresents the number of Kautz functions used in the kernels projections, the -order randomV olterra kernel expanded in the orthonormal basis is represented by the random process (;i1;:::;i)2 Z +7!B(;i1;:::;i)and the random process (;k)2 Z+7!lij(;k)is a ﬁltering of the deterministic input signal by the random Kautz functions. The Kautz functions are supposed random because their deﬁnition depends on the dynamics of the system and, as the system response is considered as a random process, it is presumed that the Kautz functions will randomly change too. Conclusively, the coefﬁcients of the kernels can be calculated using the least-squares approximation in a deterministic way37, and then, adopting Monte Carlo simulations, the process is repeated until the stochastic model converges. The Monte Carlo method was chosen because it is easier to perform when the deterministic algorithm is known23,43. More information about the random Kautz functions and the process of the random V olterra kernels estimation may be found in Villani et al. (2019)36,38. Damage detection based on novelty detection In the previous work, two points were analyzed separately in the damage detection process: the kernel coefﬁcients and the kernel contributions. Nevertheless, it is difﬁcult for practical applications to decide which feature is better. Thus, the damage detection index used here examines both features jointly. Regarding the stochastic model identiﬁed with training data in the reference condition, the damage- sensitive index may be determined in the reference status. First of all, recognizing that the V olterra series will be truncated in the third order kernel, the kernels coefﬁcients can be allocated together to be used as damage sensitive feature lin(;i1) =8 >>>< >>>:B1(;1) B1(;2) ... B1(;J 1)9 >>>= >>>;; qua(;i1=i2) =8 >>>< >>>:B2(;1;1) B2(;2;2) ... B2(;J 2;J2)9 >>>= >>>;; (4) cub(;i1=i2=i3) =8 >>>< >>>:B3(;1;1;1) B3(;2;2;2) ... B3(;J 3;J3;J3)9 >>>= >>>;; where the random process (;i1)2Z+7!lin(;i1) represents the coefﬁcients of the ﬁrst kernel, (;i1=i2)2 Z+7!qua(;i1=i2)represents the coefﬁcients of the diagonal of the second kernel and (;i1=i2=i3)2 Z+7!cub(;i1=i2=i3)represents the coefﬁcients of the main diagonal of the third kernel. Additionally, the contribution of the kernels can be calculated as Prepared using sagej.cls 4 Journal Title XX(X) ylin(;k)J1X i1=1B1(;i1)l1;i1(;k); yqua(;k)J2X i1=1J2X i2=1B2(;i1;i2)l2;i1(;k)l2;i2(;k); (5) ycub(;k)J3X i1=1J3X i2=1J3X i3=1B3(;i1;i2;i3)::: :::l3;i1(;k)l3;i2(;k)l3;i3(;k); where (;k)2Z+7!ylin(;k)is the linear contribu- tion, (;k)2Z+7!yqua(;k)is the quadratic con- tribution and (;k)2Z+7!ycub(;k)is the cubic contribution. So, to reduce the order of the classiﬁcation problem, the Principal Component Analysis (PCA)26,44–46 can be applied to the kernels contributions ylin(;k) PCA Clin(;1);:::;Clin(;npca); yqua(;k) PCA Cqua(;1);:::;Cqua(;npca);(6) ycub(;k) PCA Ccub(;1);:::;Ccub(;npca); where (;npca)2Z+7!Clin(;npca)represents the principal components of the linear contribution, (;npca)2 Z+7!Cqua(;npca)represents the principal compo- nents of the quadratic contribution, (;npca)2Z+7! Ccub(;npca)represents the principal components of the cubic contribution and npcais the number of principal components considered. The number of components was deﬁned based on the contribution of each component in the construction of the covariance matrix44. After the calculation of the kernels coefﬁcients and the principal components of the kernels contributions, the damage index can be deﬁned in the reference condition Ilin= [lin(;i1)Clin(;npca)](Ns(J1+npca)); Inlin= [qua(;i1=i2)cub(;i1=i2=i3)::: (7) :::Cqua(;npca)Ccub(;npca)](Ns(J2+J3+2npca)); whereIlinandInlin are, respectively, the linear and nonlinear indices in the reference situation. The linear and nonlinear indices will be assessed with the aim of comparison between the linear and nonlinear methodologies performance. So, in the reference condition, Ilinis aNs (J1+npca)matrix andInlin is aNs(J2+J3+ 2npca) matrix, being Nsthe number of observations used in the training phase of the reference stochastic model. Now, with the structure in an unknown (“test”) situation, a new deterministic model can be identiﬁed, and the indices may be estimated in the unknown status Ilin= [lin(i1)Clin(npca)](1(J1+npca)); Inlin= [qua(i1=i2)cub(i1=i2=i3)::: (8) :::Cqua(npca)Ccub(npca)](1(J1+npca));whereIlinandInlin are, respectively, the linear and nonlinear indices in an unknown condition. As the indices calculated in the reference condition are matrices composed by more than one single feature, the novelty detection has to be performed considering multivariate data. Therefore, the comparison between the indices calculated in the reference and an unknown condition can be made considering de squared Mahalanobis distance33,47 D2 m= [ImIm]T-1 Im[ImIm]; (9) wherem=linorm=nlin,ImandImare, respectively, the mean vector and the covariance matrix of the index calculated in the reference condition. The simple machine learning method based on Mahalanobis distance is used here because the classiﬁcation is done between two possible conditions (healthy and damaged) and the goal is to examine the performance of the V olterra kernels characteristics as damage sensitive features and not to investigate the differences between reﬁned classiﬁcation methods. Other classiﬁcation methods can be used in the future to improve the methodology depending on the real application confronted. With the squared Mahalanobis distance calculated, a hypothesis test may be proposed. In this work, this distance calculated in the reference condition is modeled with a chi- square distribution, which may be calculated based on the assumption of independence and normality in the underlying multivariate features from which the squared Mahalanobis distance is calculated48,49. Ideally, a sampling distribution of the Mahalanobis distance is desirable, but no such analytical form is known to exist, so the theoretical model is ﬁt to the (limited) empirical data obtained. This approximation has satisfactory performance, as will be shown further along. So, the hypothesis is proposed H0:D2 mX2; H1:D2 m X2;(10) whereX2is the chi-square distribution, H0is the null- hypothesis (healthy condition) and H1is the alternative hypothesis (damaged condition). Besides, the probability of a distance value calculated to be included in the theoretical chi-square distribution can be computed integrating its Probability Density Function (PDF)50 pm=F(D2 mj) =ZD2 m 0t(2)=2et=2 2=2(=2)dt; (11) wherepmis the probability of the value D2 mbelonging to the chi-square distribution, (:)is the Gamma function, and is the number of degrees-of-freedom. Finally, a sensitivity value can be determined depending on the application and probability of false alarms tolerated, and the hypothesis test can be rewritten H0:pm ; H1:pm< ;(12) where represents the sensitivity chosen for the hypothesis test. The deﬁnition of this parameter depends on the practical Prepared using sagej.cls Villani et al. 5 application and, as an experimental laboratory setup is utilized, several values will be examined to study the performance of the method. Figure 1 shows a ﬂowchart of the damage detection approach. On the left, it can be observed the training phase, with the identiﬁcation of the stochastic reference model and the estimation of the damage indices in the reference condition. On the right, it is represented the identiﬁcation of a new model in an unknown status and the calculation of the new indices. Then, the indices are correlated using the squared Mahalanobis distance, and ﬁnally, the hypothesis test is applied to classify the condition of the structure between healthy and damaged. Experimental setup The experimental setup used is presented in Fig. 2. The structure monitored is formed by a clamped-free beam, that is constructed by gluing four thin beams of Lexan together, 2:424240[mm3] each one, with the intention of emulating a damage propagation that is described further on. At the free boundary, two steel masses are afﬁxed and interact with a magnet, generating a nonlinear behavior in the system response, even in the reference condition due to added magnetic potential. Moreover, the setup includes: A National Instruments acquisition system: –CompactDAQ Chassis(NI cDAQ-9178); –A C Series Sound and Vibration Input Modules (NI-9234); –A C Series V oltage Output Module (NI-9263). Electrodynamic Transducer Labworks Inc. (ET-132); Ampliﬁer MB Dynamics SL500VCF); Load cell PCB PIEZOTRONICS (208C02); Accelerometer PCB PIEZOTRONICS (352C22). The electrodynamic transducer is employed to excite the structure with different signals and considering two levels of input (low - 1 V RMS and high - 6 V RMS). The output data is measured by the accelerometer positioned close to the free extremity of the beam, because the authors are only interested in the region of the ﬁrst mode shape of the structure (see Fig. 2). The input signal analyzed is the voltage signal applied in the electrodynamic transducer. As a single-input/single-output (SISO) model is considered, this pair of signals is enough to identify the V olterra kernels and monitor the structure health. All the acquisition parameters, signals considered, and equipment used were the same in the experiments performed considering the different structural conditions. Intrinsically nonlinear behavior The mechanical system used exhibits nonlinear operation even in the reference condition, without the presence of damage. Figure 3a shows the results obtained during the stepped sine test applied considering two levels of input. When the input applied has a low level of amplitude (1 V RMS), the output signal shows linear characteristics for both up-sweep and down-sweep inputs. However, when the input signal is at a sufﬁciently high amplitude (6 V RMS), the nonlinear phenomena can be seen with the jump presented in the test. Additionally, Fig. 3b shows the time-frequencydiagram of the system response considering a high level of a chirp input in the region of the ﬁrst mode shape of the structure. The presence of the cubic harmonic in the response conﬁrms the nonlinear characteristic of the response caused by the interaction between the magnet and the steel masses. So, the system studied presents nonlinear behavior, subject to sufﬁcient level of input, even when the structure is healthy. This characteristic is obtained because the magnetic force changes nonlinearly with distance from the end masses to the magnet. This intrinsically nonlinear behavior can be confused with the simulated damage that causes a distinct nonlinear characteristic to the system response, as will be seen next. Damage simulated The damage imposed on the structure aims to simulate a breathing crack present in the system. In this sense, four different beams were built to be used in the application of the damage detection methodology: Training beam: beam constructed with 4 intact Lexan beams and used in the training phase of the algorithm (see Fig. 4a); Test beam: beam constructed with 4 intact Lexan beams and used in the test phase of the algorithm (see Fig. 4a); Damage I: beam constructed with 3 intact and 1 cut beam (see Fig. 4a and b); Damage II: beam constructed with 2 intact and 2 cut beams (see Fig. 4a and b); The cut in the beams is positioned close to the excitation point (see in Fig. 2b). This spot was chosen to obtain the required nonlinear behavior to test the performance of the algorithm. The damage condition might be judged severely, but the position and excitation combined were deﬁned to generate a condition that is difﬁcult to detect – mainly in the condition Damage I – as will be shown further along. Figure 5 shows the time-frequency diagram of the system response considering a chirp input with a high level of amplitude and the structure in damaged conditions. In Fig. 5a it is observed the appearance of a quadratic harmonic when we compare with Fig. 3b, that is a consequence of the crack, without signiﬁcant alterations in the behavior of other components of the response (ﬁrst and third harmonics look similar to the ones observed in Fig. 3b). With the propagation of the damage (damage II) the quadratic and cubic harmonics grow up and the resonance frequency changes (see in Fig. 5b). So, it is recommended that the initial propagation of the damage has inﬂuence in the quadratic harmonic of the system and when the extension of the damage is more signiﬁcant, all linear and nonlinear components of the response are affected. The damage detection approach applied has to be able to detect the appearance of the nonlinear behavior caused by the damage, without confusing this evolution with the cubic nonlinear behavior caused by the presence of the magnet. Additionally, the problem becomes more complicated when the data variation is assumed in the problem, requiring a strategy to separate these effect as shown in the next section. Prepared using sagej.cls 6 Journal Title XX(X) Figure 1. Flowchart of the damage detection methodology proposed. (a) (b) Figure 2. Experimental apparatus; (a) photos and (b) scheme (dimensions in [mm]). 101102 30 32 34 36 38 40 Amplitude [ms-2/V] Frequency [Hz] (a) 3x 1x0.5 1 1.5 2 2.5 3 3.5 Time [s]020406080100120140Frequency [Hz] -60-50-40-30-20-1001020 dB/Hz (b) Figure 3. Intrinsically nonlinear behavior of the system. (a) Stepped sine test, where 4 ,,andrepresent, respectively, frequency up and frequency down considering low level of input, and frequency up and frequency down considering high level of input. (b) Time-frequency diagram of the system response. Data variation In order to study the performance of the strategy in the presence of data variation or other uncertainties, the distance between the magnet and the steel masses (see in Fig. 2b) wasvaried from 2 mm to 3.5 mm during the tests, repeated on different days to obtain a total of 200 experimental tests for each beam constructed. Prepared using sagej.cls Villani et al. 7 (a) (b) Figure 4. Structural conditions considered; (a) all beams constructed and (b) damaged beams (zoom). 3x 2x 1x 0.5 1 1.5 2 2.5 3 3.5 Time [s]020406080100120140Frequency [Hz] -60-50-40-30-20-1001020 dB/Hz (a) 3x 2x 1x0.5 1 1.5 2 2.5 3 3.5 Time [s]020406080100120140Frequency [Hz] -60-50-40-30-20-10010 dB/Hz (b) Figure 5. Time-frequency diagram of the system response obtained in damaged conditions; (a) damage I and (b) damage II. Figure 6 illustrates the variation of the system response during the tests. The results consider only the ﬁrst mode shape frequency range, as only this range will be examined in the damage detection process. It is clear that it is easy to classify the structural condition when the extension of the damage is even more severe (damage II), but it is not plausible to observe visually the deviation between the reference and damage I conditions when the Frequency Response Function (FRF) and the linear modal parameters are considered. These modal parameters were estimated considering a line-ﬁt method51and the experimental realizations of the FRF. These results show how hard it is to detect the presence of the damage when the data variation is considered in the analysis. In the scenario considered in this work, with the experimental data measured, the model proposed has to be able to detect the presence of the damage in the uncertain ambient state without confusions between the nonlinear behavior caused by the presence of the magnet and that one caused by the presence of the damage. In this sense, the next section presents the main results obtained with the application of the methodology based on stochastic V olterra series. Application of the proposed methodology This section matches the application of the methodology described in the section “The damage detection methodologybased on stochastic V olterra series” to detect damage admitting the experimental setup described in section “Experimental setup”. The main results obtained are shown, and the performance of the linear and nonlinear analysis are compared. Reference model identiﬁcation The ﬁrst step to the utilization of the methodology is the estimation of the stochastic reference model. As mentioned before, the ﬁrst three V olterra kernels were considered in the analysis. The number of Kautz functions and the Kautz parameters related to each V olterra kernel were deﬁned as described in the previous work Villani et al. (2018)36. So, the number of functions used here are J1= 2,J2= 4 andJ3= 6. To obtain the input/output signals used in the kernels estimation, the structure was excited admitting a chirp input signal varying the excitation frequency from 25 to 40 Hz (ﬁrst mode shape region) and considering two levels of amplitude (1 V RMS - linear system behavior and 6 V RMS - nonlinear system behavior). The chirp input is attractive because it excites the linear and nonlinear components of the system response with enough energy, leading to a better estimation of the high-order kernels19,37,52. The output signal considered (velocity vibration signal) is obtained through the integration of the acceleration signal measured by the accelerometer (see Fig. 2). The velocity signals are used because of the previous implementation Prepared using sagej.cls 8 Journal Title XX(X) 10-1100101102 20 25 30 35 40 45 Amplitude [ms-2/V] Frequency [Hz] (a)282930313233343536 0.5 1 1.5 2 2.5 Ressonance frequency [Hz] Damping ratio [%] (b) Figure 6. Variation of the data measured. (a) Frequency Response Function, where the continuous lines represent the mean values and the color zones the curves obtained with 99% of conﬁdence, B - reference, B - damage I and B - damage II. (b) Variation of linear modal parameters, where represents the reference condition, 4damage I and damage II. of the model identiﬁcation procedure considering this type of output. As done before36,37, the kernel estimation is performed in two steps, i.e., the ﬁrst kernel is identiﬁed considering the underline linear behavior of the system, and then, the second and third kernels are identiﬁed considering the nonlinear components of the system response. The deterministic model is identiﬁed several times, considering the 200 experimental realizations obtained from the healthy training beam, to construct the stochastic reference model. Figure 7 presents the veriﬁcation and validation of the model identiﬁed. Figures 7a and 7b shows the stochastic model output with 99% of conﬁdence bands, in the time domain, considering the same chirp input used in the model estimation with a high level of amplitude (6 V RMS), in comparison with experimental data measured. It is observed that the stochastic model can predict the system output considering the data variation. Additionally, Fig. 7c shows the stochastic model output with 99% of conﬁdence bands, in the frequency domain, holding a single sine input with excitation frequency close to the system resonance frequency (35Hz) and a high level of amplitude (6 V RMS), in comparison with experimental data measured. As can be seen, the stochastic model can predict the harmonic components of the response. This characteristic is unusual in the sense that as the damage considered produces variations in the nonlinear components of the response, it is expected that the model could be sensitive to these variations with adequate performance to detect the damage faster than the linear approach. With the stochastic reference model estimated, the damage detection procedure can be applied and the results obtained are shown in the next section. Damage detection performance First of all, it is interesting to observe whether the Mahalanobis distance computed considering the indices obtained in the reference condition ( IlinandInlin) belongs to the chi-square distribution, as proposed in the methodology. Figure 8 shows adjustment between the histograms obtained from the Mahalanobis distance, calculated with the indexes in the reference condition, and the chi-square theoreticaldistribution. Of course, the theoretical approximation is not perfect but satisfactory, considering the limited amount of data available in the analysis. Moreover, this approximation overcomes the problems involving the empirical estimation of distributions based on experimental data, and the Kernel Density Estimator used before, related to the needed amount of experimental data and the choices of the kernel and smoothing parameter47,53. This aspect of the indices distribution allows the application of the hypothesis test proposed based on the probability of the distance calculated belonging to the chi-square distribution. The evolution of the distance calculated from the indices with the propagation of the damage is presented in Fig. 9. It is observed that the linear index is not sensitive to the presence of the damage at the beginning of the propagation, i.e., the damage I condition (see Fig. 9a). The nonlinear index presented more sensitivity to the presence of damage, showing an adequate separation between the damage I condition and the reference condition (see Fig. 9b). Besides, the nonlinear index calculated from the data measured considering the test beam presents some outliers that will be reﬂected, probably, in false positives, depending on the threshold value used. These outliers are a consequence of the high sensitivity of the nonlinear components to structural variations. These results are not unexpected since the test and training beams are nominally the same, but not identical. The use of more than one training beam could improve the methodology performance in a real application. Another interesting aspect of the index used is the increase of the distance values with the propagation of the damage that, in the future, may be well correlated to the severity of the damage. Considering the distances determined, the hypothesis test proposed may be applied. Table 1 shows the results obtained for both linear and nonlinear indices. It is clear that the nonlinear index presented a higher capability to detect the presence of the damage. The false alarms (false detection) are also higher considering the nonlinear index. However, this value can be reduced without loss of performance to detect the damage with the decrease of the test sensitivity as perfomed, e.g., by Avenda ˜no-Valencia and Fassois (2017)31. Prepared using sagej.cls Villani et al. 9 -0.2-0.15-0.1-0.0500.050.10.150.2 2 2.2 2.4 2.6 2.8 3 3.2 3.4 Velocity [m/s] Time [s] (a)-0.2-0.15-0.1-0.0500.050.10.150.2 2.7 2.72 2.74 2.76 2.78 2.8 Velocity [m/s] Time [s] (b)10-1410-1210-1010-810-610-410-2100 20 40 60 80 100 120 Amplitude [(m/s)2/Hz] Frequency [Hz] (c) Figure 7. Veriﬁcation and validation of the reference model. (a) Model response in the time domain considering a chirp excitation. (b) Model response in the time domain considering a chirp excitation (zoom). (c) Model response in the frequency domain considering a sine excitation. B represents the 99% model response conﬁdence bands, represents the model response mean andrepresents ﬁve realizations of the experimental data. 00.050.10.150.20.25 0 10 20 30 40 50 Probability Density Function D2 lin (a)00.020.040.060.080.1 0 20 40 60 80 100 Probability Density Function D2 nlin (b) Figure 8. Comparison between the theoretical and experimental distributions. (a) Distance calculated based on the linear index in the reference condition. (b) Distance calculated based on the nonlinear index in the reference condition. The histogram represents the experimental data and the black line represents the theoretical chi-square distribution. 100101102103104105 Traning Test Damage I Damage II D2 lin Structure Condition (a)100101102103104105 Traning Test Damage I Damage II D2 nlin Structure Condition (b) Figure 9. Evolution of the Mahalanobis distance calculated with the progression of the damage. (a) Distance calculated based on the linear index. (b) Distance calculated based on the nonlinear index. represents the reference training beam, represents the reference test beam, 4represents the damage I condition beam and represents the damage II condition beam. The linear index presented problems to detect the damage I condition, even with the variation of the test sensitivity. In real applications, the value of the test sensitivity ( ) has to be determined depending on the level of security/conﬁdence required and on the previous knowledge about the monitoredstructure. In some applications, this value can also be optimized to present a better performance in the damage detection process, based on measured data obtained from damaged structures31. Prepared using sagej.cls 10 Journal Title XX(X) Table 1. Results obtained through the application of the hypothesis test. Hypothesis test sensitivity Percentage of false detection [%] Percentage of true detection [%] ( ) Training Beam Test Beam Damage I Damage II Linear Index1020.5 15.5 34.5 100 1040 0 5 100 1060 0 0 100 10120 0 0 75 Nonlinear Index1024.5 12.5 100 100 1041 7 100 100 1060 4.5 100 100 10120 1.0 96.5 100 However, a better way to analyze the damage detection capability of the methodology is computing the Receiver Operating Characteristics (ROC) curve. This curve relates the false alarm ratio with the true detection ratio obtained applying the hypothesis test for different values of sensitivity ( ). As a result, the closer to the point (0,1) is the curve, and better is the performance of the index because it presents a higher probability of detecting the damage with a low level of false alarms. This curve is used here to study the discrepancies between the performance of the linear and nonlinear approaches proposed. Figure 10 shows the curves obtained after the application of the hypothesis test using different threshold values, considering both indices (linear and nonlinear) and all structural conditions studied (training beam, test beam, the damage I and damage II) . In the ﬁgure, it can be seen that the nonlinear index presents a better performance since the red curve is close to the point (0,1). This result was expected as the linear components of the response are not sensitive to the initial propagation of the damage (see Fig. 9 and Tab. 1). Moreover, the performance of the nonlinear index (almost perfect, but this is admittedly a consequence of a ﬁnite data set), with the curve very close to the point (0,1), was achieved in an experimental laboratory application. Even with the consideration of the data variation simulated during the experimental tests, these results do not reﬂect a real-structure application. It is supposed in a real-world application a higher number of confounding effects, different types of damage occurring coincidentally and other aspects that may decrease the methodology performance. However, the capability of the features related to the nonlinear kernels to detect the presence of the damage before the ones reported to the linear kernel and the ability of the methodology to do not confuse the intrinsically nonlinear behavior of the system with the presence of damage, showed as the main result of this work, must be preserved. Final remarks The problem of damage detection in an intrinsically nonlinear system, regarding the data variation associated to uncertainties, subjected to the presence of damage that causes a nonlinear characteristic in the system response was investigated in this paper. In this sense, an initially nonlinear beam was analyzed, and the data variation was emulated by the random variation imposed in the experimental setup (variation of the distance between the magnet and the beam). The damage studied was a breathing crack that0.50.60.70.80.91.0 0 0.2 0.4 0.6 0.8 1.0 True Detection Ratio False Alarm Ratio Figure 10. Receiver Operating Characteristics (ROC) curve. represents the linear index and represents the nonlinear index. effects the system to exhibit a nonlinear operation with a distinct character of the initial one. In this condition, the methodology has to be adequate to distinguish the intrinsically nonlinear operation and the data variation to the presence of the damage. A method based on a stochastic version of the V olterra series, with the use of the random kernels coefﬁcients and contributions as damage detection features, combined with a novelty detection technique, was applied to solve the problem. Unlike what has been shown in the previous works published36,38, the kernels coefﬁcients, and the kernels contributions approach were applied together, considering a unique index to monitor the structural condition, aiming to augment the robustness of the method. Moreover, a theoretical distribution was introduced to the Mahalanobis distance computed in the reference condition to reduce the possible problems related to the use of the kernel density estimator previously assumed. Additionally, for the ﬁrst time, the methodology was examined through an experimental application considering an intrinsically nonlinear beam subordinated to the presence of damage that produces nonlinear aspects to the system response, all this considering data variation that reﬂects changes in the linear and nonlinear components of the response. The results obtained revealed that the monitoring of the nonlinear components of the system response, denoted by the high-order kernels coefﬁcients and contributions considered as damage detection features, is a helpful method to be implemented in damage detection problems when the Prepared using sagej.cls Villani et al. 11 nonlinear system response is present. Again, the nonlinear metric conﬁrmed to be more sensitive to the appearance of the damage and showed better performance considering the data variation compared with the monitoring of linear components. Finally, the use of the stochastic reference model, combined with the novelty detection technique and the hypothesis test, exhibited satisfactory performance to overcome the problem related to the data measured variation, providing the metric for detecting the presence of the damage with probabilistic conﬁdence even in an uncertain ambient. Based on the results achieved, although the application was performed considering an experimental laboratory setup, the authors believe that in real-world applications the nonlinear metric will also demonstrate a higher performance than the linear one for this kind of problem. Declaration of conﬂicting interests The authors declared no potential conﬂicts of interest with respect to the research, authorship, and/or publication of this article. Funding The authors are thankful for the ﬁnancial support provided from of S ˜ao Paulo Research Foundation (FAPESP), grant numbers 2012/09135-3, 2015/25676-2, 2017/24977-4 and 2017/15512-8, from Carlos Chagas Filho Research Foundation of Rio de Janeiro State (FAPERJ) under grants E-26/010.002.178/2015 and E-26/010.000.805/2018, Coordenac ¸ ˜ao de Aperfeic ¸oamento de Pessoal de N ´ıvel Superior - Brasil (CAPES) - Finance Code 001, and Brazilian National Council for Scientiﬁc and Technological Development (CNPq) grant number 307520/2016-1.