the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Experimental validation of a shortterm damping estimation method for wind turbines in nonstationary operating conditions
Philippe Jacques Couturier
Lars Morten Sørensen
Jon Juel Thomsen
Modal properties and especially damping of operational wind turbines can vary over short time periods as a consequence of environmental and operational variability. This study seeks to experimentally test and validate a recently proposed method for shortterm damping and natural frequency estimation of structures under the influence of varying environmental and operational conditions from measured vibration responses. The method is based on Gaussian process timedependent autoregressive moving average (GPTARMA) modelling and is tested via two applications: a laboratory threestorey shear frame structure with controllable, timevarying damping and a flutter test of a fullscale 7 MW wind turbine prototype, in which two edgewise modes become unstable. Damping estimates for the shear frame compare well with estimates obtained with stochastic subspace identification (SSI) and standard impact hammer tests. The efficacy of the GPTARMA approach for shortterm damping estimation is illustrated through comparison to shortterm SSI estimates. For the fullscale flutter test, GPTARMA model residuals imply that the model cannot be expected to be entirely accurate. However, the damping estimates are physically meaningful and compare well with a previous study. The study shows that the GPTARMA approach is an effective method for shortterm damping estimation from vibration response measurements, given that there are enough training data and that there is a representative model structure.
 Article
(5872 KB)  Fulltext XML
 BibTeX
 EndNote
