Leading edge surface erosion is an emerging issue in wind turbine blade reliability, 1 causing reduction in power performance, aerodynamic loads imbalance, increased noise emission 2 and ultimately additional maintenance costs, and if left untreated, leads to the compromise 3 of the functionality of the blade. In this work, we ﬁrst propose an empirical spatio-temporal 4 stochastic model for simulating leading edge erosion, to be used in conjunction with aeroelastic 5 simulations, and subsequently propose a deep learning model trained on simulated data, which 6 aims to monitor leading edge erosion by detecting and classifying the degradation severity. The 7 main ingredients of the model include a damage process that progresses at random times, across 8 multiple discrete states characterized by a non-homogeneous compound Poisson process, which is 9 used to describe the random and time-dependent degradation of the blade surface, thus implicitly 10 affecting its aerodynamic properties. The model allows for one, or more, zones along the span of 11 the blades to be independently affected by erosion. The proposed model accounts for uncertainties 12 in the local airfoil aerodynamics via parameterization of the lift and drag coefﬁcients curves. 13 The proposed model is used to generate a stochastic ensemble of degrading airfoil aerodynamic 14 polars, for use in forward aero-servo-elastic simulations, where we compute the effect of leading 15 edge erosion degradation on the dynamic response of a wind turbine under varying turbulent 16 input inﬂow conditions. The dynamic response is chosen a deﬁning output as this relates to the 17 output variable that is most commonly monitored under a Structural Health Monitoring (SHM) 18 regime. In this context, we further propose an approach for spatio-temporal dependent diagnostics 19 of leading erosion, namely, a deep learning attention-based Transformer, which we modify for 20 classiﬁcation tasks on slow degradation processes with long sequence multivariate time-series as 21 inputs. We perform multiple sets of numerical experiments, aiming to evaluate the Transformer 22 for diagnostics and assess its limitations. The results reveal Transformers as a potent method for 23 diagnosis of such degradation processes. The attention-based mechanism allows the network to 24 focus on different features at different time intervals for better prediction accuracy, especially for 25 long time-series sequences representing a slow degradation process. 26 Keywords: Wind turbine; Structural Monitoring; Leading edge erosion; Poisson process; aeroelas- 27 tic simulations; Lift and drag; Deep learning; Transformer; Diagnostics 28 Version September 1, 2021 submitted to Journal Not Speciﬁed https://www.mdpi.com/journal/notspeciﬁedDRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 2 of 43 Nomenclature 29 a Wind shear exponent 30 E Expected value of a random variable 31 V Variance of a random variable 32 s Turbulence 33 CD Aerodynamic drag coefﬁcient 34 CL Aerodynamic lift coefﬁcient 35 Cm Aerodynamic moment coefﬁcient 36 Cp Aerodynamic pressure coefﬁcient 37 CNN Convolutional Neural Network 38 EZ Eroding Zone 39 LE Leading Edge 40 LEE Leading Edge Erosion 41 MLP Multilayer Perceptron 42 MTS Multivariate time-series 43 NHCPP Non-Homogeneous Compound Poisson Process 44 O&MOperation and Maintenance 45 R Blade radius 46 RNN Recurrent Neural Network 47 RV Random Variable 48 SHM Structural Health Monitoring 49 Ti Turbulence intensity 50 U Mean wind speed 51 WT Wind Turbine 52 1 Introduction 53 This work is motivated by current efforts on development of a novel micro-electro- 54 mechanical monitoring system based on aerodynamic surface pressure and aero-acoustic 55 measurements for Structural Health Monitoring of wind turbine blades, as part of 56 the Aerosense project1. Sensing nodes are distributed along the span of the blade 57 delivering measurements of local sectional aerodynamic pressure distributions, aero- 58 acoustic and acceleration signals (see illustration in Figure 1). Assuming availability of 59 such a monitoring system, we here propose an approach to the modelling and diagnostics 60 of erosion of the leading edge for wind turbine blades. 61 1https://www.aerosense.ai/DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 3 of 43 Figure 1. Illustration of Aerosense, a novel MEMS-based aerodynamic surface pressure and aero-acoustic measurement system for wind turbines blades. Bogdanoff and Kozin (1985) [ 1] deﬁne cumulative damage as the “irreversible ac- 62 cumulation of damage throughout life that ultimately leads to failure”. Leading edge 63 erosion of a wind turbine blade represents such a process, albeit at slow degradation 64 rates, with effects slow to manifest and even harder to detect early on. Most diagnostic 65 techniques today rely on direct visual inspection[ 2], image processing [ 3], and statistical 66 analysis, e.g. regression or clustering of SCADA output, such as blade pitch, rotor speed 67 or electrical power [ 4,5]. Visual inspection and image acquisition imply stopping the 68 turbines, resulting in loss of power production. More importantly, they give a local 69 (both spatial and temporal) snapshot of the erosion condition without a direct link to 70 the historical progression of the degradation and how it relates to environmental and 71 operational conditions. Even though it is better positioned to discriminate the envi- 72 ronmental and operational conditions versus degradation process, long term SCADA 73 based statistical analysis invariably do not possess the necessary resolution to determine 74 accurately the severity and spatial extent of erosion on the blades. A gap therefore exists 75 in the ﬁeld of spatio-temporal modelling and diagnostics of leading edge erosion for 76 wind turbines rotor blades. In this work, we are particularly interested in leveraging 77 time series emanating from the novel Aerosense based aerodynamic surface pressure 78 monitoring system for wind turbines blades. The aim is to (1) detect the occurrence 79 of erosion and (2) classify the severity (intensity) of erosion. This could be regarded 80 as either the inverse problem of reconstruction for the lift and drag coefﬁcients, or a 81 problem of ﬁnding the optimal discriminants for erosion versus other effects from time 82 series response signals. This is a time dependent multi-class and multi-label problem. 83 To tackle this problem and build meaningful diagnostics models, realistic time-series 84 data is required for training. Our options are either collect multi-year ﬁeld data, where 85 the rotor blade are known to undergo leading erosion, or develop a numerical model of 86 leading edge erosion coupled to an aeroelastic simulator. As data availability for such 87 systems is often limited and restricted by data sharing protocols, we here opt for the 88 route of simulation, with the additional beneﬁt of delivering a model that can be used 89 for simulating such slow leading edge erosion degradation processes. 90 Towards accomplishing the goal of diagnostics and inference of erosion of the 91 leading edge of wind turbine blades, this work makes the following contributions: 92 • We develop a stochastic spatio-temporal erosion model of the leading edge of wind 93 turbine blades, which is characterized by a non-homogeneous compound Poisson 94 process across discrete states, embedded in a generator of stochastic ensemble of de- 95 grading airfoil aerodynamic polars for use in forward aero-servo-elastic simulations. 96 The coupled model is able to compute the aeroelastic non-stationary response of a 97 wind turbine, thus reﬂecting its behaviour under the effect of leading edge erosion 98DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 4 of 43 and varying turbulent input inﬂow conditions over long period of degradation 99 time horizon. 100 • We adapt a deep learning multi-variate time-series based Transformer with atten- 101 tion mechanism to detect and classify long–term and slow leading edge erosion 102 processes assuming availability of on-blade sectional monitoring data, under short 103 and long-term wind inﬂow uncertainties and inject direct aerodynamic uncertain- 104 ties on the lift and drag coefﬁcients of the airfoil sections along the span of the 105 blade. 106 A graphical abstract of this research is presented in Figure 2. Figure 2. A graphical abstract of modelling and diagnosing erosion of the leading edge of wind turbine blades. 107 The remainder of this article is organized as follows. In Section 2 we offer a review 108 of the state of the art in modelling and diagnostics of leading edge erosion. In Section 3 109 we present the details of the spatio-temporal stochastic model for leading edge erosion. 110 In section 4 we present the uncertainty modelling and the aeroelastic simulations setup. 111 In section 5 we elaborate the theory of the deep learning multi-variate time-series based 112 Transformer for diagnostics and inference. In Section 6, we illustrate the novelty and 113 principles of the proposed framework on a simulations-based application and further 114 provide discussions and list the limitations in Section 7. 115 2 Review of Prior Art 116 We ﬁrst review modelling approaches of leading edge erosion in conjunction with 117 aeroelastic simulations and subsequently assess the state-of-the-art in diagnostics meth- 118 ods. 119 2.1 Modelling leading edge erosion 120 Several studies have quantiﬁed the impact of leading edge erosion on aerodynamic 121 and aeroelastic performance of wind turbines. 122DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 5 of 43 The general thread in the ﬁrst set of studies amounts to describing the constitu- 123 tive laws, physical processes and instigators at the micro scale resulting in surface 124 erosion/degradation and the various stages and types of erosion severity e.g., [ 6–10]. 125 Such models incorporate the effect of wind speed, air density, particle size, incubation 126 time, and erosion intensity due to rainfall, snowfall, sea spray, and fog. It is unclear 127 however how the micro-scale models couple to the macroscale aspects of leading edge 128 erosion and consequently the effect on turbine performance. 129 The second set of studies quantiﬁes the impact on wind turbine performance via 130 computational approaches (aeroelastic or CFD simulations) or extensive experimental 131 campaigns in wind tunnels, coupled to uncertainty propagation schemes e.g., [ 11–16]. 132 The limitations in all these works is that they do not allow for multiple zones along the 133 span of the blades to be independently affected by stochastic erosion processes. 134 Bortolotti et al. [17] proposes a simple stochastic model to describe the extent of 135 blade span-wise degradation due to leading edge erosion via a single factor that is 136 assumed to follow a truncated beta distribution. However, the long term temporal 137 degradation process of blade leading edge is not accounted for. Dimitrov [18] describes 138 a ﬁrst attempt to introduce temporal progression to the damage growth for leading edge 139 erosion due to rainfall events indirectly via loss of annual energy production. 140 To the best of our knowledge, we ﬁnd a gap in existing literature for comprehensive 141 models that encompass both stochastic spatial and temporal elements to describe the 142 erosion process of the leading edge of wind turbine blades and the effects on their 143 performance over long periods of degradation time horizons. 144 2.2 Diagnosing leading edge erosion 145 We here review state-of-the-art in leading edge erosion diagnostics and discuss 146 potential approaches from recent work in probabilistic and machine learning research to 147 solve the problem of spatio-temporal dependent diagnostics. The review is organized 148 along two axes: (1) whether or not a forward model is known and to what extent it 149 is used in training and testing and (2) whether or not the learning is supervised or 150 unsupervised, i.e., whether or not the training relies on access to matched ground truth 151 and measurement pairs. 152 Different approaches exist to diagnosing the extent of leading edge erosion on wind 153 turbine blades. Traditionally, manual human inspection is used, although modern visual 154 approaches are being developed to utilize drones paired with computer vision [ 19]. 155 Visual inspection, however, requires the turbine to be shutdown and is highly dependent 156 on weather conditions. Our aim is to use data-driven continuous monitoring methods to 157 perform diagnostics without resorting to turbine shutdown. We review different machine 158 learning and probabilistic approaches that have been put forth and which, given data 159 recorded by sensors placed on the wind turbine (e.g. novel aerodynamic pressure sensors, 160 SCADA, accelerometers, strain gauges, etc.), would be able to estimate the extent of 161 leading edge erosion on wind turbine blades. In this context, the diagnostics task can 162 be described as a multivariate time-series (MTS) classiﬁcation/clustering problem; an 163 algorithm is tasked with predicting the degradation status/class of a blade by analysing 164 multiple time-dependent signals. We therefore include MTS methods, which have 165 not been speciﬁcally applied to wind turbine blades in the past, allowing us to explore 166 promising machine learning approaches from diverse ﬁelds unrelated to structural health 167 monitoring: activity recognition, ﬁnance, anomaly detection, sound (speech, music) 168 analysis, seismology, meteorology, human health signal (EEG and ECG) interpretation, 169 etc. If the forward model of the degradation process is known, the diagnostics problem 170DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 6 of 43 becomes a supervised classiﬁcation task, seeing that each MTS data point is associated 171 with its degradation ground truth. Data-driven MTS classiﬁcation is a ﬁeld which is 172 currently generating a lot of attention, with many promising methods that have been 173 recently developed [ 20]. Within these approaches, we can differentiate between those 174 that use the raw time-series data [ 21], known as end-to-end approaches, and those which 175 ﬁrst extract features that are then passed to a standard classiﬁcation algorithm [22]. 176 Extracting features from time-series data can be accomplished in many ways, with 177 domain dependant feature engineering or statistical methods often being the ﬁrst choices. 178 Algorithms based on Empirical Mode Decomposition (EMD) are implemented for rep- 179 resenting MTS as a superposition of simpler well-behaved components called Intrinsic 180 Mode Functions (IMFs). EMD approaches are suitable for damage detection and degra- 181 dation state recognition in the analysis of nonlinear and non-stationary signals. However, 182 they have shown low performance in highly-stochastic signals with spike pulses and 183 jumps [ 23,24]. Distance-based metrics, which evaluate the similarity between time- 184 series data form another popular class of conventional feature construction methods 185 [25]. Alternatively, bag-of-pattern models [ 26] [27] can also be used to generate features 186 from sequential inputs. These methods work by dividing the time-series sample into 187 subsequences that are then turned into an ensemble of discrete identiﬁers and ﬁnally 188 a word count histogram, performed on the ensemble, is used to generate a feature 189 vector. In [ 28] and [ 29], feature vectors are built from MTS data by passing the inputs 190 through numerous random convolutional kernels. This approach, named ROCKET, has 191 a relatively low computational cost, yet, when combined with a standard ridge classiﬁer, 192 it achieves state-of-the-art accuracy on the UCR time-series archive dataset [ 30]. Other 193 methods turn the time-series inputs into feature images [ 31], which are then fed into 194 Convolutional Neural Networks (CNNs) [ 32] designed for image classiﬁcation, thus 195 taking advantage of the recent advances in computer vision performance. Although 196 feature engineering and feature extraction can potentially lead to accurate diagnostics 197 with good explainability, we are interested in implementing end-to-end methods that 198 use raw MTS data. This is due to the fact that feature-based methods often require large 199 amounts of cumbersome human intervention, which limits generalization. Robustness 200 is another potential downside, where changing the input, due to sensor malfunction 201 for instance, might require re-selecting relevant features, while it can be factored into 202 end-to-end approaches, provided it is part of the training data. 203 Conventional (’shallow’) data-driven classiﬁcation methods do not often use un- 204 transformed MTS as inputs, due to the difﬁculty associated with extracting relevant 205 information from noisy, complex data. In contrast, deep learning methods often excel 206 at learning hidden discriminative features, making them particularly well suited for 207 end-to-end approaches [ 33]. CNNs have been applied to end-to-end MTS classiﬁca- 208 tion in the context of prognostic health management [ 34], surgeon skill classiﬁcation 209 [35] and other tasks [ 36]. CNNs are limited by their receptive ﬁeld, making their use 210 limited to relatively short time-series. Recurrent neural networks (RNNs) [ 37] are, by 211 design, effective at solving tasks relating to sequential data. Long Short Term Memory 212 (LSTM) Networks [ 38], currently the most commonly used type of RNNs, have been 213 used to successfully classify clinical measurements [ 39] and action recognition tasks [ 40]. 214 However, RNNs are also limited to relatively short input sequences due to memory 215 requirements and are difﬁcult to train. In recent years, Transformers [ 41] have attracted a 216 lot of attention in the natural language processing community and have mostly replaced 217 RNNs. These neural networks leverage the power of attention mechanisms to learn 218 the relevant interdependencies within sequences. Transformers have increasingly been 219 applied in other domains, such as computer vision [ 42], time-series forecasting [ 43], and 220 multivariate time-series classiﬁcation [ 44]. Graph Neural Networks [ 45], a class of deep 221 data-driven methods which utilizes graph-structured data, have also been applied to 222 end-to-end MTS classiﬁcation [ 46]. This is done by ﬁrst constructing an adjacency matrix 223DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 7 of 43 from the untransformed MTS data. In [ 47], spatio-temporal Graph Transformer Neural 224 Networks (GTNN) are implemented in order to effectively capture both dynamic spatial 225 and temporal trends. We could envisage using a GTNN model to leverage temporal and 226 spatial correlations in the degrading sections of the three blades on a wind turbine. 227 As previously mentioned, deep end-to-end models are in general better suited 228 for generalisation than feature-extraction approaches. Moreover, they can be further 229 optimised for transferability [ 48], for instance by applying transfer learning methods, 230 such as using domain adaption layers [49]. 231 When the forward model of the degradation process is unknown, other methods 232 for ground truth gathering must be considered, such as hand-labeling individual sam- 233 ples; a process which is often prohibitively time-consuming for large datasets. Instead, 234 assuming that the labels are unknown, one can use unsupervised MTS clustering ma- 235 chine learning methods. In this type of approach, the data-driven algorithm learns to 236 automatically cluster MTS samples into groups that exhibit similar degradation states 237 [50]. Most MTS clustering methods are based on feature-extraction. However, it can be 238 noted that in an unsupervised learning context, feature-extraction has a considerable 239 impact on clustering performance [ 51], and it is therefore difﬁcult to rely on standard 240 domain-agnostic methods to construct features. In response to this difﬁculty, an emerg- 241 ing trend has been to learn the optimal feature representation of MTS data for clustering 242 using neural networks [52], [44]. 243 It is also worth mentioning unsupervised MTS anomaly detection methods [ 53], 244 as these binary clustering approaches could be potentially extended to diagnostics. 245 Anomaly detection for MTS is often based on similarity metrics [ 26,54,55], although 246 these approaches are limited when it comes to capturing non-linearity and contextual 247 anomalies (i.e. behavior that is abnormal only given the time at which it occurs). Neural 248 networks have proven to be effective at detecting non-linear anomalies [ 56]. LSTMs 249 were used to detect web trafﬁc anomalies in [ 57], and [ 58] proposed the use of CNNs 250 to detect abnormal behavior at multiple scales. General Adversarial Networks [ 59], a 251 learning framework where two neural networks compete in a zero-sum game, have also 252 been applied to anomaly detection [ 60]. Contrastive learning [ 61] is another approach 253 to unsupervised clustering from multivariate time series signals for fault detection 254 and diagnostics in the context of changing operating conditions and novel fault types. 255 Franceschi et al. [52] showed the scalability of contrastive learning via triplet losses on 256 long time series without labeling and without structural assumptions on non-temporal 257 data. It is expected however that the performance of contrastive learning to suffer 258 under highly stochastic input and operating conditions as is the case in wind turbines 259 applications. 260 In the context of data-driven wind turbine blade condition monitoring, a number of 261 methods have been reported in the literature [ 62]. In [ 63] vibrational frequency response 262 functions were used as input features to neural networks, in order to detect changes in the 263 structural health of a blade. LSTMs were trained to classify faults, based on multivariate 264 sensor recordings of a model wind turbine test rig in [ 64], while [ 65] implemented a 265 CNN to detect possible wind turbine blade breakages based on SCADA data. In [ 66], 266 Support Vector Machines were trained on features built from active acoustic recordings, 267 to diagnose healthy and damaged states of a model wind turbine. Another feature- 268 extraction approach was combined with decision trees in [ 67] to diagnose different fault 269 scenarios for model wind turbine. Finally, ice buildup on wind turbine blades was 270 predicted with random forests trained on SCADA data in [ 68], while Weijtjens et al. 271 [69] used Gaussian Processes trained on data emanating from low-cost tower-mounted 272 sensors to do so. 273DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 8 of 43 Considering their extensive use in other domains, Transformers bear great potential 274 for SHM damage detection. Some researchers have made use of attention mechanisms 275 similar to those found in Transformers for bearing remaining useful life prediction 276 [70], but to the best of the authors’ knowledge, only in [ 71] have proper Transformer 277 models been used in this context. In this work, however, an end-to-end approach 278 is not used; instead, the model is fed a pre-processed input composed of time series 279 features such as Mel Frequency Cepstral Coefﬁcients (MFCCs) and Short Time Fourier 280 Transformation. We choose Transformers as the base of our end-to-end approach for 281 diagnostics. However, we implement additional modiﬁcations so as to be able to handle 282 multivariate sequences on the order of several 10’s of thousands of data points while 283 still being able to track long term degradation progression, reﬂected by newly sampled 284 sequences with possible discontinuities. Few, if any, of the existing MTS learning models 285 discussed above are capable of this. 286 3 Modelling Erosion of the Leading Edge 287 3.1 Non-Homogeneous Compound Poisson Process 288 Leading edge erosion (LEE) is caused by environmental variations in the blade 289 surface temperature oscillations, moisture, UV radiation, raindrops, sand, hailstones 290 or other particles impacting the leading edge of the blade. It could also be initiated by 291 surface cracks due to global strain from blade ﬂexing [ 8], errors during manufacturing 292 in terms of paint deposits, consistency of the gel coat thickness and bonding strength, 293 or initiating surface damage during blade handling in transport or installation. This 294 causes the surface material to be removed from the blade surface, leaving a rough proﬁle 295 that degrades the aerodynamic performance and, on the long term and if left untreated, 296 impacts the structural integrity of the blade. Blade erosion may be more predominant 297 in the sectional regions, where the local ﬂow velocity is higher as this induces more 298 energetic impacts with soiling or eroding elements. This is further exacerbated by the 299 relatively greater impact the outboard region has on overall rotor performance. In 300 addition, the Reynolds number and inﬂow turbulence inﬂuence the airfoil transition 301 location and subsequently the aerodynamic performance of the wind turbine blades. As 302 a result, it is not plausible to simulate all possible stages of LE erosion and uncertain 303 inﬂow conditions; we therefore opt for a more pragmatic approach. Erosion initially 304 initiates on the blades in the form of pits near the leading edge on the pressure side. 305 These pits develop gradually over time into gouges, then steadily grow in their size and 306 density to coalesce as delaminations. Gaudern [7]identiﬁes ﬁve categories of severity 307 based on proprietary inspection reports and photographs of Vestas turbines that have 308 been operating for up to 5 years. In Sareen et al. [6]we ﬁnd an analysis of the effect of 309 erosion at the leading edge on the aerodynamic performance of the DU96W180 310 airfoil, and categorized erosion from pits to delamination in different types and stages, 311 as shown in Table 1. Each combination of type and stage of erosion has associated 312 with it aerodynamic polar curves ( CL,CD,CM), derived via measurements in a wind 313 tunnel [ 6,7,11]. Throughout the paper we adopt the labeling of the leading edge erosion 314 severity shown in Table 2. 315 We propose a two-step empirical spatio-temporal stochastic model of LEE. Step 1 316 of the model includes a damage process that occurs at random times emulating the non- 317 stationary time of arrival of degradation and its magnitude, with 10 discrete states (what 318 we call erosion severity as shown in Table 2) based on Non-Homogeneous Compound 319 Poisson Process (NHCPP) to describe the random degradation of the blade surface 320 and implicitly its aerodynamic properties [ 6]. In step 2 of the model, the local leading 321 edge degradation manifests by physical changes in the aerodynamic polars. This step 322DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 9 of 43 Table 1: Types and stages of leading edge erosion with the number of pits (P), number of gouges (G), and magnitude of leading edge delamination (DL) [6]. Type A Type B Type C Stage 1 100P — — Stage 2 200P 200P/100G — Stage 3 400P 400P/200G 400P/200G/DL Stage 4 — 800P/400G 800P/400G/DL+ Stage 5 — —- 1600P/800G/DL++ Table 2: Class labels adopted in this research for the severity categorization of leading edge erosion. Severity class Type & stage 0 Clean airfoil 1 Type A, Stage 1 2 Type A, Stage 2 3 Type A, Stage 3 4 Type B, Stage 2 5 Type B, Stage 3 6 Type B, Stage 4 7 Type C, Stage 3 8 Type C, Stage 4 9 Type C, Stage 5 generates an ensemble of stochastic aerodynamic polars for a given erosion severity of 323 LEE in NHCPP by parameterizing the lift and drag coefﬁcients curves [72]. 324 The stochastic degradation process of a wind turbine blade LEE is given by a 325 non-homogeneous Poisson process (NHCPP) over a ﬁnite time horizon. The NHCPP 326 fN(t),t0ghas rate of arrivals l(t)on the time interval t2[0,T]. A sample path 327 of a NHCPP emulates two main properties in the proposed LEE: (1) random time of 328 arrival, 0t1t2...: these are the time instances corresponding to change in 329 the degradation level of an airfoil aerodynamic properties (so-called shocks), and (2) 330 magnitude of change fYn,n1gin the degradation level of an airfoil aerodynamic 331 properties at the nthtime of arrival. In practice, for a probabilistic assessment, many 332 NHCPP paths are sampled. The model incorporates one or more zones (adjacent or 333 disjoint) along the span of the blade(s) affected by erosion, but are assumed to undergo 334 independent degradation processes. The NHCPP modify the deﬁnition of a Poisson 335 process so that it can incorporate time-dependent rate. The non-homogeneous Poisson 336 process is suitable given the fact that the rate of LEE is not constant in time (time- 337 dependent arrival rates); on the short time scales (hourly), it is dependent on the rotor 338 speed due to variations in the wind speed, while on the longer time scales it is dependent 339 on seasonal variations such as rain intensity, and changes in environmental temperatures. 340 We thus infer that the rate of occurrence of shocks is periodic, as shown in Figure 5. 341 In addition, LEE consists of limited and ﬁnite number of stages resulting in varied 342 and non-stationary incremental effects (severity) on the airfoils aerodynamic polars. It 343 is in fact pragmatic to assume that the severity of LEE has a ﬁnite number of stages, 344 because these could be categorised (e.g. via ﬁeld observations), and their impact on 345 aerodynamic properties emulated and quantiﬁed via measurements in wind tunnels 346 as demonstrated in [ 6]. The impact of LEE on the lift and drag coefﬁcients of the 347 NACA 64618 airfoil [11,14,73] for 10 severity classes is shown in Figures 3 and 4. 348 The non-homogeneous Poisson process is similar to an ordinary Poisson process, 349 except that the average rate of arrivals per period l(t)is allowed to vary with time 350 (Figure 5). A classic deﬁnition of the Non-Homogeneous Poisson Process is formulated 351DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 10 of 43 -30 -20 -10 0 10 20 30 Angle of Attack [deg.]-1.5-1-0.500.511.5Lift Coefficient [-]Clean LEE Severity 1 LEE Severity 2 LEE Severity 3 LEE Severity 4 LEE Severity 5 LEE Severity 6 LEE Severity 7 LEE Severity 8 LEE Severity 9 Figure 3. Lift coefﬁcients for NACA 64618 for 10 severity classes of leading edge erosion. -15 -10 -5 0 5 10 Angle of Attack [deg.]0.0050.010.0150.020.025Drag Coefficient [-] Figure 4. Drag coefﬁcients for NACA 64618 for 10 severity classes of leading edge erosion. as follows [ 74]. Let l(t):[0,¥)![0,¥)be an integrable function; the counting process 352 fN(t),t0gis called a non-homogeneous Poisson process with rate l(t)if all the 353 following conditions hold: 354 1. N(0)=0 355 2.8t,s0, and 0ut,Nt+sNtis independent of Nu 356 3.8t,s0,P(Nt+sNt=0)=1l(t)s+o(s) 357 4.8t,s0,P(Nt+sNt=1)=l(t)s+o(s), and 358 5.8t,s0,P(Nt+sNt2)=o(s) 359 where, sis a a very short interval of time. Here o(s)is a function that is negligible compared to s, ass!0. Assuming o(s) =g(s), then: lim x!0g(s) s=0 One approach for generating Non-Homogeneous Poisson processes is the “process 360 analogue” of acceptance-rejection called thinning, which is the scheme we adopt in our 361 model. The procedure is as follows: 362DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 11 of 43 1. Choose lusuch that l(t)lu8t2[0,T] 363 2. Initialize t=0 and I=0 364 3. Generate u1U(0, 1) 365 4. Set t t1 lulogu1366 5. If t>T, stop; Else go to next 367 6. Generate u2U(0, 1), independent of u1368 7. If u2l(t) lu, set I=I+1;S(I)=t 369 8. Go to Step 3 370 In the above procedure, the ratiol(t) luis the so-called thinning probability . In our 371 model luis chosen as the maximum of the rate function l(t)(e.g. see Figure 5). For long 372 degradation periods, it might be necessary to break [0,T]into small intervals and pick a 373 lufor each interval in order to avoid a high rejection rate. This is known as piece-wise 374 thinning. At the end of this procedure, the event times S(I)(arrival times) and the 375 counting process are yielded according to the Non-Homogeneous Process. Finally, for a 376 non-homogeneous Poisson process with rate l(t), the number of arrivals in any interval 377 is a Poisson random variable; however, its parameter can depend on the location of the 378 interval. More speciﬁcally, we can write: 379 N(t+s)N(t)PoissonZt+s tl(a)da (1) 0 6 12 18 24 30 36 42 48 54 60 66 72 78 84 90 96 102 108 114 120 126 132 138 144 150 156 162 168 174 180 186 192 198 204 210 216 222 228 234 240 Time [months]00.050.10.150.20.250.3Rate of Occurence, (t) Figure 5. Cyclic rate function l(t)for the non-homogeneous Poisson process, representing the effect of weekly and seasonal variations on the rate. l(t)represents the rate of change in damage of LE erosion per given period. Introducing the compounding effect to the Non-Homogeneous Poisson Process 380 yields the so-called Non-Homogeneous Compound Poisson Process (NHCPP). The Com- 381 pound Poisson Process, which is a very useful extension of a Poisson process, replaces 382 the unit jumps of a Poisson process with random jump sizes [ 75]. The jumps magnitude 383 at the ntharrival time has value Yn,n1 attached to it. The successive Ynare assumed 384 to be independent and identically distributed positive random variables each having 385 the same distribution, and are also assumed independent of the underlying Poisson 386 counting process of shocks. We assume that Y=(Y1,Y2, . . .)is a sequence of indepen- 387 dent, identically distributed, real-valued random variables, and that Yis independent 388 of the underlying Poisson process. The compound Poisson process associated with 389 the given Poisson process fN(t),t0gand the sequence Yis the stochastic process 390 Z=fZt,t0g, where: 391 Zt=N(t) å n=1Yn (2) The damage is always positive, i.e., that P(Yn0)=1, and the damage accumulates 392 additively according to Ztin Equation 2. In our model Ztrepresents the cumulative 393DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 12 of 43 amount of erosion incurred on the blade section at time N(t). In our model, we cate- 394 gorised 10 erosion severity levels reﬂected in the lift and drag coefﬁcients, as shown 395 in Figures 3-4 (other categorization could be speciﬁed as needed). The ﬁnal level of 396 degradation leads to the loss of function of the blade, not in the sense of a catastrophic 397 loss (e.g. rupture), but rather implying that the blade loses its ability to generate lift 398 efﬁciently. At this stage, we assume that once the blade section has reached level 9 of 399 erosion, the process stops under the assumption that any further degradation becomes 400 limited, since all the surface protective coatings and paint have already delaminated. 401 We adopt a classical approach whereby Y1,Y2, ...,Ynare each exponentially distributed 402 according to density: 403 gY(y)=mexpmy(3) where mis the mean jump (shock) magnitude. It follows that the sum Y1+Y2+...+Yn 404 is Gamma distributed G(n)(z). For details, we refer the interested reader to [ 76,77]. 405 Formally, we adopt a truncated exponential distribution for Ynin Equation 3. The logic 406 being that large erosion magnitude jumps (e.g. large jumps in erosion severity) in one 407 shock event are physically hard to justify. As a result, the upper truncation limit for 408 the shock is set to a magnitude of 4 in our model. Furthermore, as more events arrive 409 between [0,T], and damage is compounded Y1+Y2,+..., we adjust the upper truncation 410 limit for the shock magnitude from 4 to 1 when fI3^ ZmaxåN(t)=I n=1Yn 3g. 411 Finally, sampled jumps of magnitude Yn>0fortn<Tare not allowed once the 412 compounded damage has reached its highest compounded severity class ( Zmax=9) 413 before the end of time horizon Thas been reached. This is reasonable as any additional 414 shocks (e.g. rain or hail) cannot not incur any additional erosion, as there are not any 415 more surface materials that can be easily eroded. 416 Finally, we introduce wind speed and blade location dependencies to the model. For 417 a given time of arrival of a shock, our model associates higher damage (1) in the event 418 when the shock is concurrent with higher inﬂow wind speeds due to higher momentum 419 (and thus kinetic energy) of the impacting particles anywhere on the blade and (2) in 420 the outboard sections of the blade due to the local higher relative speed of the ﬂow. We 421 require the mean jump (shock) magnitude mto comprise a linearly increasing function 422 of blade radial location and quadratic with wind speed. The ﬁnal outcome is thus an 423 NHCPP degradation paths as shown in Figure 6, where 5 zones along the span of blade 424 2 are chosen to undergo leading edge degradation. Even though the exact pattern and 425 location of surface changes during operation is a random process, we see that our model 426 tends to provide higher damage rates of the leading edge on the outboard sections of the 427 blade compared to the inboard sections. 428 3.2 Stochastic aerodynamic Polars model 429 In this section, we generate a stochastic ensemble of degrading airfoil aerodynamic 430 polars for use in forward aero-servo-elastic simulations. The outputs of the NHCPP 431 are the degradation paths describing different degradation zones along the span of the 432 rotor blades. Each arriving shock at a given time along these paths corresponds to a 433 jump of a certain magnitude in the leading edge erosion severity. A jump in the erosion 434 severity corresponds implicitly to a degradation of the aerodynamic properties of the 435 airfoils affected by erosion along the span of the blade(s). However, uncertainties in 436 the degrading airfoil aerodynamics lift and drag coefﬁcients curves still remain to be 437 accounted for. The stochastic aerodynamic Polars model, that we here describe, accounts 438 for such uncertainties by parameterizing the lift and drag coefﬁcients curves. Consider- 439 able effort and capital is invested in predicting the static aerodynamic lift CL, drag CD440 and moment CMcoefﬁcients of airfoils as accurately as possible under various operating 441 conditions. The lift and drag coefﬁcients are then used by the wind turbine designers 442DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 13 of 43 0 16 32 48 64 80 96 112 128 144 160 176 192 208 224 240 Time [months]-1012345678910LEE Severity ClassBlade #2, 5 Eroding Zone(s) Eroding Zone = 1 Eroding Zone = 2 Eroding Zone = 3 Eroding Zone = 4 Eroding Zone = 5 Figure 6. A sample of the NHCPP of the leading edge erosion on 5 zones along the span of blade 2 over a 20-year period (240 months). Section 5 is the outermost section along the span of the blade. as input to aero-servo-elastic simulations to predict power performance, extreme and 443 fatigue loads in addition to stability margins in normal and extreme operating and stand 444 still conditions. An airfoil’s static aerodynamic data are almost exclusively derived 445 from measurements acquired in wind tunnel tests. However, an airfoil section on the 446 wind turbine rotor operates in 3-dimensional, un-steady and turbulent inﬂow under the 447 guidance of a control system, none of which are accounted for in static wind tunnel tests. 448 Some aspects of uncertainty in airfoil data (surface roughness, 3D corrections, effect of 449 Re numbers, wind tunnel measurements or geometric distortions) have been studied in 450 [72]. The general consensus is that uncertainties in airfoil data do affect a wind turbine’s 451 performance and structural loading. Consequently, our proposed LEE model also takes 452 into account the inherent uncertainty in airfoil static lift and drag coefﬁcients during 453 the erosion period via a stochastic model of static airfoil lift and drag polar curves. The 454 details of this model can be found in Abdallah et al. [72], with a short summary presented 455 herein for conciseness. 456 A well known result from Blade Element Moment theory links the aerodynamic lift 457 (CL) and drag ( CD) coefﬁcient to axial ( a) and tangential ( a0) induction factors, as follows: 458 a 1a=s0(CLcosf+CDsinf) 4 sin2f a0 1a=s0(CLsinfCDcosf) 4lrsin2f(4) where s0is the rotor solidity, Fis the angle of the incoming relative wind with the rotor 459 plane, l=wr V0,wis the rotor speed, ris the radial distance from the rotor center, and V0is 460 the freestream wind. Our approach to affecting aerodynamic induction is by introducing 461 a stochastic model of the lift and drag coefﬁcients curves in BEM. 462 Several sources of uncertainties affect the lift and drag coefﬁcients. These uncertain- 463 ties associated with assessment of airfoil characteristics in wind tunnels, uncertainties 464 due to 3D ﬂow correction, uncertainties stemming from surface roughness, uncertainties 465 related to the blade geometric distortions in manufacturing and handling, uncertainties 466 related to the blade geometric distortions when deﬂected under load, uncertainties 467 due to the effects of Reynolds number, uncertainties associated with extending airfoil 468DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 14 of 43 aerodynamic characteristics to post stall, and ﬁnally uncertainties stemming from the 469 validation of airfoil data by ﬁeld full scale prototype test. It is not possible to quantify 470 the joint distribution of all these random variables and as a result a more simpliﬁed 471 approach is chosen via a stochastic model as proposed in [ 72]. The stochastic model 472 consists in parameterizing the lift coefﬁcient curve by the slope in the linear range¶CL ¶a, 473 the point indicating the start of the trailing edge separation (AoA TES,CL,TES), the point 474 of maximum lift (AoA max,CL,max)and the point where the stall recovery is initiated 475 (AoA SR,CL,SR). The drag coefﬁcient is parameterized by a bias at low angles of attack, 476 and by the point where max drag coefﬁcient occurs at AoA =90where it exhibits the 477 largest variations. Figure 7 shows samples of the stochastic lift and drag coefﬁcients for 478 NACA 64618 for an erosion of severity class 5. These perturbations result in modiﬁed CL479 and CDcurves that maintain the primary characteristics of the original curve, but differ 480 in both magnitude and feature location, emulating possible uncertainties affecting the 481 aerodynamic properties of the airfoil section under real operation. 482 -200 -100 0 100 200-1.5-1-0.500.511.5Lift Coefficient [-] -15 -10 -5 0 5 10 15 Angle of Attack [deg.]00.010.020.030.040.05Drag Coefficient [-]Generated Original Figure 7. Samples of the generated stochastic lift and drag coefﬁcients of airfoil NACA 64618 corresponding to erosion severity class 5. Putting it all together, the NHCPP coupled to the generator of stochastic ensemble of 483 aerodynamic polars yields realistic degrading aerodyanmic polars over the degradation 484 time horizon, for every degrading zone and airfoil section on the rotor blades. For 485 example, Figure 8 shows the degrading CLover a 20-year period (240 months) for 486 forNACA 64618 airfoil. Figures 9-11 represent the corresponding degrading max lift 487 coefﬁcient, degrading angle of attack at max lift coefﬁcient and degrading lift coefﬁcient 488 ata=5, respectively. 489DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 15 of 43 Figure 8. NHCPP generated degrading static CLforNACA 64618. This is used as input to the aeroelastic simulations for Blade 1, section 5. For each pre-selected degrading section on a blade a set of such curves are generated. 0 16 32 48 64 80 96 112 128 144 160 176 192 208 224 240 Time [months]0.60.811.21.41.61.82CL, max [-] Figure 9. Degradation of the max lift coefﬁcient for Blade 1, section 5. 3.3 Overview of the Algorithm 490 The main ingredients of the model include a (1) spatio-temporal stochastic damage 491 process and (2) a generator of stochastic ensemble of degrading airfoil aerodynamic 492 polars. The algorithm is summarized in Algorithm 1. 493 4 Uncertainty Modelling and Aeroelastic Simulations 494 Two categories of uncertain random variables (RV) are considered: wind inﬂow 495 and aerodynamic effects. We have elaborated on the aerodynamic uncertainties in the 496DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 16 of 43 0 16 32 48 64 80 96 112 128 144 160 176 192 208 224 240 Time [months]681012141618 at CL, max [deg.] Figure 10. Degradation of angle of attack at max lift coefﬁcient for Blade 1, section 5. 0 16 32 48 64 80 96 112 128 144 160 176 192 208 224 240 Time [months]00.20.40.60.811.21.4CL at = 5 deg. [-] Figure 11. Degradation of the lift coefﬁcient at a=5for Blade 1, section 5. previous section. In this section, we detail the wind inﬂow uncertainties and the aeroe- 497 lastic simulations setup. The conventional wind turbine structural damage computation 498 process utilizes aeroelastic load simulations for wind turbines under normal and extreme 499 operation and wind inﬂow turbulence conditions, whereby the extreme and fatigue load 500 cycles determined over a short period of time at each mean wind speed are extrapo- 501 lated over the full expected lifetime. Such a process does not consider the changes in 502 inﬂow and environmental changes over long time periods, indeed relevant to LE erosion 503 process. 504DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 17 of 43 Algorithm 1: Spatio-temporal stochastic model of leading edge erosion. Input :fCL,CDg: static lift and drag coefﬁcients Result:fCL0MC,CD0MCg: degrading lift and drag coefﬁcients 1Set the blades Biaffected by LE erosion: Bi,i2f1, 2, 3g; 2Set the eroding zones EZi,jforBi:Zi,j,j2f1, 2, 3, 4, 5g; 3Set the erosion time horizon [0,T]; 4Associate airfoils Ai,j,kto erosion intensity range N= [0, 9],k2f1,Kg; 5foriter=1, . . . , ndo 6 Generate turbulent inﬂow samples for [0,T]:U,Ti,a,Y, andS; 7 foreach blade B ido 8 foreach Zone EZ i,jdo 9 foreach Airfoil A i,j,kdo 10 //NHCPP 11 while t=0, ..., Tdo 12 do: sample arrival times N(t); 13 do: sample erosion jumps Yn; 14 do: compound the erosion jumps, Zt; //Inject smooth temporal transition 15 do:[ytrans] =GeneralizedLogisiticFunction () 16 end 17 //Generate stochastic airfoil data samples 18 foreach Erosion severity Z tin a given NHCPP path do 19 do: Associate severity Zt!f CL,CDgofAi,j,k; 20 do:[CL0MC,CD0MC] =Generate _Stochastic _Aero _Data(); 21 end 22 end 23 end 24 end 25end 4.1 Inﬂow RVs and their Stochastic Models 505 The variation in the structural dynamic response of wind turbines is signiﬁcantly 506 dependent on the turbulent inﬂow wind ﬁeld conditions, including the mean wind speed, 507 turbulence, wind shear, and inﬂow skewness. In accounting for these inﬂuences, we 508 introduce the following RVs in the simulations setup: Mean wind speed, U, turbulence 509 intensity, Ti, wind shear, a, horizontal inﬂow skewness, Y, and the vertical inﬂow 510 skewness, S. 511 The mean wind speed follows a Weibull distribution UWBL (AU,KU), truncated 512 to[425]m/s, with parameters speciﬁed as follows: 513 E(U)=8.5,where A U=2E(U)pp KU=2.0(5) The conditional dependence between the turbulence sUand the mean wind speed Uis 514 deﬁned in the Normal Turbulence Model described in the wind turbine design standard 515 [78]. Here, we elect to use a reference ambient turbulence intensity Ire f=.16(the 516 expected value of the turbulence intensity at 15 m/s is called Ire f). This dependency is 517 given by the local statistical moments of suLN msU,s2 su as: 518 E sujU =Ire f(0.75u+3.8) V sujU = 1.4Ire f2(6)DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 18 of 43 As a result, the turbulence intensity is expressed as: Ti=su u. 519 The wind proﬁle above ground level is expressed using the power law relationship, 520 which deﬁnes the mean wind speed uat height Zabove ground as a function of the 521 mean wind speed uh, measured at hub height Zhas reference: 522 u uh=Z Zha (7) where ais a constant called the shear exponent. The conditional dependence between 523 the wind shear exponent aN ma,s2 a and the mean wind speed Uis given by [79]: 524 E ajU =0.088(ln(u)1) V ajU =1 u2 (8) We deﬁne a custom conditional dependence between the inﬂow horizontal skew- 525 nessYand the mean wind speed U, truncated to [11, 11 ]deg.,YN mY,s2 Y : 526 E YjU =ln(u)3 V YjU =15 u2 (9) The conditional dependence between the inﬂow vertical skewness Sand the mean 527 wind speed Uand turbulence intensity Ti, truncated to [6, 6]deg., is designed such as: 528 SN mS,s2 S :( mS=2,s2 S=1.0, if Ti10% & a0.1. mS=1.5,s2 S=1.0, otherwise.(10) The blade leading edge degradation scenario occurs over a 240-months period; 529 consequently, we introduce an additional element of uncertainty by allowing a time- 530 varying aspect of the inﬂow conditions, namely, the mean wind speed E(U), the shape 531 parameter KUof the Weibull distribution, and the reference ambient turbulence intensity 532 Ire fas shown in the example in Figure 12. 533 We choose to sample the wind inﬂow and the aerodynamic RVs using the Sobol 534 Quasi-Random sequences, which are designed to generate samples that are uniformly 535 distributed over the unit hypercube, i.e., as uniformly as possible over the multi- 536 dimensional input space [ 80]. In total, we sample 1200 joint samples of the wind 537 inﬂow RVs over the degradation period, as shown in Figures 13 - 14. A sample of U, 538 sU,a,YandSis used to generate a realization of inﬂow turbulent wind ﬁeld time 539 series, combined with a sample of stochastic and degrading CLand CDform the input to 540 theOpenFAST aero-servo-elastic environment to simulate the corresponding aeroelastic 541 response of the wind turbines structure under continuous erosion of the leading edge of 542 the blades. 543 4.2 Aero-servo-elastic Simulations Setup 544 Our coupled aero-servo-elastic simulations with the proposed NHCPP leading 545 edge degradation model setup is based on the OpenFAST simulator. OpenFAST [ 81], 546 is a coupled aero-hydro-servo-elastic analysis tool for modeling wind turbines. The 547DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 19 of 43 0 1 2 3 4 5 6 7 8 9 10 11 Weather period66.577.588.599.51010.511Mean wind speed, E[U] [m/s] 1.71.751.81.851.91.9522.052.12.152.2 Weibull shape parameter, KUE[U] KU (a) 0 1 2 3 4 5 6 7 8 9 10 11 Weather period0.080.10.120.140.160.18Reference turbulence, Iref (b) Figure 12. Long term inﬂow variations: (a) Varying mean wind speed and shape parameter KUof the Weibull distribution. (b) Varying reference turbulence Ire f. primary use of OpenFAST is to run nonlinear time-domain simulations. The simulations 548 in OpenFAST considered the reference NREL three-bladed up-wind, horizontal-axis 549 wind turbines with 126mrotor diameter, 5MW rated power and hub height of 90m. 550 The layout of the airfoils along the span of the blade are shown in Figure 15. We also list 551 some of the more important properties of the simulated wind turbine in Table 3. 552 OpenFAST is a wind-turbine-speciﬁc nonlinear time-domain simulator that em- 553 ploys a combined modal and multibody dynamics formulation, adopting limited degrees 554 of freedom (DOF). Large structural elements (blades and tower) are characterized by 555 properties such as stiffness and mass per unit length, which account for ﬂexibility char- 556 acteristics. OpenFAST models the turbine dynamics using 24 degrees of freedom (DOFs), 557 including two blade-ﬂap modes and one blade-edge mode per blade, two fore-aft and 558 two side-to-side tower bending modes, nacelle yaw, the generator azimuth angle and the 559 compliance in the drive train between the generator and hub/rotor. The aerodynamic 560 model is based on the unsteady Blade Element Momentum theory, including skew 561DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 20 of 43 Figure 13. Joint samples of the environmental sources of uncertainties for degradation period 0. E(U)=8.5m/s,KU=2.0,Ire f=0.09. Figure 14. Joint samples of the environmental sources of uncertainties for degradation period 5. E(U)=9.0m/s,KU=1.95, Ire f=0.12. inﬂow, dynamic stall and generalized dynamic wake. The Blade aerodynamic polars are 562 provided as a-priori input and are used as lookup tables or for interpolation. In our case, 563 the aerodynamic polars correspond to the degrading lift and drag coefﬁcients generated 564 in the NHCPP , such as the one shown in Figure 8. The stochastic input wind ﬁeld 565 uses the Kaimal turbulence model [ 82]. A common practice is to generate a signiﬁcant 566 number of stochastic simulations for various operating and environmental conditions 567 in order to cover variability on aeroelastic fatigue and extreme load analysis due to the 568 stochastic nature of the input wind ﬁeld and variability in weather patterns over long 569 time horizons ( >1-year). 570DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 21 of 43 0 5 10 15 20 25 30 35 40 45 50 55 60 Blade span [m]01234567Chord [m]DU40 DU35 DU30 DU25 DU21 NACA64618 Figure 15. Distribution of airfoils along the span of the blade. Example of the NREL5MW reference wind turbine. Table 3: Properties of the NREL 5-MW reference wind turbine. Number of blades 3 Rotor diameter 126 m Hub height 90 m Rated power 5 MW Cut-in wind speed 3 m/s Cut-out wind speed 25 m/s Control Variable Speed, Collective Pitch Variable speed from cut-in to cut-out wind speed Variable pitch from cut-in to cut-out wind speed Rated wind speed 11.4 m/s Cut-in and rated RPM 6.9 12.1 RPM In our simulation setup, we assign a maximum of 5 independent zones that could 571 undergo leading edge erosion per blade. The location (normalized by the blade length) 572 of the degradation zones along the span of the three blades are: 573 EZ:[0.65, 0.68 ],[0.70, 0.80 ],[0.82, 0.85 ],[0.87, 0.90 ],[0.92, 1.00 ] We simulate 20 years (240 months) worth of wind turbine operation with blades LEE. 574 We generate one sample every (approximately) six days, resulting in a total of 1200 575 simulated samples. Every sample is a 600-second multi-variate time series with 228 576 sensor signals each, sampled at 100Hz. Few examples of the time series signals are 577 shown in Figures 16- 18. 578 100 200 300 400 500 600 700 Time [s]05101520Wind Speed [m/s] Figure 16. Sample hub height wind speed. (Black) beginning and (Red) end of degradation period.DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 22 of 43 100 200 300 400 500 600 700 Time [s]00.20.40.60.811.2CL [-] Figure 17. Sample lift coefﬁcient at the airfoil proﬁle corresponding to 95% of the span of blade 1. (Black) beginning and (Red) end of degradation period. 100 200 300 400 500 600 700 Time [s]20004000600080001000012000Blade Root Flap BM [kNm] Figure 18. Sample ﬂapwise bending moment at the root of blade 1. (Black) beginning and (Red) end of degradation period. 4.3 Retained Output from Aero-servo-elastic Simulations 579 As mentioned in the beginning of the paper, the motivation for this research is the 580 Aerosense research project, where we are developing a novel MEMS-based aerodynamic 581 surface pressure and aero-acoustic smart measurement system for wind turbines blades, 582 for use in Structural Health Monitoring tasks. In this simulated setup, we assume one 583 Aerosense node located at 0.96 of the blade length, located within eroding zone 5 (near 584 the tip of the blade). From our aero-servo-elastic simulations, we elect to retain those 585 output signals which are expected to emanate from the Aerosense measurements nodes, 586 namely, lift and drag coefﬁcients, angle of attack, and wind speed. 587 Table 4: Retained sensors output of aero-servo-elastic simulations. Sensor name Description Time Time steps of the simulations Wind 1VelX X-direction wind velocity at hub-height B1N9Cl Lift force coefﬁcient at Blade 1, Aerosense Node at 0.96 R B1N9Cd Drag force coefﬁcient at Blade 1, Aerosense Node at 0.96 R B1N9Alpha Angle of attack at Blade 1, Aerosense Node at 0.96 R 5 Diagnosing LEE via Transformers 588 The simulation pipeline described above allows us to generate data which replicate 589 the sensory output of Aerosense nodes. Thus, we are able to produce datasets to train 590 a deep learning method, with the end goal being to detect and estimate the extent of 591 leading edge erosion on wind turbine blades. We choose to implement a modiﬁed 592 version of a popular class of sequential deep learning models: Transformers [41]. 593DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 23 of 43 To motivate the diagnostics task, we consider the case where CLtime series signals 594 emanate from two sections along the span of the blade. Figure 19 shows a scatter plot 595 of time series of CLat sections 0.96Rvs0.75Rof the blade for wind speeds varying 596 between 616m/sand no LEE (clean blade). Figure 20 shows a similar plot for a 597 ﬁxed wind speed 11m/sand evolving severity of leading edge erosion at section 0.96R. 598 An important property is that the relational dependency of the CLbetween the two 599 sections will change over time either due to changes in the severity of the leading edge 600 erosion or due to short-term variation in inﬂow and operating conditions. The problem 601 is thus to detect leading edge erosion and its severity, which could be regarded as a 602 problem of ﬁnding the optimal discriminants for erosion versus other effects (inﬂow, 603 operational conditions, aerodynamic uncertainties, etc.) from multivariate time-series 604 response signals. 605 Figure 19. CLfor wind speeds varying between 6 16m/sand no blade erosion. Figure 20. CLfor ﬁxed wind speed U=11m/sand evolving severity of leading edge erosion. The following are aspects we wish to build into the diagnostics method for LEE 606 classiﬁcation, taking into account the nature/circumstance of an operating wind turbine 607 in the ﬁeld: 608DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 24 of 43 • Labeled data for LEE are hard to acquire and scarce; as a result, any method must 609 be designed not to suffer from overﬁtting under scarce labeled data availability. 610 • Diagnostics shall be done with remote streaming of sensor data. Human interven- 611 tion and turbine down-time should be alleviated. 612 • The data used in the diagnostics method shall replicate the sensory output of 613 a single micro-electro-mechanical-systems based aerodynamic and aero-acoustic 614 measurement node located near the tip of one wind turbine blade. 615 • The method shall be capable of ingesting 10 minute long multivariate time-series 616 (industry standard SCADA recording length), sampled at 100 Hz, resulting in 617 60’000 data points long sequences. A sampling rate of 100 Hz may seem excessive, 618 but the aim is to capture even small turbulence scales. 619 • Exploits a supervised pre-training scheme to initialize a Transformer-based network 620 by utilizing large-scale labeled data from aeroelastic simulations of the turbine and 621 a degradation model as presented earlier in the paper. 622 • Once pre-training is completed, the pre-trained model can be further adapted to 623 the speciﬁc classiﬁcation task by ﬁne-tuning the model parameters on small-scale 624 labeled data. In this way, the general knowledge and representations about LEE 625 can be transferred to a label-scarce task, thereby improving the generalization 626 performance of the model as well as reducing the risk of overﬁtting. 627 • Physics-constraints shall be built-in the loss or likelihood function. 628 • The output predictions should be probabilistic in nature. 629 5.1 Experiments and Datasets 630 We utilize multivariate time-series data generated via the coupled simulation 631 pipeline previously described. The term multivariate is used in the sense that we 632 are monitoring multiple signals pertaining to lift, drag, angle of attack and inﬂow ve- 633 locity emanating from an Aerosense node located at one section at the tip of the blade. 634 Multivariate could also be thought of in the sense that we can take advantage of the pres- 635 ence of multiple Aerosense stations at different locations along the span of the blade(s), 636 however, this is beyond the scope of the current work. After some initial testing, we 637 have opted to utilize the four variables shown in Table 4 as inputs to the Transformer. 638 For instance, no performance beneﬁt was found when including simulated acceleration 639 data. 640 Three sets of experiments, each with different data requirements, are performed 641 aiming to answer the following questions: 642 1. Are we able to learn the LEE severity classes from MTS data, in the general machine 643 learning sense, with balanced data classes and no prior knowledge? 644 2. In a continuous monitoring context, are we able to diagnose jumps in LEE severity and 645 therefore identify the degradation path that the system takes? 646 3. Are we able to do so in a realistic setting, with all previously described sources of uncertainty 647 present in the simulations? 648 Thus, the goal of the ﬁrst set of experimental tasks is to train the neural network 649 on a generalized dataset with reduced uncertainty and where each degradation class 650 appears in equal measure. This is not representative of a real degradation path, where 651 the inherent stochasticity of the NHCPP may result in some classes appearing very 652 brieﬂy, but we aim to avoid biasing the network towards a random predominant class. 653 Moreover, this set of experiments also aims to understand how data availability affects 654 classiﬁcation performance, and how separable the classes are. Two datasets, correspond- 655 ing to experiments 1.1 and 1.2, are generated: one with the full 10 severity classes and one 656DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 25 of 43 with a reduced amount of severities (levels 0, 1, 6 and 9). For this set of experiments, no 657 information on previous states is used. In each of the sub-experiments, 4800 aero-elastic 658 simulations are gathered, then split into train, validation and test sets in a 70/20/10 659 split. 660 The second set of experiments aims to assess whether using full degradation paths 661 in a continuous monitoring setup to train the diagnostics method is a suitable strategy. 662 Here, the datasets are comprised of full NHCPP degradation paths, albeit with reduced 663 uncertainty. Due to the stochastic degradation, the datasets do not have balanced classes. 664 Assuming that a continuous monitoring system is in place, we therefore have access to 665 the previous degradation states. This allows for degradation monotonicity to be enforced 666 in the prediction, thus constraining the output to physically possible solutions (i.e., it 667 is physically impossible to have a state which is less degraded than previous states - 668 except for the case of direct service intervention for repair and maintenance, which is out 669 of scope in this work). Another objective is to evaluate whether severity grouping is a 670 viable strategy. Indeed, grouping the LEE stages by type (see Table 1) could be beneﬁcial 671 as a coarse predictor, if it proves to be sufﬁciently accurate. Thus, in experiment 2.1, all 672 10 degradation levels are used as labels, while in experiment 2.2 the levels are grouped 673 by type. Three full NHCPP degradation paths (new sample every 6 days, over a 20 year 674 period) with 1200 simulations each are separated into training and validation sets in a 675 80/20 split. A full degradation path, with 1200 data points, is reserved as the test dataset 676 for ﬁnal evaluation. 677 The goal of the third set is to evaluate how the model would perform in a realistic 678 continuous monitoring situation. Here, the dataset is also made of full NHCPP simu- 679 lations, but in this case, all possible stochasticity is turned on. Again, we assume that 680 we have access to the previous degradation states. We also evaluate class grouping, 681 separating the experiment into two. Three full NHCPP degradation paths (new sample 682 every 6 days, over a 20 year period), each consisting of 1200 MTS data points, are divided 683 into training and validation sets based on a 80/20 split. The test dataset is comprised of 684 one full degradation path, also with 1200 data points unused in training or validation. 685 Table 5 summarizes the differences between the datasets of the three sets of experi- 686 ments, including the partition of data between the different subsets. 687 Table 5: Dataset properties of the different diagnostic experiments. ParameterExperiments 1.1 1.2 2.1 2.2 3.1 3.2 NHCPP severities 0-9 0, 1, 6, 9 0-9 0-9 0-9 0-9 Severity type grouping - - - 0, A, B, C - 0, A, B, C Stochastic degradation - - X X X X Inﬂow turbulence X X X X X X Aerodynamic uncertainty - - - - X X Weather variability - - - - X X Num. training samples 3360 3360 2880 2880 2880 2880 Num. validation samples 960 960 720 720 720 720 Num. testing samples 480 480 1200 1200 1200 1200DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 26 of 43 5.2 Transformer Architecture 688 Here, we introduce the architecture of our Transformer neural network, used to 689 classify long sequence multivariate time-series data and infer the degradation status of 690 the leading edge. 691 Transformers are models that rely on self-attention mechanisms to highlight and 692 learn dependencies within sequences. Typically, self-attention is implemented by parsing 693 the input sequence into into key (K), query (Q), and value (V) vectors. The attention 694 weights on the values are then obtained by taking the scaled dot products of the query 695 with all keys, and then by applying a softmax function. Using the classical notation [ 41], 696 this writes: 697 Attention (Q,K,V) =softmax (QKT pdk)V (11) where dkis the dimension of the key and query vectors and is used to scale the dot 698 product. In practice, hscaled dot-products are used in parallel, within Multi-Head 699 attention layers. A Multi-Head model is able to attend jointly to information gathered 700 from different parts of the input, providing an improvement in representation over a 701 single attention mechanism. The ﬁnal attention vector is obtained by concatenating each 702 of the dot-product results: 703 MultiHead (Q,K,V) =Concat (head 1, ...,head h)WO(12) with 704 head i=Attention (QWQ i,KWK i,VWV i) (13) where the projections are parameter matrices WQ i2Rdmodeldk,WK i2Rdmodeldk,Vi2 705 Rdmodeldvand WO2Rhdvdmodel. 706 Multi-Head attention layers are combined with fully connected networks and layer 707 normalization to form the base architecture of each layer of the Transformer stack. We 708 use a novel time-windowing Transformer model which is based on the patch system 709 used in the Vision Transformer (ViT) model [ 42]. This architecture, as shown in Figure 21, 710 aims to address the quadratic attention bottleneck, an obstacle which limits the length 711 of sequences a Transformer can use. Indeed, the computational complexity of the self- 712 attention layers in Transformer models scales as O(L2), where L is length of the input 713 sequence. Given our requirements, the input time-series have a length of L=600000, 714 which is prohibitively resource intensive for standard Transformer models. To alleviate 715 this issue, we propose to divide the data into N windows along the time-dimension, 716 before passing each window through the learnable input embedding. This allows us to 717 use the full input time-series without resorting to downsampling which would entail 718 losing information. 719 Before dividing the MTS into windows, each of the individual channels/variables 720 is normalized in time. Then, the input embedding encodes each time-series window 721 into a vector of size dmodel , which is the latent size used throughout the Transformer. 722 The type of input embedding used has a considerable impact on performance and many 723 embedding types can be envisaged. We choose to use a learnable linear embedding for 724 simplicity, based on initial testing and the literature [ 42]. A class token is concatenated 725 to the embedded sequence, such that its state at the output of the Transformer is used 726 to infer the degradation severity. Class tokens are commonly used in natural language 727 processing tasks [ 83], and are an effective way to retrieve categorical information from 728 Transformers. Traditionally, positional encodings are then added to the sequence in order 729 to instil directional information into the model. Contrary to standard NLP Transformer 730 models, we do not use a positional encoding. Two positional encoding methods were 731DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 27 of 43 tested (classical sine functions and simple linear encoding), however, both of these 732 methods hindered classiﬁcation performance. 733 The classiﬁcation head consists of a single hidden layer Multilayer Perceptron 734 (MLP) that processes the class token at the output of the Transformer. It outputs a vector 735 with a size equal to the number of possible degradation classes, which, when passed 736 through a softmax activation, gives the likelihood scores of each class. For each class 737 output xi, the softmax likelihood writes: 738 softmax (xi) =exp(xi) åjexp(xj))(14) If we aim to predict the degradation at multiple zone along the blade, we can 739 simply stack multiple MLPs, one for each zone.Table 6 summarizes the values used in 740 the Transformer architecture. 741 Table 6: Parameters of the Transformer architecture. Number of windows, n 300 Window size, w 200 Internal Transformer dim., dmodel 256 Transformer stack size, s 6 Num. self-attention heads, h 8 Self-attention head dim., dhead 64 Output MLP dim., dMLP 2048 5.3 Loss functions 742 We ﬁrst train and evaluate the Transformer model on the datasets with equal 743 amounts of degradation classes, in the sense of a traditional classiﬁcation problem. Thus, 744 the objective function for this ﬁrst task is the standard cross-entropy loss: 745 LCE(y, ˆy) =n å i=1yilog(ˆyi) (15) where yis a one-hot label vector indicating the correct class, ˆyis a vector containing the 746 predicted softmax probabilities for each class and n is the number of classes. 747 In the second and third set of experiments, it is assumed that the state of the system 748 is known for the previous sampling period. We use this information to enforce physical- 749 ity; we add a second term to the loss function that penalizes a predicted degradation 750 class that is lower than the known previous class. This second term is the margin ranking 751 loss: 752 LMR(ˆc, cprev) =max(0,ˆccprev+m) (16) where ˆcis the predicted class, cprevis the degradation class of the previous known state, 753 and mis the margin hyperparameter. Overall, the objective function for the second and 754 third experiments is therefore: 755 Ltot=LCE+aLMR (17) where ais a hyperparameter used to balance the two components. 756DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 28 of 43 Figure 21. Time-windowing architecture of the novel Transformer model. The input multivariate time-series is split along the time dimension into windows which are then individually passed through the learnable linear embedding. The learnable class token that is added to the sequence is used in the classiﬁcation MLP head to predict the level of degradation. 5.4 Training regime 757 During training, we use a batch size of 20 along with the Adam optimizer [ 84] with 758 parameters b1=0.9,b2=0.999 ,e=108and a learning rate of lr=5105to perform 759 stochastic gradient descent. After 40 epochs of training, the learning rate decays by a 760 factor of 0.2. 761 To ensure that the trained model is able to generalize well, and does not overﬁt, 762 dropout with a 0.3 rate is used in the attention, fully-connected and embedding layers. 763 Furthermore, training and validation metrics (loss and accuracy) are monitored through- 764 out the training sequence, allowing us to halt training without overﬁtting. Figure 22 765 shows the training curves of experiment 1.1, highlighting that training for around 60 766 epochs is sufﬁcient in this case. 767DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 29 of 43 Figure 22. Training and validation accuracy for experiment 1.1 Beyond the 60th epoch no real gain in validation accuracy is made while training accuracy trends to 100%, which indicates overﬁtting. 6 Results 768 For each diagnostic experiment, a Transformer model is trained on the training 769 dataset, monitored on the validation dataset and ﬁnal evaluation takes place on the test 770 dataset. All model hyperparameters are static throughout the experimetnts, only the 771 number of training epochs varies, as we choose the epoch with the best validation per- 772 formance. We report in the following sections the results for all three sets of experiments 773 gathered on the previously un-seen test dataset. 774 6.1 Experiment set 1 775 The ﬁrst set of experiments aims to assess the potential of the Transformer for diag- 776 nostics given ideal conditions: a balanced data set where each individual degradation 777 class appears in equal measure. Table 7 reports the testing accuracy scores gathered for 778 experiments 1.1 and 1.2. We note a signiﬁcant difference between the results of the two 779 experiments, where the 4 class setup outperforms the 10 class setup by 30% in terms 780 of test accuracy. Knowing that Transformers are data-intensive models that respond 781 well to large datasets, this discrepancy can be partly explained by accounting for the 782 number of simulations per individual class. Indeed, as both training sets have the same 783 total number of samples, there are 2.5 times the amount of samples per severity class in 784 experiment 2.2. To factor out the disparity in the amount of data per class, we generate a 785 large version of training dataset of task 1.1, with a total of 10176 training samples and 786 2544 validation samples. In this experiment (1.1 Large), the number of training samples 787 per class is larger than for experiment 1.2 (1017 versus 840), yet the four class model 788 still outperforms the 10 class model by over 16%. In general it is advantageous to train 789 on larger amounts of severity classes from a risk management perspective, as it allows 790 for more optimal maintenance intervention schemes. However, as highlighted above, a 791 ﬁne-grained approach requires more training data in order to diagnose accurately with 792 high conﬁdence and could result in higher amounts of false positives, which could then 793 lead to an increase in Operation and Maintenance ( O&M) costs. Thus, there is a trade-off 794 between risk minimization via ﬁne-grained diagnostics of LEE severity and increased 795 costs due to higher amounts of false positives and requirements for larger training 796 datasets, which can be optimized to meet risk management policies and budgetary 797 speciﬁcations. 798DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 30 of 43 Table 7: Testing classiﬁcation accuracy of experiments 1.1 and 1.2. Testing Accuracy [%] Experiment All Predictions >70% Likelihood Exp. 1.1 66.67 72.65 Exp. 1.1 Large 80.63 85.56 Exp. 1.2 96.04 97.22 (a)Exp. 1.1: 10 levels of LEE severity. (b)Exp. 1.2: 4 levels of LEE severity. Figure 23. Confusion matrices of the predicted severity classes, evaluated on the test datasets of the ﬁrst set of experiments. The difference in accuracy can also be explained by the separability of the individual 799 severity classes. We show in Figure 23 the confusion matrices generated on the test 800 sets for experiments 1.1. and 1.2. It can be noted that the Transformer model is able to 801 successfully differentiate between the lowest and highest of erosion severities, as demon- 802 strated by the lack of misclassiﬁed samples with large differences in severity. Similar 803 LEE erosion classes exhibit higher amounts of wrongly classiﬁed samples. This aligns 804 with our expectations: small increases in the level of erosion only marginally increases 805 the roughness of the leading edge, which has a minor impact on blade aerodynamics. 806 It is therefore more challenging to distinguish between similar LEE severities based on 807 integrated pressure quantities. 808 6.2 Experiment set 2 809 Here, we assess the use of stochastic degradation paths as diagnostic training points. 810 Table 8 shows accuracy scores gathered on an un-seen degradation path for experiments 811 2.1 and 2.2. Figure 24 plots the predicted states of each MTS data point alongside the 812 true degradation path, for the test dataset which contains a full NHCPP path with up to 813 10 possible LE severity classes. In addition, we show the prediction conﬁdence of each 814 model output based on its softmax likelihood. Figure 25 displays the median predicted 815 severity using a three month rolling window. 816 These results show that this approach is viable and that we are able to successfully 817 detect the jumps in degradation severity with some limitations. The main issue with 818 training on stochastic degradation paths is the data imbalance. Indeed, the training 819 dataset contains a majority of points belonging to classes 0, 8 and 9, while classes 5, 6 and 820 7 are underrepresented. This is reﬂected in the results, where the accuracy is high for 821 class 0, but low for classes 5 and 6. Moreover, prediction conﬁdence is in general lower 822DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 31 of 43 Table 8: Testing classiﬁcation accuracy of experiments 2.1 and 2.2. Testing Accuracy [%] Experiment All Predictions >70% Likelihood Exp. 2.1 65.42 78.41 Exp. 2.2 67.75 69.48 for underrepresented classes. A possible solution to address this is to train on grouped 823 severity classes. We evaluate this approach in experiment 2.2 and plot the resulting 824 predictions in Figures 26 and 27. Overall, this approach yields higher prediction accuracy, 825 but it is still limited. Grouping artiﬁcially creates a somewhat balanced set and reduces 826 class sparsity, however, it also leads to increased intra-class heterogeneity. As a result, 827 learning to set the boundaries between the classes becomes more difﬁcult and dependant 828 on the number and the separability of the groups. This is reﬂected in the disparity of the 829 minimum 70% conﬁdence accuracy results between experiments 2.1 and 2.2. Indeed, 830 the low intra-class heterogeneity of the 10 class approach results in low conﬁdence 831 predictions for ambiguous samples, while the same ambiguous samples are incorrectly 832 classiﬁed with high a likelihood in the grouped approach due to the network accounting 833 for high intra-class variance. There is a balance to be found between the accuracy gain 834 and increase in intra-class heterogeneity induced by class grouping, depending on the 835 number of groups. Should grouping be desired, unsupervised clustering approaches 836 should be envisaged, in order to ﬁnd the optimal number of groups. 837 In experiments 2.1 and 2.2, non-physical predictions are penalised via the addi- 838 tional loss term (see equation 16). The result of using this extra component is apparent 839 in Figures 24 and 26: there are noticeably fewer misclassiﬁed points below the true 840 degradation path than above it (11.5% vs 23.08% for experiment 2.1, 6.67% vs 25.58% for 841 experiment 2.2). 842 Figure 24. Actual degradation path and predicted erosion severity classes for experiment 2.1. Only predictions with a likelihood greater than 70% are shown.DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 32 of 43 Figure 25. Actual degradation path and median predicted erosion severity classes for experiment 2.1 using a three month rolling window and upper rounding. Figure 26. Actual degradation path and predicted erosion severity classes for experiment 2.2. Only predictions with a likelihood greater than 70% are shown. Figure 27. Actual degradation path and median predicted erosion severity classes for experiment 2.2 using a three month rolling window and upper rounding.DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 33 of 43 6.3 Experiment set 3 843 In experiments 3.1 and 3.2, we assess how the diagnostics method performs on more 844 uncertain degradation paths, with multiple sources of variability. Not only are there 845 large class imbalances due to the stochastic degradation, but aerodynamic uncertainty 846 (see Figure 7) and long-term weather ﬂuctuations makes it challenging to distinguish 847 between different LEE severities. Weather variability is not commonly used for long-term 848 aero-elastic modelling in the wind energy industry, and is somewhat unrealistic, but 849 the goal is to understand the limitations of the Transformer model by making inference 850 extremely challenging. We report in Table 9 the accuracy scores gathered on an un-seen 851 degradation path for experiments 3.1 and 3.2. Overall, classiﬁcation accuracy is much 852 lower than previous experiments, but using only high conﬁdence predictions yields a 853 larger boost in accuracy. This can be explained by the inherent stochasticity and the high 854 intra-class variance of this dataset which makes high likelihood predictions rarer. 855 Table 9: Testing classiﬁcation accuracy of experiments 3.1 and 3.2. Testing Accuracy [%] Experiment All Predictions >70% Likelihood Exp. 3.1 35.00 45.75 Exp. 3.2 54.67 65.22 Figures 28 and 29 show how challenging this ﬁnal task is, where we see many 856 non-physical misclassiﬁed points and an overall low prediction conﬁdence. Here, the 857 data imbalance has a big impact on classiﬁcation performance. As the training set 858 does not contain any samples with a LEE severity class 9, the network is unable to 859 predict this severity on the given test set. Given these harsh training conditions, it 860 is unsurprising that performance is underwhelming. To mitigate this, class grouping 861 by LEE type is again tested in experiment 3.2 (see Figures 30 and 31). Compared to 862 experiment set 2, the beneﬁt of grouping classes is clear here: performance is improved 863 by almost 20%. This indicates that class grouping is a viable approach if intra-class 864 heterogeneity is high by default and some classes are completely underrepresented. 865 Nevertheless, even the results for grouped classes are lackluster, indicating that other 866 strategies should be examined to deal with the very imbalanced datasets induced by the 867 stochastic degradation paths. For instance, one could consider conditional retraining 868 methods. By modifying the architecture so that it outputs conﬁdence metrics based on 869 how similar a sample is to the training set, as implemented in [ 85], one could then launch 870 re-training procedures if multiple test samples that are very dissimilar to the training 871 space are encountered. 872DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 34 of 43 Figure 28. Actual degradation path and predicted erosion severity classes for experiment 3.1. Only predictions with a likelihood greater than 70% are shown. Figure 29. Actual degradation path and median predicted erosion severity classes for experiment 3.1 using a three month rolling window and upper rounding. Figure 30. Actual degradation path and predicted erosion severity classes for experiment 3.2. Only predictions with a likelihood greater than 70% are shown.DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 35 of 43 Figure 31. Actual degradation path and median predicted erosion severity classes for experiment 3.2 using a three month rolling window and upper rounding. 6.4 Transferability and Curriculum Learning 873 Here, we evaluate the transferability of the different experiments, and assess the 874 use of curriculum learning [ 86] as a method to improve accuracy on difﬁcult results. In 875 a curriculum learning training setup, the network is progressively exposed to harder 876 datasets. In our implementation, we ﬁrst train the model on the dataset from experiment 877 1.1, then progressively add the training data from experiments 2.1 and 3.1. The beneﬁts 878 of this approach should be twofold: (1) the data-intensive Transformer model is exposed 879 to a larger dataset, and (2) the weights of the model are ﬁne-tuned on more challenging 880 samples, which helps to avoid local minima during gradient descent. Tables 10 and 11 881 show the results of the different transferability and curriculum learning experiments. Table 10: Transferability and curriculum learning results. Test Accuracy [%] Training Set(s) Set 1.1 Set 2.1 Set 3.1 Set 1.1 – 71.25 28.25 Set 2.1 65.00 – 29.42 Set 3.1 30.00 39.25 – Set 1.1, then 2.1 77.71 83.50 20.33 Set 1.1, then 2.1 & 3.1 77.08 80.50 37.08 882 Table 11: Transferability and curriculum learning results, for high conﬁdence predictions. >70% Likelihood Test Accuracy [%] Training Set(s) Set 1.1 Set 2.1 Set 3.1 Set 1.1 – 80.21 30.19 Set 2.1 77.11 – 32.86 Set 3.1 35.58 40.97 – Set 1.1, then 2.1 81.35 87.29 20.52 Set 1.1, then 2.1 & 3.1 84.42 88.53 50.22 We remark the following: 883DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 36 of 43 • The high accuracy of transferability tests between experiments 1.1 and 2.1 highlight 884 the similarity of these datasets. Moreover, the balanced dataset of experiment 1.1 885 improves prediction quality for set 2.1. 886 • Although both sets 1.1 and 2.1 contain all degradation classes, networks trained on 887 these sets do not transfer well to test set 3.1, undoubtedly due to the uncertainty 888 included in this experiment. 889 • There is better transferability between sets 2.1 and 3.1, than between sets 1.1 and 3.1. 890 This can be explained by the extra loss component which enforces monotonicity in 891 the second and third sets of experiments. 892 • The ﬁrst curriculum learning experiment (set 1.1 then 2.1) yields a large performance 893 boost for test sets 1.1 and 2.1. This is to be expected due to the data-intensive nature 894 of Transformers and the similarity between these sets. 895 • The ﬁrst curriculum learning experiment reduces transferability to set 3.1. This is 896 an indication of a loss of capacity to generalize to ambiguous data. 897 • The best test results on set 3.1 are obtained in the second curriculum learning 898 experiment. This suggests that pretraining on sets with reduced stochasticity 899 followed by ﬁne-tuning on uncertain data is an effective approach. 900 • The second curriculum learning experiment yields the best high conﬁdence accuracy 901 for sets 1.1 and 2.1. Thus, adding difﬁcult, stochastic data points of experiment 3.1 902 to the training dataset helps with regularization, enabling the model to construct a 903 better internal representation of each severity class. 904 The curriculum learning test which includes the data from all experiments, leads to 905 the highest high-conﬁdence accuracy for the 2.1 test set. As this set is representative of a 906 degradation path with standard aero-elastic simulation practices, we plot the results in 907 Figures 32 and 33. 908 Figure 32. Results for test dataset 2.1 showing actual degradation path and predicted erosion severity classes for the curriculum learning experiment combining all datasets. Only predictions with a likelihood greater than 70% are shown.DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 37 of 43 Figure 33. Results for test dataset 2.1 showing actual degradation path and median predicted erosion severity classes for the curriculum learning experiment combining all datasets, using a three month rolling window and upper rounding. 7 Limitations and Discussion 909 All assumptions and modelling choices, constants, choice of probability distribu- 910 tions and their parameters in the proposed NHCPP model for leading edge erosion 911 degradation could in principle be calibrated [ 87,88]. This, however, will prove very 912 challenging given the data-scarce nature of measurements of erosion processes on op- 913 erating wind turbines in the ﬁeld over long time horizons. In addition, the impact of 914 ﬁeld data quantity and quality on model calibration and its implication on the generated 915 NHCPP will need to be closely evaluated. One interesting extension to our NHCPP 916 macro degradation model, is to further couple it to a micro model that describes the 917 relative velocity between particles (including their density and size) and the blade, which 918 leads to an expression for the mass ﬂow of particles impinging the blade leading edge 919 per unit of area [ 10]. Subsequently, a matricial expression for the number of particles 920 impacting a unit surface is derived to describe the dependency of erosion severity from 921 the atmospheric conditions at the site and the operational settings of the wind turbine. 922 Similar to the NHCPP , the Gamma process has also been traditionally used for 923 modeling deterioration. A Gamma process has independent non-negative increments, 924 but the shocks follow a Gamma distribution. While in a Poisson process there is a ﬁnite 925 number of jumps in a ﬁnite time interval, in a Gamma process there are inﬁnite numbers 926 of jumps in a ﬁnite time interval. For this reason, the former are suitable for modelling 927 usage such as damage due to sporadic shocks, whereas the Gamma processes are suitable 928 for describing continuous (’graceful’ - a succession of tiny increments), monotonic and 929 gradual deterioration. The Gamma process can be regarded as a compound Poisson 930 process of Gamma-distributed increments in which the Poisson rate tends to inﬁnity and 931 increment sizes tend to zero in proportion [89]. 932 One of the main outcomes of the diagnostic experiments is the fact that the Trans- 933 former model is data-intensive. Indeed, our results point to the fact that more data almost 934 always aids in classiﬁcation performance. Furthermore, training on balanced datasets 935 is always desirable, but obtaining these is difﬁcult in real-life monitoring scenarios as 936 it would require gathering data from many systems undergoing different degradation 937 paths. To overcome this difﬁculty, the coupled NHCPP aeroleastic simulation setup 938 can be used to replicate realistic degradation and produce balanced training sets which 939DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 38 of 43 can then be supplemented with real ﬁeld data, in a similar manner to the curriculum 940 learning experiments. 941 As mentioned in the literature review, feature-extraction learning methods are often 942 preferred for MTS classiﬁcation tasks owing to their higher precision. In this context, 943 We tested the use of the state-of-the-art MINIROCKET method [ 29] with 10k kernels, 944 combined with a standard linear classiﬁer at the output for our diagnostics problem. 945 This yielded an accuracy of only 55.21% on the test set of experiment 1.1 and 68.96% 946 on the test set of experiment 1.2. This highlights the fact that very long sequences are 947 challenging to deal with, while many methods found in the literature are not optimized 948 for such tasks. Furthermore, feature extraction methods may encounter difﬁculty ﬁnding 949 optimal features, due to the fact that this problem features a very slow degradation 950 rate with no ﬂagrant traits from one stochastic sample to another. On the contrary, 951 Transformers are well suited for this kind of input, as the self-attention mechanisms are 952 effective at ﬁltering out unnecessary data while focusing on the parts of the signal that 953 are important. 954 While the loss term added to penalize non-monotonicity of the degradation process 955 does help to reduce non-physical predictions, it is not strictly enforced. Another option 956 would be to enforce this condition in the softmax output layer. Although this method 957 would most certainly improve accuracy, it would lead to error accumulation when 958 implemented in the context of a continuous monitoring setup. With our approach 959 erroneous predictions can be ignored by using a rolling window for instance. The 960 proposed model could, however, be further improved by additional tuning of the 961 objective function. One such improvement involves incorporation of a inﬂow-dependant 962 multiplier that reduces the importance of samples that have highly unlikely inﬂow 963 conditions. Furthermore, while the leading edge degradation process is inherently time- 964 dependant, time was not factored into the ﬁnal diagnostics approach. Initial tests which 965 included a global time stamp in the model input, showed the Transformer overﬁtting 966 to this variable. A probabilistic loss term that penalizes unlikely severities given the 967 amount of operational time, could be a potential solution to this. 968 Another possible improvement could be the modiﬁcation of the training procedure 969 to a bi-directional procedure, similar to the BERT model [ 83]. In this approach, a pairwise 970 input would be fed into the Transformer, and during the masked pre-training phase, the 971 model would try to predict the class of either the ﬁrst or second MTS at random. This 972 could lead to a better representation of the data and allow the use of known reference 973 states for comparative diagnostics. 974 Our proposed approach does not take positional information into account, yet, as 975 shown in Figure 20, this may allow to more clearly distinguish the LEE severity classes 976 along the span of the blade. Inserting positional information into our Transformer model 977 can be achieved in different ways; one possible approach is to encode a graph structure 978 in the Transformer architecture, as implemented in [90]. 979 8 Conclusion 980 In this work we tackle the problem of leading edge erosion of wind turbine blades 981 on a two-fold front. Firstly, we deal with the development of an appropriate model 982 for simulation of such degradation processes and, secondly, we propose a monitoring- 983 driven method for diagnostics of such damage processes. On the ﬁrst front, we propose 984 a stochastic spatio-temporal empirical model for modelling leading edge erosion degra- 985 dation based on Non-Homogeneous Compound Poisson Process and couple this to the 986 non-linear time marching OpenFast aeroelastic wind turbine computer simulator. The 987DRAFT Version September 1, 2021 submitted to Journal Not Speciﬁed 39 of 43 coupled model is able to compute the aeroelastic non-stationary dynamic response of a 988 wind turbine reﬂecting its behaviour under the effect of leading edge erosion, varying 989 inﬂow conditions and aerodynamic uncertainties over a long period of degradation time 990 horizon. On the diagnostic front, we adapt a deep neural network, namely a Transformer, 991 to allow for use on very long sequence multivariate time-series. This allows us to solve 992 the problem of spatio-temporal diagnostics of leading edge erosion on wind turbine 993 blades, using the data emanating from the Non-Homogeneous Compound Poisson 994 Process setup; a scheme that intends to emulate data recorded by aeroacoustic sensors 995 placed on wind turbine blades. We showed that the diagnostics model effectively cap- 996 tures the temporal trends induced by long-term degradation of the leading edge. An 997 attractive feature of this method is that it is well-suited for spatio-temporal degradation 998 problems with a very long time horizon. 999 Acknowledgments 1000 This work is funded by the BRIDGE Discovery Program of the Swiss National 1001 Science Foundation and Innosuisse, project number 40B20_187087 . The authors also 1002 acknowledge the support of the European Research Council via the ERC Starting Grant 1003 WINDMIL (ERC-2015-StG #679843) on the topic of Smart Monitoring, Inspection and 1004 Life-Cycle Assessment of Wind Turbines. 1005