A novel operational modal analysis (OMA) method (Ebbehøj et al., 2023) for shortterm modal damping estimation for structures under the influence of varying environmental and operational conditions (EOCs), such as wind turbines in operation, is tested with a controlled laboratory experiment and a wind turbine flutter test.
The dynamic properties of wind turbines (i.e. natural frequencies, damping ratios, and mode shapes) can be sensitive to changing EOCs (AvendañoValencia et al., 2017; Bogoevska et al., 2017). Aeroelasticity, active control, material properties, and nonlinear damping mechanisms (e.g. friction) are examples of phenomena and factors which can cause EOC variability (M. O. Hansen et al., 2006; Wang et al., 2022; Chen and Duffour, 2018). EOC variability can act on short and long timescales relative to the fundamental frequency of the given structure. For example, changing temperatures of wind turbine towers affecting their natural frequencies work on the order of hours and days (Hu et al., 2015), which in this context are considered longterm effects. By contrast, complex aero–servo–elastic interactions of multimegawatt wind turbines can cause shortterm variability of damping in particular over a few minutes due to dependencies on e.g. rotor speed, blade pitch angle, and wind speed (M. O. Hansen et al., 2006; M. H. Hansen et al., 2006).
Estimating aeroelastic (or operational) damping accurately is essential for improving the design of multimegawatt wind turbines, as it is a crucial design parameter for modelling fatigue and aeroelastic instabilities, i.e. stallinduced vibrations (SIVs) and vortexinduced vibrations (VIVs) (Veers et al., 2023). Improved aeroelastic damping estimation may therefore enable designs associated with less risk or less material usage. However, “a precise evaluation of aeroelastic damping remains an elusive goal in some operating conditions”, as stated in Grand challenges in the design, manufacture, and operation of future wind turbine systems (Veers et al., 2023, p. 1090). The combination of nonstationary EOCs and EOCsensitive damping complicates the task of obtaining precise aeroelastic damping estimates for wind turbines.
OMA covers a broad class of outputonly system identification methods for estimating modal parameters for structures in operating conditions where the input (i.e. forcing or excitation) is not measured. Standard OMA techniques, for instance, the covariancedriven stochastic subspace identification (COVSSI) (Peeters and Roeck, 2001; Brincker and Ventura, 2015), typically assume the system is linear and time invariant (LTI) with input resembling stationary white noise and requiring long time measurements. Brincker and Ventura (2015) suggest a minimum measurement time of $\frac{\mathrm{10}}{f\mathit{\zeta}}$, where fζ is the lowest natural frequency–damping ratio multiple of a given structure. This translates to a measurement time requirement of approximately 28 min for a mode with a natural frequency of 0.2 Hz and 3 % damping ratio, which is in the range of a multimegawatt wind turbine mode. Therefore, the capabilities to track shortterm variations in damping are limited for these methods.
When identifying modal damping from outputonly measurements, the effect of input on the measured output must be accounted for since both input and damping govern the nearresonant vibration amplitudes. This can be done by either eliminating the effect (i.e. averaging it out) or modelling it (Au, 2017). Standard OMA methods rely on the former approach. For instance, in COVSSI covariance, Toeplitz matrices of the output measurements are computed and used as equivalent free response approximates of an assumed LTI system (Peeters and Roeck, 2001). However, a considerable amount of data is required (minimum measurement time of $\frac{\mathrm{10}}{f\mathit{\zeta}}$) for the covariance Toeplitz matrices to be adequately estimated. Consequently, this approach involves a tradeoff between temporal resolution and estimate accuracy, which can be a limiting factor in the context of shortterm variability.
Nonstationary autoregressive moving average (ARMA) time series models offer an avenue for accounting for nonstationary input and timevarying system characteristics, including modal parameters. ARMA models closely resemble the mathematical structure of discretetime equations of motion, where the autoregressive (AR) part plays the role of the lefthand (homogeneous) side, and the AR model coefficients carry information on the modal parameters. Similarly, the moving average (MA) part resembles the input, thus filtering out the effect of the excitation from the modal parameters. In timedependent autoregressive moving average (TARMA) models, the model coefficients and corresponding modal parameters can vary in time. One type is the smoothness priors TARMA (SPTARMA) model, where the model coefficients are modelled as autocorrelated stochastic (random) variables, for which evolution in time is constrained by smoothness priors (Poulimenos and Fassois, 2006; Spiridonakos et al., 2010; Kitagawa and Gersch, 1996). The model coefficients of the SPTARMA model are estimated locally in time, e.g. with a Kalman filter. The SPTARMA models may be capable of tracking general nonstationary signals but have limited capabilities in tracking abrupt or shortterm changes (Poulimenos and Fassois, 2006). Functional series TARMA (FSTARMA) models constitute TARMAtype models whose model coefficients are represented by predefined basis functions, allowing AR and MA coefficients to evolve deterministically in time. FSTARMA models are fitted to data by estimating global projection coefficients for each basis function, which minimizes the model prediction errors. If the basis functions can capture the timevarying nature of the system, FSTARMA models can track abrupt and shortterm changes in the system. However, modal parameter trajectories in time for complex systems under EOC variability might not lend themselves to be represented by a reasonable number of timedependent basis functions, resulting in an intractable number of parameters to be estimated.
One approach to capture the effects of changing EOCs on vibrating structures is to embed measured environmental and operational variables (EOVs) into the model. Various approaches have been proposed for this. Multimegawatt wind turbines pose a particular challenge due to the intricate aero–servo–elastic interactions. Bogoevska et al. (2017) expand SPTARMA model residuals with a polynomial chaos expansion (PCE) to account for longterm EOC variability in order to improve the accuracy of structural health monitoring (SHM) of wind turbines. AvendañoValencia et al. (2017) introduce a linearparametervarying AR (LPVAR) model to capture the shortterm dynamics, and Gaussian process (GP) regression is used to account for longterm variability associated with changing wind speed in terms of 10 min averages, which is extended from a singleoutput (univariate) model to a multipleoutput (multivariate) model by AvendañoValencia et al. (2020) and AvendañoValencia and Chatzi (2020), which e.g. enables mode shape identification.
The present work concerns shortterm damping (and natural frequency) estimation based on outputonly measurements for structures influenced by shortterm varying EOCs, where short term is of the order of seconds. The methods mentioned above are based on models conditioned on (e.g. 10 min) statistics, which limits their ability to track shortterm changes. The nonstationary and EOVdependent GPTARMA model (introduced in Ebbehøj et al., 2023) is therefore conditioned on EOV time series. The GPTARMA model combines an FSTARMA model where the basis functions may depend on multiple EOV time series with Gaussian processes by modelling the projection coefficients for the basis functions as Gaussian rather than deterministic variables to allow for better representation of unaccounted disturbances. The capabilities of the GPTARMA model to track EOC variability are limited by how well the basis functions capture the nonstationarities of the response (e.g. slow or fast variations) and fundamentally by the measurement sampling rate. While the GPTARMA model may be nonlinear with respect to EOVs, it is linear with respect to the response it is modelling, i.e. representing an equation of motion that is linear in the dependent variables. Consequently, the model cannot capture strongly nonlinear system properties. However, weakly nonlinear effects on the effective natural frequencies and damping ratios may be approximated if these nonlinear effects correlate with operational states represented by the EOVdependent basis functions.
Verification and validation are integral parts of establishing any new method or model in structural dynamics, and it is essential for outputonly damping estimation methods due to the latent and elusive nature of damping. This work contributes, in particular, to experimental validation of the GPTARMA approach for shortterm damping estimation, suitable for application to wind turbines. The method is validated using vibration measurements from two distinctly different experimental setups: a laboratory shear frame with abruptly changing damping realized with electromagnetic dampers and a fullscale 7 MW wind turbine prototype deliberately driven to flutterlike instabilities (measurements published by Volk et al., 2020).
The paper is structured as follows: Sect. 2 summarizes the details necessary for using the GPTARMA model for shortterm modal parameter estimation from Ebbehøj et al. (2023), including the GPTARMA model definition and estimation, a model structure identification scheme, a model validation procedure, and downstream analysis of an estimated GPTARMA model. Section 3 presents the laboratory shear frame test setup; experimental procedures; and related analysis, results, and discussion. In Sect. 4, analysis of the fullscale multimegawatt wind turbine instability measurements is performed, and the results are presented and discussed.
This section summarizes the procedures for estimating shortterm damping ratios (and natural frequencies) from outputonly measurements using a GPTARMA model, which is introduced and presented in detail in Ebbehøj et al. (2023). Necessary details for using the method are given, including the GPTARMA model definition, procedures for estimating model parameters, the selection of an appropriate model structure, the extraction of modal parameters, and their uncertainties. The entire procedure is summarized in Table 2.
2.1 The GPTARMA model
The GPTARMA model introduced in Ebbehøj et al. (2023) can be used to model a nonstationary, discretetime (displacement/velocity/acceleration) response y_{t}∈ℝ influenced by EOCs, which can be described using m EOVs ${\mathit{\xi}}_{t}\in {\mathbb{R}}^{\mathrm{1}\times m}$, where subscript t denotes the time index, and y_{t} and ξ_{t} are defined for $t=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}N$. The GPTARMA model is closely related to an FSTARMA model (Poulimenos and Fassois, 2006; Spiridonakos and Fassois, 2014), for which the ARMA coefficients evolve in time as trajectories spanned by EOVdependent basis functions. The GPTARMA model extends FSTARMA by modelling the basis function coefficients as Gaussian variables:
where e_{t} is a normally and independently distributed (NID) zeromean innovation process with variance ${\mathit{\sigma}}_{e,t}^{\mathrm{2}}$, which may be time varying; a_{i}(ξ_{t}) and c_{i}(ξ_{t}) are the ith EOVdependent AR and MA coefficients; and n_{a}(n_{c}) denotes the AR (MA) model order. The AR and MA coefficients a_{i}(ξ_{t}) and c_{i}(ξ_{t}) are linear combinations of basis functions g_{j}(ξ_{t}) and h_{j}(ξ_{t}) and the associated Gaussian projection coefficients ${u}_{{\mathrm{a}}_{i,j}}$ and ${u}_{{\mathrm{c}}_{i,j}}$. The complete set of AR (MA) basis functions for the full time series $t=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}N$ constitutes a functional subspace defined as
where ${\mathit{g}}_{j}\left(\mathit{\xi}\right)=\left[{g}_{j}\right({\mathit{\xi}}_{\mathrm{1}}),\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{g}_{j}({\mathit{\xi}}_{N}){]}^{T}\in {\mathbb{R}}^{N\times \mathrm{1}}$ and ${\mathit{h}}_{j}\left(\mathit{\xi}\right)=\left[{h}_{j}\right({\mathit{\xi}}_{\mathrm{1}}),\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{h}_{j}({\mathit{\xi}}_{N}){]}^{T}\in {\mathbb{R}}^{N\times \mathrm{1}}$ are the jth basis functions for time index $t=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}N$, and $\mathit{\xi}\in {\mathbb{R}}^{N\times m}$ contains the EOVs at the corresponding time indices. The basis functions g_{j} and h_{j} may be members of different orthogonal function families (e.g. Legendre polynomials or trigonometric functions) and depend on different EOVs (e.g. rotor speed, wind speed, and pitch angle). A basis function type consists of a function family and an EOV. The functional subspace ℱ_{AR} (ℱ_{MA}) may consist of k_{a} (k_{c}) basis function types, each with a basis order of ${p}_{{\mathrm{a}}_{k}}$ (${p}_{{\mathrm{c}}_{k}}$).
For example, ℱ_{AR} could consist of a constantvalued bias vector g_{1}; first and secondorder Legendre polynomials in wind speed v, g_{2}(v) and g_{3}(v); and a firstorder Legendre polynomial in rotor speed Ω, g_{4}(Ω). This would correspond to k_{a}=2 basis function types (neglecting the trivial constantvalued bias type) with basis orders of ${p}_{{\mathrm{a}}_{\mathrm{1}}}=\mathrm{2}$ and ${p}_{{\mathrm{a}}_{\mathrm{2}}}=\mathrm{1}$ for the Legendre polynomials in wind speed and rotor speed, respectively. The orthogonality among basis functions in the functional subspaces should be ensured, using e.g. modified Gram–Schmidt orthogonalization (Stewart, 2013). In cases where EOVs contain highfrequency (i.e. much higher than the frequency of the highest mode) scatter, it can be beneficial to zerophase lowpass filter (e.g. using MATLAB's filtfilt
function) the EOVs to prevent the highfrequency scatter from propagating to the modal parameters (Ebbehøj et al., 2023).
ARMA models are closely linked with discretetime equations of motion (EOMs). The AR part resembles the lefthand side of discretetime EOMs, which means the AR coefficients carry the physical characteristics of the system it models, i.e. natural frequencies and damping ratios. The MA part resembles the righthand side of discretetime EOMs as it can capture the effect of stochastic excitation on the measured response y_{t}. Consequently, ℱ_{AR} should be composed such that it can represent the timevarying nature of the natural frequencies and damping ratios, and ℱ_{MA} should be selected to account for nonstationary stochastic excitation; examples of this are given in Sects. 3.2 and 4.1.
Equations (1)–(3) can be expressed as follows in the more compact regression form:
where the regression vector ${\mathbf{\Phi}}_{t}={\mathbf{\Phi}}_{t}\left({\mathit{\xi}}_{t}\right)\in {\mathbb{R}}^{({n}_{\mathrm{a}}{p}_{\mathrm{a}}+{n}_{\mathrm{c}}{p}_{\mathrm{c}})\times \mathrm{1}}$ contains the regressed response measurement y_{t−i} and innovation e_{t−i} multiplied by appropriate basis function values, and the Gaussian projection coefficients are collected in $\mathit{\theta}\in {\mathbb{R}}^{({n}_{\mathrm{a}}{p}_{\mathrm{a}}+{n}_{\mathrm{c}}{p}_{\mathrm{c}})\times \mathrm{1}}$. The GPTARMA model for the full time series can be written in regression form by employing the stacked response and innovation vectors $\mathit{y}={\left[{y}_{\mathrm{1}+{n}_{m}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{y}_{N}\right]}^{T}\in {\mathbb{R}}^{(N{n}_{m})\times \mathrm{1}}$ and $\mathit{e}={\left[{e}_{\mathrm{1}+{n}_{m}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{e}_{N}\right]}^{T}\in {\mathbb{R}}^{(N{n}_{m})\times \mathrm{1}}$:
where ${n}_{m}=\mathrm{max}\left(\right[{n}_{a},{n}_{c}\left]\right)$ denotes the maximum model order, and $\mathbf{\Phi}\in {\mathbb{R}}^{({n}_{\mathrm{a}}{p}_{\mathrm{a}}+{n}_{\mathrm{c}}{p}_{\mathrm{c}})\times (N{n}_{m})}$ is the regression matrix.
For a given data set $\mathcal{D}=\mathit{\{}\mathit{y},\mathit{\xi}\mathit{\}}$, the GPTARMA model is fully characterized by $\mathcal{M}=\mathit{\{}\mathcal{S},\mathcal{P}\mathit{\}}$, where $\mathcal{S}=\mathit{\{}{\mathcal{F}}_{\mathrm{AR}},{\mathcal{F}}_{\mathrm{MA}},{n}_{\mathrm{a}},{n}_{\mathrm{c}}\mathit{\}}$ contains the model structure, and $\mathcal{P}=\mathit{\{}{\mathit{\mu}}_{\mathit{\theta}},{\mathbf{\Sigma}}_{\mathit{\theta}},{\mathbf{\Sigma}}_{e}\mathit{\}}$ contains the timevarying innovation variance ${\mathbf{\Sigma}}_{e}=\mathrm{diag}\left(\right[{\mathit{\sigma}}_{e,\mathrm{1}+{n}_{m}}^{\mathrm{2}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{\mathit{\sigma}}_{e,N}^{\mathrm{2}}\left]\right)\in {\mathbb{R}}^{(N{n}_{m})\times (N{n}_{m})}$ and the hyperparameters for the Gaussian projection coefficients, i.e. the mean ${\mathit{\mu}}_{\mathit{\theta}}\in {\mathbb{R}}^{({n}_{\mathrm{a}}{p}_{\mathrm{a}}+{n}_{\mathrm{c}}{p}_{\mathrm{c}})\times \mathrm{1}}$ and covariance ${\mathbf{\Sigma}}_{\mathit{\theta}}\in {\mathbb{R}}^{({n}_{\mathrm{a}}{p}_{\mathrm{a}}+{n}_{\mathrm{c}}{p}_{\mathrm{c}})\times ({n}_{\mathrm{a}}{p}_{\mathrm{a}}+{n}_{\mathrm{c}}{p}_{\mathrm{c}})}$.
2.2 Model parameter estimation
A procedure for maximum likelihood (ML) estimation of the hyperparameters and innovation variance in 𝒫 (via the expectation maximization algorithm) presented in Ebbehøj et al. (2023) (see also AvendañoValencia et al., 2017) is summarized in this section.
With the model parameter θ and response y in Eq. (6) modelled as random variables, the probability density function (PDF) of their joint distribution governs the probability of the two (model parameters and response) occurring simultaneously:
where $p\left(\mathit{\theta}\right\mathit{y},\mathbf{\Phi},\mathcal{P})$ denotes the posterior distribution of the model parameters and $p\left(\mathit{y}\right\mathbf{\Phi},\mathcal{P})$ the marginal likelihood of the response. For the GPTARMA model in Eq. (6), these distributions are Gaussian and can be specified as (Ebbehøj et al., 2023) (see also Murphy, 2023, Chap. 2, and AvendañoValencia et al., 2017)
where
where 𝔼[X] is the expected value of X, ${\mathbf{\Sigma}}_{\mathit{\epsilon}}\in {\mathbb{R}}^{(N{n}_{m})\times (N{n}_{m})}$ is the covariance matrix of the prior (i.e. before observing the actual response) onestepahead prediction errors, and the innovation variance for each time index is collected in ${\mathbf{\Sigma}}_{e}=\mathrm{diag}\left({\mathit{\sigma}}_{\mathrm{1}+nm}^{\mathrm{2}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{\mathit{\sigma}}_{N}^{\mathrm{2}}\right)\in {\mathbb{R}}^{(N{n}_{m})\times (N{n}_{m})}$. The mean and covariance of the posterior parameter distribution, $\widehat{\mathit{\theta}}\in {\mathbb{R}}^{({n}_{\mathrm{a}}{p}_{\mathrm{a}}+{n}_{\mathrm{c}}{p}_{\mathrm{c}})\times \mathrm{1}}$ and ${\widehat{\mathbf{P}}}_{\mathit{\theta}}\in {\mathbb{R}}^{({n}_{\mathrm{a}}{p}_{\mathrm{a}}+{n}_{\mathrm{c}}{p}_{\mathrm{c}})\times ({n}_{\mathrm{a}}{p}_{\mathrm{a}}+{n}_{\mathrm{c}}{p}_{\mathrm{c}})}$, are also referred to as the maximum a posteriori (MAP) estimates of μ_{θ} and Σ_{θ}, respectively.
The marginal likelihood of the response in Eq. (9) also serves as the marginal likelihood of the hyperparameters $\mathcal{L}\left(\mathcal{P}\right\mathit{y},\mathbf{\Phi})$, meaning that the hyperparameters in 𝒫 can be estimated by maximizing the marginal hyperparameter likelihood. This forms a ML optimization problem, which can be formulated in terms of the loglikelihood as (AvendañoValencia et al., 2017; Murphy, 2023; Rasmussen and Williams, 2006):
where the vector $\mathit{\epsilon}=[{\mathit{\epsilon}}_{\mathrm{1}+nm},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{\mathit{\epsilon}}_{N}{]}^{T}\in {\mathbb{R}}^{(N{n}_{m})\times \mathrm{1}}$ contains the prior onestepahead prediction errors computed at time index t as
Solutions to the ML optimization can be approximated using the generalpurpose expectation maximization (EM) algorithm, which constitutes an expectation step (E step) and a maximization step (M step), which are iterated until convergence. During the E step, the posterior (expected) mean and covariance of the model parameters, $\widehat{\mathit{\theta}}$ and ${\widehat{\mathbf{P}}}_{\mathit{\theta}}$, are evaluated using Eqs. (10a)–(10d) based on the hyperparameters from the previous EM iteration. In the M step, the hyperparameters are updated using the model parameter posterior mean and covariance from the E step by the explicit update equations as follows (AvendañoValencia et al., 2017):
where superscript (i) denotes the designated variable from the ith EM iteration. The timevarying innovation variance ${\widehat{\mathbf{\Sigma}}}_{e}\in {\mathbb{R}}^{(N{n}_{m})\times (N{n}_{m})}$ can be estimated using a 2K+1 sample moving window:
where ${e}_{\mathit{\tau}}={y}_{\mathit{\tau}}({\mathbf{\Phi}}^{(i\mathrm{1})}{)}^{T}{\widehat{\mathit{\theta}}}^{(i\mathrm{1})}$, and ${\widehat{\mathit{\mu}}}_{e,t}$ is the innovation mean for the window with time index t. Although the latter is assumed to be zero, it is included in Eq. (14) because it has been observed to improve the accuracy of damping ratio estimates in cases where the response measurements are influenced by stochastic excitation with a timevarying mean value. The timeinvariant innovation variance can likewise be estimated using Eq. (14) with $K=N/\mathrm{2}$.
The EM algorithm has been shown to converge to a local likelihood maximum, not necessarily global (e.g. Bishop, 2006). Consequently, accurate initial values of the hyperparameters 𝒫^{(0)} are required. This can be obtained by estimating the model parameters for the corresponding FSTARMA model, where ${u}_{{\mathrm{a}}_{i,j}}$ and ${u}_{{\mathrm{c}}_{i,j}}$ are scalars rather than Gaussian variables, and using these parameter estimates as the initial hyperparameter values 𝒫^{(0)}. The corresponding FSTARMA model parameters can be estimated using e.g. the twostage least squares (2SLS) or the twostage weighted least squares (2SWLS) method (Spiridonakos and Fassois, 2014; Poulimenos and Fassois, 2006). Convergence of the EM algorithm is implied by consistent and small changes in $\mathrm{ln}\mathcal{L}\left(\mathcal{P}\right\mathit{y},\mathbf{\Phi})$ and in the hyperparameters 𝒫 over EM iterations.
2.3 Model validation
Once a GPTARMA model is estimated, it is important to validate that it adequately represents the observations (Poulimenos and Fassois, 2006; Spiridonakos and Fassois, 2014; Madsen, 2007). Model validation in this present work (as in Ebbehøj et al., 2023) consists of residual analysis and crossvalidation. Residual analysis tests whether the model residual e_{t} (i.e. the innovations) satisfies the NID assumption, i.e. whether it resembles white noise. A standard whiteness test checks whether the autocorrelation function (ACF) of the residuals resembles that of white noise, i.e. significantly uncorrelated at time lags τ>0. However, this is not generally applicable to the case of nonstationary residuals, as these would be correlated through the timevarying innovation variance. An approach to partially circumvent this issue is to standardize the residuals using the estimated timevarying innovation variance ${\widehat{\mathit{\sigma}}}_{e,t}^{\mathrm{2}}$ (Fouskitakis and Fassois, 2002):
Given that e_{t} is a zeromean white noise sequence with a timevarying variance that is adequately approximated by ${\widehat{\mathit{\sigma}}}_{e,t}^{\mathrm{2}}$, the standardized residual z_{t} is stationary. Using the standardized residual z_{t} for the ACF test renders the test sensitive to the accuracy of the timevarying innovation variance estimate. It should thus be supplemented by a test that is insensitive to such an estimate.
An alternative whiteness test, fully applicable to the nonstationary case, is a simple sign test. Sign changes in a zeromean white noise sequence are expected to occur (on average) every other time step and are governed by the binomial distribution, which is approximated by the Gaussian distribution for large N_{seq} (Madsen, 2007):
where N_{seq} is the number of samples in the residual sequence, and the mean and variance are given by ${\mathit{\mu}}_{\mathrm{wn}}=({N}_{\mathrm{seq}}\mathrm{1})/\mathrm{2}$ and ${\mathit{\sigma}}_{\mathrm{wn}}^{\mathrm{2}}=({N}_{\mathrm{seq}}\mathrm{1})/\mathrm{4}$, respectively. Whether a residual sequence $\mathit{\{}{e}_{t}:t=\mathrm{1}+{n}_{m},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}N\mathit{\}}$ resembles white noise can thus be tested by checking whether the number of sign changes adheres to Eq. (16) under some significance level.
Crossvalidation is performed by splitting a data set of a single recording in a training and a test set containing 75 % and 25 % of the data points, respectively. The model is solely estimated using the training set and is subsequently tested in terms of residual tests and whether the orders of magnitude of prediction errors are the same for the training and test set. If the prediction errors of the training set are much smaller than those of the test set, it suggests the model is overfitted; i.e. the model excessively represents the measured realization of the stochastic response rather than the underlying system. However, it is only applicable to the nonstationary case if the two sets have comparable characteristics.
2.4 Model structure identification
In this section a procedure for identifying a suitable model structure, i.e. $\mathcal{S}=\mathit{\{}{\mathcal{F}}_{\mathrm{AR}},{\mathcal{F}}_{\mathrm{MA}},{n}_{\mathrm{a}},{n}_{\mathrm{c}}\mathit{\}}$, is summarized (see Ebbehøj et al., 2023, for details). A suitable model structure is sufficiently complex to capture the underlying system producing the response measurements while avoiding overfitting. To identify a suitable model structure, a range of candidate models with different model structures is estimated and compared in terms of predictive performance and capability of capturing the vibration modes of interest.
To compare the predictive performance (i.e. the prior onestepahead prediction errors) of the candidate models, the residual sum of squares normalized by the series sum of squares (RSS/SSS) and the Bayesian information criteria (BIC) are used. The RSS/SSS is given by
where ε is the prior prediction error defined in Eq. (12). The BIC is given by
where d is the number of independently adjusted parameters. RSS/SSS and BIC measure the prediction errors, but BIC also penalizes model complexity, i.e. the number of model parameters. The required model complexity for capturing the vibration modes of interest is identified using frequency stabilization diagrams, as is commonly used for time domain OMA methods (Brincker and Ventura, 2015; Peeters and Roeck, 2001; AvendañoValencia et al., 2017). Typically, the predictive performance converges at lower model orders than the model orders required to capture the modes of interest. This means the convergence of the predictive performance is typically a necessary but insufficient condition.
A simple backward regression scheme is employed to identify a suitable model structure, i.e. starting with high model orders, ${n}_{\mathrm{a}}^{*}$ and ${n}_{\mathrm{c}}^{*}$, and functional subspaces of high dimensionalities, ${\mathcal{F}}_{\mathrm{AR}}^{*}$ and ${\mathcal{F}}_{\mathrm{MA}}^{*}$, ensuring sufficient model complexity to capture the nonstationary response. The complexity of the initial model is then reduced in four stages: first, models with lower AR orders are estimated, and the most suitable AR order is identified. Then, the optimal MA order (given the optimal AR order) is identified in a similar fashion. The optimal basis orders for each AR and MA basis function type are identified in a similar manner. These stages are repeated for different combinations of AR and MA basis function types. The procedure is summarized in Table 1.
2.5 Estimating modal parameters and their uncertainties
In this section the necessary results for estimating “frozen” modal parameters (excluding mode shapes) from an estimated GPTARMA model and approximating the corresponding uncertainties are summarized. The frozen properties of a timevarying system represent the LTI properties at frozen time t^{′}; i.e. the timevarying system is represented by an infinite sequence of LTI systems (Poulimenos and Fassois, 2006). See Ebbehøj et al. (2023) for more details and also Poulimenos and Fassois (2006), AvendañoValencia and Fassois (2014), Spiridonakos et al. (2010), AvendañoValencia et al. (2017), and Yang and Lam (2019) on the computation of frozen modal parameters and their uncertainty.
The frozen modal parameters can be computed from the timevarying AR coefficients since an equivalent discretetime statespace model with system matrix ${\mathbf{L}}_{t{}^{\prime}}$ can be formulated from the AR coefficients $\mathit{\left\{}{a}_{i}\right({\mathit{\xi}}_{t{}^{\prime}}):i=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{n}_{\mathrm{a}}\mathit{\}}$ at each time step t^{′}. The first step is to compute the discretetime eigenvalues of the eigenvalue problem:
where t^{′} is the frozentime index, and ${\mathit{v}}_{i}={\mathit{v}}_{i,t{}^{\prime}}$ and ${\mathit{\rho}}_{i}={\mathit{\rho}}_{i,t{}^{\prime}}$ are the ith right eigenvector and eigenvalue of ${\mathbf{L}}_{t{}^{\prime}}$:
where I is the identity matrix with dimensions as indicated by the subscript.
The ith natural frequency ${f}_{i,t{}^{\prime}}$ (in Hz) and damping ratio ${\mathit{\zeta}}_{i,t{}^{\prime}}$ for frozentime index t^{′} can be computed by (Yang and Lam, 2019; Poulimenos and Fassois, 2006; Spiridonakos et al., 2010; AvendañoValencia et al., 2017)
is the equivalent ith continuoustime eigenvalue.
The uncertainty in the modal parameter estimates can be approximated by propagating the estimated AR coefficient uncertainty (quantified by ${\widehat{\mathbf{P}}}_{\mathit{\theta}}$) through the steps required to compute modal parameters from AR coefficients. These steps constitute obtaining discretetime eigenvalues in Eq. (19), transforming the discretetime eigenvalues to continuoustime eigenvalues in Eq. (21b), and computing the modal parameters from the continuoustime eigenvalues in Eq. (21a). The uncertainty propagation can be achieved by either Monte Carlo simulation or analytically using the firstorder delta method. The former is straightforward to implement but is computationally costly, whereas the latter offers a much smaller computational burden at the price of a more elaborate implementation.
The analytic uncertainty propagation method is only valid for Gaussian PDFs with small variances. The necessary results are stated below; see Ebbehøj et al. (2023) for a brief introduction and Yang and Lam (2019) for derivations and further details. The uncertainty in the ith set of natural frequencies and damping ratios can be quantified by the following variance:
where ${\mathbf{\Sigma}}_{{f}_{i},{\mathit{\zeta}}_{i}}\in {\mathbb{R}}^{\mathrm{2}\times \mathrm{2}}$, vec(⋅) is the vectorization operator transforming an M×N matrix $\mathbf{A}=[{\mathit{a}}_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{\mathit{a}}_{N}]$ to a column vector $\mathrm{vec}\left(\mathbf{A}\right)={\left[{\mathit{a}}_{\mathrm{1}}^{T},{\mathit{a}}_{\mathrm{2}}^{T},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{\mathit{a}}_{N}^{T}\right]}^{T}\in {\mathbb{R}}^{MN\times \mathrm{1}}$, and the vectorized posterior variance matrix $\mathrm{vec}\left({\mathbf{P}}_{{\mathbf{L}}_{t{}^{\prime}}}\right)$ corresponding to ${\mathbf{L}}_{t{}^{\prime}}$ is
where ${\mathit{\sigma}}_{{a}_{i}}^{\mathrm{2}}\left({\mathit{\xi}}_{t{}^{\prime}}\right)=\sum _{j=\mathrm{1}}^{{p}_{\mathrm{a}}}{\mathit{\sigma}}_{{a}_{i,j}}^{\mathrm{2}}{g}_{j}({\mathit{\xi}}_{t{}^{\prime}}{)}^{\mathrm{2}}$, and 0 is a zero matrix with dimensions as indicated by the subscript. Note that ${\mathit{\sigma}}_{{a}_{i}}^{\mathrm{2}}\left({\mathit{\xi}}_{t{}^{\prime}}\right)$ is not necessarily Gaussian, as the individual terms of the sum can be dependent. The outputs of the three functions are
where $\mathit{\rho}=[{\mathit{\rho}}_{\mathrm{1}},{\mathit{\rho}}_{\mathrm{2}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{\mathit{\rho}}_{{n}_{\mathrm{a}}}{]}^{T}$, $\mathit{\eta}=[{\mathit{\eta}}_{\mathrm{1}},{\mathit{\eta}}_{\mathrm{2}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{\mathit{\eta}}_{{n}_{\mathrm{a}}}{]}^{T}$, and $\mathit{\beta}=\mathit{\left\{}\right[{f}_{i,t{}^{\prime}},{\mathit{\zeta}}_{i,t{}^{\prime}}{]}^{T}:i=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{n}_{\mathrm{a}}\mathit{\}}$. The partial derivative of F_{1,i} with respect to $\mathrm{vec}\left({\mathbf{L}}_{t{}^{\prime}}\right)$ is
where a^{⊤} denotes the complex conjugate transpose of a, ⊗ denotes Kronecker's product, and w_{i} is the ith left eigenvectors of ${\mathbf{L}}_{t{}^{\prime}}$ satisfying
where w_{i} is a column vector. The partial derivatives of functions F_{2,i} and F_{3,i} are
where hats denote estimated values.
The procedure for computing natural frequencies and damping ratios and their uncertainties from AR coefficients is summarized in Table 2.
2.6 Considerations for practical implementation
In this section some considerations that may be useful for the practical application of the GPTARMA approach are listed. These are based on using the GPTARMA model for various applications, including the wind turbine blade response of an operational wind turbine during an extreme coherent wind gust (BHawC simulation; see Ebbehøj et al., 2023), field measurements from Siemens Gamesa's fleet, and the applications in the present paper:

Reduce bandwidth of response (output). Lowpass filter and subsequently down sample the response prior to GPTARMA modelling such that the reduced bandwidth only contains the modes of interest; i.e. unimportant response components at higher frequencies are filtered out. This reduces the computational cost of estimating the model and may help reduce the risk of overfitting in cases with low amounts of training data.

Filter out highfrequency content of EOV signals by lowpass filtering. This may allow for a better fit, as modal properties usually do not vary at high frequencies (many times per second) and may prevent highfrequency noise propagating from AR coefficients to modal parameters.

Use zerophase filtering for lowpass filtering of response and EOV signals. Consequently, the response and EOV signals are not phase shifted by filtering.

Add a small amount of artificial zeromean NID noise to the (lowpassfiltered) response signal prior to GPTARMA modelling. Adequate RMS noise levels relative to the signal are typically between 1 % and 5 %. This is an example of jittering and data augmentation in general. It can have a regularizing effect, i.e. reduce the risk of overfitting, and is commonly used in the fields of e.g. image analysis and for machine learning methods (Iwana and Uchida, 2021; Bishop, 1995; Ding et al., 2015; Erenler and Serinagaoglu Dogrusoz, 2019). According to Bishop (1995), adding noise to training data is equivalent to Tikhonov regularization.
This section presents an experimental setup of a shear frame (SF) structure with controlled timevarying damping and the experimental procedures used to generate response measurements for testing and validating the GPTARMA approach. In Ebbehøj et al. (2023), the method is verified using responses from two simulated cases: a simple twomode model response and an edgewise blade response during normal operation obtained with BHawC (Siemens Gamesa's inhouse aeroelastic code). This section provides an experimental validation based on measurements of a simple structure conducted in a controlled setting.
3.1 Experimental setup and procedures
An experimental test rig with a threestorey SF structure was used to test modal parameter estimation methods for structures with timevarying damping. The structure has a pair of voltagecontrolled electromagnetic dampers (EMDs) at each floor. Two types of tests were conducted: a baseexcitation test, where the SF structure was excited with a shaker table, and a hammer test, where the structure was excited with a (measured) hammer impact. Figure 1 illustrates the experimental setup and instrumentation, and Fig. 2 shows the SF structure. The floors are quadratic ($\mathrm{110}\times \mathrm{110}\times \mathrm{10}$ mm) and are constructed from aluminium with 2 mm thick copper plate inlays at the top and bottom to increase the effect of eddy current damping through EMDs (A in Figs. 1 and 2). The beams separating the floors are $\mathrm{150}\times \mathrm{18}\times \mathrm{0.5}$ mm and are realized by four steel rulers (C), which are fixed to each floor.
As can be seen in Fig. 1b, the shaker table (used for baseexcitation tests) is driven by an amplified, digital signal using a National Instruments (NI) relay module (NI 9481) and LabView as an interface with the computer. The EMDs are powered by a 24 V signal, which is controlled in an on–off manner using the NI relay module, which is controlled using FlexLogger software. The displacement response at each floor is measured by laser displacement sensors and sampled by an NI data acquisition (DAQ) module (NI 9215) via FlexLogger, along with relevant control signals, e.g. control signals for the EMDs and amplifier.
The input signal for the shaker table used for the baseexcitation test is a stochastic normally distributed broadband signal (i.e. lowpassfiltered white noise) to achieve stochastic excitation of the first three modes of the SF structure. The test is run for 60 min, all signals are collected using a sampling frequency of 50 Hz, and all EMDs are turned on and off during the test with periods of constant voltage ranging from 240 to 540 s.
For the baseexcitation test, the displacement response at the third floor is used and referred to as y(t) in Sect. 3, and the measured voltage over the EMDs v(t) is used as the EOV to predict changes in damping caused by the EMDs. Figure 3 shows the displacement response y(t) and the corresponding spectrogram and spectrum, along with voltage over the EMDs v(t), indicating when the EMDs are on (24 V) and off (0 V). From the spectrum it can be seen that the first three modes have frequencies of about 1.8, 5.3, and 7.8 Hz. The measurements sampled at 50 Hz are lowpass filtered (stop band frequency of 8.05 Hz) to prevent aliasing and are subsequently downsampled to 16.67 Hz to facilitate faster estimation of the GPTARMA model, without losing relevant information about the first three modes. In addition, artificial zeromean NID noise with an RMS level of 1 % is added to the downsampled signal.
To validate modal estimates obtained from the outputonly baseexcitation tests, the modal parameters are also estimated by standard impact hammer tests (e.g. Ewins, 2000; Halvorsen and Brown, 1977). The instrumentation for the impact hammer tests consists of an accelerometer (E) placed on the first floor, an impact hammer with a force transducer, and a B&K 3160 front end for data acquisition. The test is conducted by impacting each floor multiple times and measuring the resulting acceleration response. The accelerometer and force transducer signals are sampled by the B&K front end and analysed using B&K Pulse and ME'scope software to estimate modal properties via local frequency response function (FRF) curve fitting. The accelerometer is only mounted on the first floor during impact hammer tests and is thus dismounted prior to baseexcitation tests. During hammer tests, the shaker table is fixed. Hammer tests are conducted with all EMDs on and all EMDs off, and each test is conducted three times.
3.2 Model structure identification
A suitable model structure for modal analysis based on the thirdfloor response y(t) is identified using the procedure summarized in Table 1. In this controlled experimental setting, the excitation can be considered stationary since the shaker table is driven by a stationary NID signal. The functional subspace ℱ_{MA} representing the effect of the excitation therefore only consists of a constantvalue bias vector, and the innovation variance is assumed constant (corresponding to $K=N/\mathrm{2}$ in Eq. 14). To capture the steplike damping variations caused by the EMDs turning on and off, a firstorder Legendre polynomial in the voltage supplied to magnetic dampers (i.e. the EOV) is included in the functional subspace ℱ_{AR}, along with a constantvalued bias vector to account for the naturally occurring damping. If the voltage v(t) had changed gradually (i.e. not binary) from 0 to 24 V, ℱ_{AR} would likely need to include higherorder Legendre polynomials to model the voltage–damping relation. The functional subspaces for the AR and MA parts are populated as
where ${\mathit{A}}_{\mathrm{0}}\in {\mathbb{R}}^{N\times \mathrm{1}}$ and ${\mathit{C}}_{\mathrm{0}}\in {\mathbb{R}}^{N\times \mathrm{1}}$ are constantvalued bias vectors, $\mathit{v}\in {\mathbb{R}}^{N\times \mathrm{1}}$ contains measurements of the voltage over the EMDs for times $t=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}N$, and g_{2}(v) is a firstorder Legendre polynomial in v. The voltage signal is filtered with a centred moving average filter (better suited for step signals than regular lowpass filters) with a window width of 1 s to prevent measurement noise from propagating to the modal parameters (see Sect. 2.1).
The AR and MA model orders n_{a} and n_{c} are selected such that the predictive capabilities of the model in terms of the RSS/SSS and BIC (see Sect. 2.4) are converged, and the model captures modes of interest. In Fig. 4 RSS/SSS and BIC are seen to converge at about n_{a}=7 and n_{c}=2, constituting lower boundaries for the AR and MA model orders. The next step is to assess the model orders required to capture the modes of interest.
Inspecting the frequency stabilization diagram in Fig. 5, the three stable poles corresponding to the three natural frequencies of the SF can be seen to be stabilized at model orders n_{a}=11 and n_{c}=8, which are selected as the model orders used in downstream analysis. For the MA model order, the poles arguably already converge at n_{c}=4, but the frequency of the third mode changes slightly until fully stabilizing at n_{c}=8. This minor consideration can be taken because the amount of training data is large relative to the model complexity; i.e. the risk of overfitting is small. In this case the driving model selection criterion is the ability to capture the modes of interest, as is typically the case (see Sect. 2.4).
3.3 Model validation
In this section the model structure identified in Sect. 3.2 is validated using the procedure presented in Sect. 2.3, based on the model residual e_{t} $(t=\mathrm{1}+{n}_{m},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}N)$. The standardized residual z_{t} and the sample ACF can be seen in Fig. 6. Judging by the time series plot, the standardized residuals for both the training and the test set appear to be stationary stochastic. This observation is supported by the ACF, for which significant correlations only exceed the 95 % confidence limits for white noise at 3 % and 4 % of the time lags (i.e. less than 5 %) for the training and test data set. Furthermore, the sign test shows that the number of sign changes in the residual sequence is well within the 95 % confidence interval of sign changes for a sequence of Gaussian white noise of the same length.
The above whiteness tests, based on both standardized residuals and sign changes, suggest that the present GPTARMA model is valid and well suited for downstream analysis.
3.4 Results: modal parameter estimates
This section presents modal parameter estimates computed from the validated GPTARMA model of the thirdfloor displacement response y(t) and modal parameter estimates based on hammer tests. The GPTARMA model is estimated using 2700 s of data, corresponding to 44 820 samples (at a sampling frequency of 16.67 Hz), and about 4860, 14 310, and 21 060 oscillations of the first to third mode. Figure 7 shows the predicted response against the measured response, GPTARMA estimates of natural frequencies and damping ratios with 95 % confidence intervals, and SSI estimates for comparison. The SSI algorithm used is correlation driven (Peeters and De Roeck, 1999; Brincker and Ventura, 2015), and the stable poles corresponding to the three modes are selected manually from frequency stabilization diagrams. The SSI estimates are also based on the displacement response of the third floor y and are computed in windows of constant voltage over the electromagnetic dampers (eight windows seen in the plot). In addition, SSI estimates based on measurements in 30 s windows are shown for the sixth and seventh windows (2043 to 2523 s).
The predicted response in Fig. 7 is observed to model the measured response very well. The voltage over the electromagnetic dampers v(t) is overlaid in the plot, and 0 volt means no added damping from the electromagnetic dampers. A slight difference in response levels can be seen by comparing sections of the response with and without the electromagnetic dampers activated.
Good agreement for natural frequency estimates between GPTARMA and SSI is observed. The modal parameter estimates are also summarized in Table 3. Furthermore, the widest confidence intervals for the GPTARMA natural frequency estimates are ± 0.11, 0.03, and 0.02 Hz for the first, second, and third modes and seem representative of the spread of the SSI estimates, as seen in Table 3.
The figure and table show good agreement between SSI and GPTARMA damping estimates, where the best agreement is observed for the second mode. The recommended minimum measurement time (Brincker and Ventura, 2015) is dictated by the third mode (with EMDs off) and is about ${T}_{\mathrm{min}}=\frac{\mathrm{10}}{{f}_{\mathrm{3}}{\mathit{\zeta}}_{\mathrm{3}}}=\mathrm{214}$ s (based on the hammer test estimates in Table 3); i.e. the eight constantvoltage windows of minimum 240 s should be sufficiently long for adequate SSI estimates. The GPTARMA damping ratio estimates in Fig. 7 illustrate the idea of letting the AR coefficients, and consequently the modal parameters, be represented by EOVdependent basis functions. In this case, the voltage over the EMDs dictates the instantaneous changes in the damping estimates. In comparison, the 30 s SSI estimates (2043 to 2523 s) are seen to be inconsistent, and many deviate considerably from the other estimates. This illustrates the need for dedicated methods for shortterm damping estimation, such as the GPTARMA approach. Using the GPTARMA approach comes with the cost of specifying basis functions capable of representing (i.e. correlating with) the underlying causes of nonstationarity in the response, for which prior knowledge of the system is helpful. That is not needed for standard OMA methods such as SSI, although the validity of the LTI assumptions must be examined.
The widest damping estimate confidence intervals for the first to third mode are ± 6.7, 0.5, and 0.4 %; i.e. the damping estimate for the first mode is associated with the highest uncertainty in the three modes, as for the natural frequencies. This might reflect the fact that the training data contain more oscillation periods of the higher modes. To help elucidate this hypothesis, Fig. 8 shows the convergence of the mean damping ratio estimates and corresponding confidence intervals for all three modes in terms of oscillation periods. The results in Fig. 7 correspond to the estimates based on the largest training data set in Fig. 8. The figure indicates that only the third mode might have converged in terms of the confidence intervals, and the first mode in particular would seem to benefit from being estimated from more data. However, the mean estimates do not change much with the increase in training data. This suggests that the mean estimates might be representative despite the wide confidence intervals.
Table 3 shows that the GPTARMA and SSI natural frequency estimates agree well with the estimates from the three hammer tests. As for the damping ratios, the three hammer test estimates compare well to the GPTARMA and SSI estimates in the sense that they are well within the same order of magnitude. The hammer test estimates cannot be expected to agree entirely with GPTARMA or SSI estimates because the estimates are based on data from two fundamentally different experimental tests in terms of e.g. excitation (impulse vs. stochastic), the amount of energy input, and the temporal and spatial distribution of the energy input.
The SF test results presented in this section experimentally validate the efficacy of the GPTARMA method in providing representative shortterm damping estimates and illustrate its efficacy in the shortterm case compared to a traditional OMA method.
In this section the GPTARMA approach is applied to edgewise blade response measurements. Specifically, the method is used to estimate shortterm, linear equivalent modal damping of edgewise rotor modes, which are deliberately driven to flutterlike instabilities corresponding to negative damping values.
The experiments were conducted with an SWT7.0154 prototype wind turbine located at the DTU Wind Test Centre in Østerild, Denmark. The measurements collected in December 2018 were published by Volk et al. (2020), following up an instability field validation study (Kallesøe and Kragh, 2016). The edgewise blade modes were driven to flutterlike instabilities by momentarily allowing the rotor to run 10 % above the expected stability limit by changing the wind turbine controller parameters. As in Volk et al. (2020), the response used for the analysis is an edgewise blade bending response, obtained using a fibre Bragg optical strain sensor positioned 72 m outboard on the blade (i.e. close to the blade tip), measuring the bending response in the rotating frame of reference, which is referred to as y(t) in Sect. 4. All frequencies are normalized by the first edgewise blade frequency, and rotor speeds are normalized by the critical rotor speed of the first edgewise mode as in Volk et al. (2020).
Figure 9 shows the time series, spectrum, and spectrogram of the measured edgewise blade response y(t) along with normalized rotor speed. The figure shows a clear relation between high rotor speeds and the occurrence of “unstable” (high response levels) modes, which appear at normalized frequencies of 1.0 and 2.5. Strong response of these modes is especially prevalent at 150–230 s and 340–420 s.
The amount of available training data from this test is quite small relative to the model complexity (i.e. the number of parameters to be estimated) needed for the GPTARMA model to represent the complex response and underlying system. The bandwidth of y(t) has therefore been reduced compared to the measurement data, and artificial zeromean NID noise with a 5 % RMS of y(t) is added to the response y(t). These steps have been taken to use the limited measurements as efficiently as possible and to avoid overfitting (see Sect. 2.6).
4.1 Model structure identification
The model structure is identified following the procedure in Table 1. The functional subspace ℱ_{AR} used to represent the AR coefficients consists of a firstorder Legendre polynomial in the rotor speed Ω(t) along with a constantvalued bias vector. This allows the model to capture the effect of the varying rotor speed on the modal damping and the presumably timeinvariant contribution from the structural damping. The MA coefficients are represented by harmonic functions of the azimuth angle ψ(t) to account for the nonstationary 1P effect from the changing direction of gravity in the rotating frame. Thus, the functional subspaces for the AR and MA parts are defined as
where ${\mathit{A}}_{\mathrm{0}}\in {\mathbb{R}}^{N\times \mathrm{1}}$ and ${\mathit{C}}_{\mathrm{0}}\in {\mathbb{R}}^{N\times \mathrm{1}}$ are realvalued bias vectors, $\mathbf{\Omega}\in {\mathbb{R}}^{N\times \mathrm{1}}$ and $\mathit{\psi}\in {\mathbb{R}}^{N\times \mathrm{1}}$ contain rotor speed and azimuth angle measurements for $t=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}N$, g_{2}(Ω) is a firstorder Legendre polynomial in the rotor speed Ω, and h_{2}(ψ)=sin (ψ) and h_{3}(ψ)=cos (ψ) are the harmonic functions in the azimuth angle ψ. Thus, the EOVs used in the model are the rotor speed and azimuth angle, $\mathit{\xi}=[\mathbf{\Omega},\mathit{\psi}]\in {\mathbb{R}}^{N\times \mathrm{2}}$.
Both EOVs are zerophase lowpass filtered with a cutoff frequency of 0.3 Hz (see Sect. 2.6), and the innovation variance is estimated using Eq. (14) with a window length of 3.33 s; i.e. the GPTARMA model can estimate damping (and natural frequency) variations down to 3.33 s. The AR (MA) model order n_{a} (n_{c}) is selected based on the predictive capability (minimizing BIC and RSS/SSS) and ability to capture the modes of interest (see Sect. 2.4), which in this case have normalized frequencies of about 1 and 2.5. As in Sect. 3.2, the model structure is dictated by the ability to capture the modes of interest; i.e. the predictive capability converges at a lower model complexity. Figure 10 shows stabilization diagrams with respect to the AR and MA orders. The selected model orders are n_{a}=17 and n_{c}=5, as these are the lowest model orders at which the poles can be considered stable. RSS/SSS and BIC converge at n_{a}=5 and n_{c}=4 (plots not included).
4.2 Model validation
The model structure identified in Sect. 4.1 is validated in this section, based on analysing the model residual ${e}_{t}(t=\mathrm{1}+{n}_{m},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}N)$ as presented in Sect. 2.3. Figure 11 shows the time series, estimated ACF, and spectrum of the standardized residual z_{t} of the estimated GPTARMA model with the model structure identified in Sect. 4.1. The time series of the standardized residuals for the training data does not appear as stationary white noise since the amplitude (i.e. variance) exhibits timevarying behaviour, especially at 200 and 400 s. This is because the estimated innovation variance ${\widehat{\mathit{\sigma}}}_{e,t}^{\mathrm{2}}$ is larger at those instances, which are the times of the measured instabilities. Comparing the standardized residuals for the test data to those of the training data, the timevarying behaviour is much less distinct for the test data. This can be explained by all instabilities being contained by the training set. The ACFs also indicate that the standardized residuals do not resemble white noise, as the ACFs for both the training and the test data exceed the 95 % confidence level for white noise at 13 % of the time lags. The spectrum shows a distinct peak at a normalized frequency of 0.72, which coincides with 3 times the (average) rotor speed, indicating that the model does not adequately capture the 3P effect. However, the 3P frequency is well separated from the frequencies of the modes of interest, so it is not likely to affect modal parameter estimates. The spectra of the standardized residuals for the training and test data sets can be seen to differ, as for the time series. Because the instabilities are only present in the training data, the response characteristics of the two sets are different, which limits the validity of crossvalidation. In unison with the other whiteness measures, the number of sign changes in the residual sequences for the training and test data sets is not within the 95 % confidence intervals of the number of sign changes for the corresponding ideal white noise sequences. This supports the interpretation that the model residuals do not resemble white noise.
The residual analysis implies that the NID assumption of the innovations is violated; i.e. the model does not completely represent the statistical structure of the response. However, the frequencies at which the residuals have the largest magnitudes do not coincide with the frequencies of the two modes of interest. Thus, the modal parameters computed from the GPTARMA based on the present data set are not expected to be entirely accurate but may still offer some insight. More available measurement data (ideally including more instability measurements) for model training could improve the model accuracy by enabling more robust and accurate model parameter estimates. In addition, a singleoutput model (as the actual GPTARMA model) cannot account for the whirling effect since the frequencies of the forward and backwardwhirling modes coincide in a response measured in the rotating frame. Such modelform error might cause correlated model residuals.
4.3 Results: modal parameter estimates
Figure 12 shows the predicted response of the GPTARMA model and the corresponding modal parameter estimates with uncertainties. For direct comparison to Volk et al. (2020), the damping estimates are reported in terms of the logarithmic decrement (logdec) δ, which (for small damping) is related to the damping ratio ζ by δ=2πζ. The estimated normalized frequencies correspond well with peaks in magnitude seen in the spectrum and spectrogram of response y in Fig. 9, and the average 95 % confidence intervals for the first and second modes are ± 0.01 and 0.02. Both edgewise modes can be seen to become “unstable” (i.e. negatively damped) and most noticeable between 128–187 s and 341–405 s. The time instants of the mean damping estimates crossing zero damping coincide well with the amplitude envelope changing from exponentially increasing to decreasing or vice versa; i.e. the mean damping estimates are qualitatively meaningful. The average 95 % confidence intervals of the damping estimates for the first and second modes are ± 5.2 and 4.8 % logdec, corresponding to about ± 0.8 and 0.8 % in terms of critical damping.
The uncertainty associated with damping estimates may seem considerable, as the confidence intervals cover both positive and negative damping values during the flutterlike instabilities. The uncertainty might be reduced by addressing the two limiting factors alluded to in Sect. 4.2, namely the sparse training data and the inability of the singleoutput GPTARMA model to represent the whirling effect.
Figure 13 shows the damping estimates in Fig. 12 as functions of normalized rotor speed rather than time, which enables assessment of stability limits for the identified modes. The GPTARMA estimates are plotted against two sets of damping estimates from Volk et al. (2020), consisting of experimental estimates (based on the same response data) and predictions computed using HAWCStab2 (Hansen et al., 2018) based on a numerical model of the tested wind turbine. The experimental damping estimates are obtained using logarithmic decrement fits (i.e. exponential fits of the response envelope) of bandpassfiltered (near each resonance frequency) response signals. The exponential fits are performed on four sections of the data, two sections with an exponentially increasing (negative damping) amplitude and two with an exponentially decreasing (positive damping) amplitude (see Volk et al., 2020, for details).
The figure shows that the GPTARMA estimates agree quite well with the exponential fit estimates, especially near the critical rotor speeds, but the GPTARMA damping estimates tend to be higher than the exponential fit estimates. Table 4 shows the critical rotor speed estimated by each of the three approaches. In terms of the critical rotor speed, the HAWCStab2 predictions are slightly more conservative than the two experimental results. However, comparing the ${\mathit{\delta}}_{i}\left(\stackrel{\mathrm{\u203e}}{\mathrm{\Omega}}\right)$ curve slopes in Fig. 13, the HAWCStab2 results predict less rotor speed sensitivity, i.e. are less conservative with respect to how quickly the instabilities occur. But these discrepancies are small relative to the GPTARMA damping estimate uncertainties and the unquantified (but inevitable) uncertainties in the exponential fit estimates.
(Volk et al., 2020)(Volk et al., 2020)A recently proposed approach based on a Gaussian process timedependent autoregressive moving average (GPTARMA) model for shortterm damping (and natural frequency) estimation from outputonly vibration response measurements for vibrating structures influenced by environmental and operational variability has been experimentally tested and validated with two distinctly different experimental setups: a laboratory shear frame structure with timevarying damping properties achieved with electromagnetic dampers and a fullscale 7 MW wind turbine prototype which was deliberately driven to flutterlike instabilities. The primary idea of the GPTARMA approach is to condition the model parameters on measured time series of environmental and operational variables, which may enable shortterm tracking of system parameters like timevarying damping and natural frequencies.
An experimental setup consisting of a shear frame structure equipped with electromagnetic dampers was presented and shown to effectively realize a system with abruptly changing damping. Shortterm natural frequencies and damping ratios were estimated using the GPTARMA model and were shown to compare well to SSI and hammer test estimates in cases where the system was time invariant. Uncertainties were observed to be larger for the first mode compared to the second and third modes, but this could be explained by the first mode being trained on effectively less training data. GPTARMA damping estimates were compared to shortterm SSI estimates based on windows of 30 s measurements. The shortterm SSI estimates were observed to be inconsistent and deviating from the remaining estimates, which illustrated the effectiveness of the GPTARMA method for shortterm damping estimation relative to traditional OMA methods. The laboratory test validated the efficacy of the GPTARMA approach in shortterm damping (and natural frequency) estimation, given a sufficient amount of training data and a representative model structure.
The GPTARMA model was also tested using edgewise blade deflection measurements from a fullscale 7 MW wind turbine prototype during a flutter test. The first and second edgewise backwardwhirling modes were found to exhibit flutterlike instabilities, in agreement with a previous study. The mean damping estimates were considered qualitatively meaningful, as the GPTARMA model predicted negative damping for two modes coinciding in time and frequency with exponentially increasing vibration amplitude. The mean damping estimates also compared quite well with estimates from a previous study obtained from the same data. The estimated stability limits, i.e. the rotor speeds at which the damping becomes zero, showed quite good agreement with a previous study. However, the model validation implied that the model residuals did not resemble white noise, meaning that the GPTARMA model trained on the available data cannot be expected to be entirely accurate. The correlated model residuals and uncertainties of the damping estimates could potentially be reduced by training the GPTARMA model on more data and extending the GPTARMA model to a multipleoutput model to represent whirling modes better.
The GPTARMA approach appears to be an effective way of estimating shortterm damping based on outputonly measurements, given enough training data and a representative model structure. The use of GPTARMA models for analysing transient instabilities has been showcased. SSI and other standard OMA methods are easier to implement and apply than the GPTARMA approach since it does not require much prior knowledge of the system; i.e. these should be used for applications where the LTI assumptions are valid but may be inadequate for applications with considerable shortterm EOC variability.
Laboratory measurements are available upon request. Wind turbine test data are confidential. Code might be shared upon request.
Conceptualization and methodology: KLE, PC, and JJT; validation: KLE, PC, LMS, and JJT; investigation and visualization: KLE and LMS; writing (original draft): KLE and LMS; supervision and writing (review and editing): PC and JJT; software and formal analysis: KLE.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This work is partially funded by the Innovation Fund Denmark (grant 015300179B). The shaker table used in the laboratory setup was donated by COWIfonden (grant A155.01).
This research has been supported by the Innovationsfonden (grant no. 015300179B) and the COWIfonden (grant no. A155.01).
This paper was edited by Maurizio Collu and reviewed by two anonymous referees.
Au, S.K.: Operational Modal Analysis: Modeling, Bayesian Inference, Uncertainty Laws, Springer, Singapore, ISBN 9789811041181, https://doi.org/10.1007/9789811041181, 2017. a
AvendañoValencia, L. D. and Chatzi, E. N.: Multivariate GPVAR models for robust structural identification under operational variability, Probabil. Eng. Mech., 60, 103035, https://doi.org/10.1016/j.probengmech.2020.103035, 2020. a
AvendañoValencia, L. D. and Fassois, S. D.: Stationary and nonstationary random vibration modelling and analysis for an operating wind turbine, Mech. Syst. Signal Process., 47, 263–285, 2014. a
AvendañoValencia, L. D., Chatzi, E. N., Koo, K. Y., and Brownjohn, J. M.: Gaussian process timeseries models for structures under operational variability, Front. Built Environ., 3, 69, https://doi.org/10.3389/fbuil.2017.00069, 2017. a, b, c, d, e, f, g, h, i
AvendañoValencia, L. D., Chatzi, E. N., and Tcherniak, D.: Gaussian process models for mitigation of operational variability in the structural health monitoring of wind turbines, Mech. Syst. Signal Process., 142, 106686, https://doi.org/10.1016/j.ymssp.2020.106686, 2020. a
Bishop, C.: Pattern Recognition and Machine Learning, in: 1st Edn., Springer, New York, ISBN 9780387310732, 2006. a
Bishop, C. M.: Training with Noise is Equivalent to Tikhonov Regularization, Neural Comput., 7, 108–116, https://doi.org/10.1162/neco.1995.7.1.108, 1995. a, b
Bogoevska, S., Spiridonakos, M., Chatzi, E., DumovaJovanoska, E., and Höffer, R.: A datadriven diagnostic framework for wind turbine structures: A holistic approach, Sensors, 17, 720, https://doi.org/10.3390/S17040720, 2017. a, b
Brincker, R. and Ventura, C. E.: Introduction to Operational Modal Analysis, John Wiley and Sons, ISBN 9781118535141, https://doi.org/10.1002/9781118535141, 2015. a, b, c, d, e
Chen, C. and Duffour, P.: Modelling damping sources in monopilesupported offshore wind turbines, Wind Energy, 21, 1121–1140, 2018. a
Ding, Y., Xue, H., Ahmad, R., Chang, T. C., Ting, S. T., and Simonetti, O. P.: Paradoxical effect of the signaltonoise ratio of GRAPPA calibration lines: A quantitative study, Magn. Reson. Med., 74, 231–239, https://doi.org/10.1002/mrm.25385, 2015. a
Ebbehøj, K. L., Tatsis, K., Couturier, P., Thomsen, J. J., and Chatzi, E.: Shortterm damping estimation for timevarying vibrating structures in nonstationary operating conditions, Mech. Syst. Signal Process., 205, 110851, https://doi.org/10.1016/j.ymssp.2023.110851, 2023. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p
Erenler, T. and Serinagaoglu Dogrusoz, Y.: ML and MAP estimation of parameters for the Kalman filter and smoother applied to electrocardiographic imaging, Med. Biol. Eng. Comput., 57, 2093–2113, https://doi.org/10.1007/s11517019020186, 2019. a
Ewins, D. J.: Modal testing: Theory, Practice and Application, in: 2nd Edn., Research Studies Press, ISBN 0863802184, 2000. a
Fouskitakis, G. N. and Fassois, S. D.: Functional series TARMA modelling and simulation of earthquake ground motion, Earthq. Eng. Struct. Dynam., 31, 399–420, 2002. a
Halvorsen, W. G. and Brown, D. L.: Impulse technique for structural frequency response testing, Sound Vibrat., 11, 8–21, 1977. a
Hansen, M. H., Thomsen, K., Fuglsang, P., and Knudsen, T.: Two methods for estimating aeroelastic damping of operational wind turbine modes from experiments, Wind Energy, 9, 179–191, 2006. a
Hansen, M. H., Henriksen, L. C., Tibaldi, C., Bergami, L., Verelst, D., Pirrung, G. R., and Riva, R.: HAWCStab2 User Manual, http://hawcstab2.vindenergi.dtu.dk/download (last access: September 2023), 2018. a
Hansen, M. O., Sørensen, J. N., Voutsinas, S., Sørensen, N., and Madsen, H. A.: State of the art in wind turbine aerodynamics and aeroelasticity, Prog. Aerosp. Sci., 42, 285–330, 2006. a, b
Hu, W.H., Thöns, S., Rohrmann, R. G., Said, S., and Rücker, W.: Vibrationbased structural health monitoring of a wind turbine system Part II: Environmental/operational effects on dynamic properties, Eng. Struct., 89, 273–290, https://doi.org/10.1016/j.engstruct.2014.12.035, 2015. a
Iwana, B. K. and Uchida, S.: An empirical survey of data augmentation for time series classification with neural networks, PLOS ONE, 16, e0254841, https://doi.org/10.1371/journal.pone.0254841, 2021. a
Kallesøe, B. S. and Kragh, K. A.: Field validation of the stability limit of a multi MW turbine, J. Phys.: Conf. Ser., 753, 042005, https://doi.org/10.1088/17426596/753/4/042005, 2016. a
Kitagawa, G. and Gersch, W.: Smoothness priors analysis of time series, Springer, New York, ISBN 0387948198, ISBN 1461207614, ISBN 9780387948195, ISBN 9781461207610, 1996. a
Madsen, H.: Time series analysis, in: vol. 72, Chapman and Hall, ISBN 9780429195839, ISBN 142005967x, ISBN 9781420059670, ISBN 1322623546, ISBN 0429195834, ISBN 1420059688, ISBN 9781420059687, https://doi.org/10.1201/9781420059687, 2007. a, b
Murphy, K. P.: Probabilistic Machine Learning: Advanced Topics, MIT Press, ISBN 9780262048439, 2023. a, b
Peeters, B. and De Roeck, G.: Referencebased stochastic subspace identification for outputonly modal analysis, Mech. Syst. Signal Process., 13, 855–878, 1999. a
Peeters, B. and Roeck, G. D.: Stochastic system identification for operational modal analysis: A Review, J. Dynam. Syst. Meas. Control. Trans. ASME, 123, 659–667, 2001. a, b, c
Poulimenos, A. G. and Fassois, S. D.: Parametric timedomain methods for nonstationary random vibration modelling and analysis  A critical survey and comparison, Mech. Syst. Signal Process., 20, 763–816, 2006. a, b, c, d, e, f, g, h
Rasmussen, C. E. and Williams, C. K. I.: Gaussian Processes for Machine Learning, in: vol. 7, MIT Press, ISBN 026218253X, 2006. a
Spiridonakos, M. D. and Fassois, S. D.: Nonstationary random vibration modelling and analysis via functional series timedependent ARMA (FSTARMA) models – A critical survey, Mech. Syst. Signal Process., 47, 175–224, 2014. a, b, c
Spiridonakos, M. D., Poulimenos, A. G., and Fassois, S. D.: Outputonly identification and dynamic analysis of timevarying mechanical structures under random excitation: A comparative assessment of parametric methods, J. Sound Vibrat., 329, 768–785, 2010. a, b, c
Stewart, P.: Gram – Schmidt orthogonalization: 100 years and more, Numer. Lin. Algebra Appl., 20, 492–532, https://doi.org/10.1002/nla.1839, 2013. a
Veers, P., Bottasso, C. L., Manuel, L., Naughton, J., Pao, L., Paquette, J., Robertson, A., Robinson, M., Ananthan, S., Barlas, T., Bianchini, A., Bredmose, H., Horcas, S. G., Keller, J., Madsen, H. A., Manwell, J., Moriarty, P., Nolet, S., and Rinker, J.: Grand challenges in the design, manufacture, and operation of future wind turbine systems, Wind Energ. Sci., 8, 1071–1131, https://doi.org/10.5194/wes810712023, 2023. a, b
Volk, D. M., Kallesøe, B. S., Johnson, S., Pirrung, G. R., Riva, R., and Barnaud, F.: Large wind turbine edge instability field validation, J. Phys.: Conf. Ser., 1618, 052014, https://doi.org/10.1088/17426596/1618/5/052014, 2020. a, b, c, d, e, f, g, h, i, j, k
Wang, Z., Yang, D. H., Yi, T. H., Zhang, G. H., and Han, J. G.: Eliminating environmental and operational effects on structural modal frequency: A comprehensive review, Struct. Control Heal. Monit., 29, 1–24, 2022. a
Yang, J. H. and Lam, H. F.: An innovative Bayesian system identification method using autoregressive model, Mech. Syst. Signal Process., 133, 106289, https://doi.org/10.1016/j.ymssp.2019.106289, 2019. a, b, c
 Abstract
 Introduction
 Methods
 Laboratory test: shear frame structure with timevarying damping
 Fullscale wind turbine test: instability experiment
 Conclusions
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Methods
 Laboratory test: shear frame structure with timevarying damping
 Fullscale wind turbine test: instability experiment
 Conclusions
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References