Computationally Efficient Uncertainty Minimization in Wind Turbine Extreme Load Assessments

During the last few decades, significant efforts have been made to scale up wind turbine sizes in order to harvest more energy from wind. The substantial increases in component size, weight, and length, however, add to the loads imposed on a turbine's structure [1], which cause degradation and failure. The key for achieving reliable operations during a turbine's desired lifetime is a quantitative estimation of the expected loads at the design stage, so that appropriate design parameters can be chosen to guarantee a system's reliability level as well as to reduce the cost of manufacturing. Load underestimation can lead to premature failures of major functional components, which in turn lead to high operations and maintenance costs [2–4], whereas load overestimation can result in overly conservative design.

Evaluating loads at the design stage typically uses physics-based load simulators and short-term field testing data. The short-term field data collected during the design stage are limited, unlike operational data available during service time [5,6]. Therefore, it has become common to use aeroelastic simulation tools to mitigate the data scarcity issue. The U.S. Department of Energy's National Renewable Energy Laboratory (NREL) has developed aeroelastic wind turbine simulators to help designers estimate loads on the turbine components [7,8].

The International Electrotechnical Commission (IEC) requires estimating the loads associated with important design load cases (DLCs) specified in the design standard, IEC 61400-1 [9]. Among the DLCs specified, DLC 1.1 requires assessing the extreme load, the largest load level that a turbine can resist during normal operations [9]. Two major approaches for estimating the extreme load with a wind turbine simulator are the Monte Carlo sampling (MCS) and the load extrapolation [10–13]. Taking the first approach, Moriarty [14] used grid computing with 60 desktops to simulate a turbine's load responses during 5 years of operations. Recently, Manuel et al. [15] simulated almost 100 years of operations with a Linux computer cluster with 1024 cores.

Although the advancement in computing power relieves the computational burden of MCS to some extent, the large number of simulation runs required in MCS is still computationally expensive. Moreover, as simulation models become more complicated to represent complex stochastic systems, each simulation run becomes computationally more demanding, which could make MCS an infeasible option in practice. Another disadvantage of MCS is that the results are highly uncertain, implying that the estimation from a set of simulation runs could largely differ across other sets.

Considering the computational limitations of MCS, extrapolating the outcomes from a limited number of simulation runs has been suggested in the literature. The basic concept of the load extrapolation is that it predicts the trajectory toward the target load exceedance probability (or probability of exceedance (POE)), using the simulation outputs. Several studies use statistical extrapolations in the binning framework [10–13]. The binning method divides the wind condition into a finite number of bins, e.g., 3–4 m/s, 4–5 m/s, $…$⁠, 24–25 m/s for wind speeds. Then, it obtains the load outputs from a small number of simulations in each bin and empirically approximates the conditional load distribution using the data belonging to each bin. The POE is estimated by summing the approximated conditional distribution weighted by the probability of each bin.

In the binning method, the statistical extrapolation has been made by approximating the load distribution with a prespecified distribution (e.g., generalized extreme value (GEV) distribution), which is called a reduced model, surrogate model, or metamodel [16]. This approach does not require a significant amount of data as MCS does. However, it does not guarantee the “unbiasedness” of estimation because the estimation is made on the metamodel and not on the simulation model itself. As discussed in Ref. [5], this metamodel-based approach is useful when the extreme load is estimated using field testing data at the design stage and the available data size is limited. When the extreme load is estimated using a simulator, however, it is preferable to estimate the load directly on the simulation model and not on the metamodel that approximates the simulation model, to obtain an unbiased estimation.

The challenges lie in how to estimate the extreme load in a way that is both computationally efficient and feasible without relying on the load extrapolation and in how to minimize the estimation uncertainty. We tackle these challenges based on the stochastic importance sampling (SIS) methods proposed in Ref. [17]. Wind turbine loads highly depend on stochastic operating conditions, and not all weather conditions need to be treated equally [18]. Wind turbines may experience greater loads in certain conditions. For example, in-plane blade tip deflections tend to increase as wind speed increases, whereas flapwise bending moments at blade roots tend to peak around the rated wind speed, say 10–15 m/s. Based on such observations, the SIS methods in Ref. [17] guide the simulation efforts to focus on the critical weather conditions, so that high load events will occur more frequently.

Our approach to assess the extreme load is to integrate the theory of order statistics with the SIS methods. The novelty of the proposed approach mainly lies in that we apply SIS to the load estimation problem, as opposed to the POE estimation problem. We use the order statistics to make a theoretically sound and practical connection between the SIS-based POE estimation and the load estimation, enabling the transfer of the variance reduction in the POE estimation to the load estimation. Our implementation results using the NREL's simulator demonstrate the superiority of the proposed approach over alternative methods in three aspects: (a) the proposed approach can substantially reduce the estimation variations over MCS; (b) it can estimate large extreme load levels with a much smaller computational resource, compared to MCS; and (c) it provides an unbiased POE estimation unlike the load extrapolation approach.

This paper is organized as follows: Section 2 describes the problem. Section 3 presents the proposed methodology. Sections 4 and 5 discuss the implementation details and results, respectively. Section 6 concludes and suggests future research.

where l_T is the resistance level. Following the wind industry practice and the international standard, IEC 61400-1 [9], Y denotes the largest load response observed in a 10 min time series (e.g., 10 min maximum flapwise bending moment at a blade root), and X is a 10 min average wind speed used as an input to the simulator. The wind speed, X, lies within a specified range of the cut-in speed, x_in m/s, to the cut-out speed, x_out m/s. $f(·)$ denotes the probability density function of the 10 min average wind speed, X. We use a truncated Rayleigh distribution for X (to be detailed in Sec. 4) [14,17]. When multiple environmental factors such as turbulence intensity and significant wave height for offshore turbines are included in Eq. (2), X becomes a multivariate input vector, and $f(·)$ denotes a joint density function of the input vector [19].

where “$inf$” implies the infimum.

Note that the extreme load level depends on the target POE, P_T. When P_T is small, a large extreme quantile, l_T, needs to be estimated, and vice versa. When a wind turbine is designed to resist the resulting load level, l_T, the anticipated load exceedance probability is smaller than, or equal to, the target level, P_T; this is why l_T is also called a design load [20]. The extreme load considered in this study is the largest load that a turbine can withstand during operations between the cut-in and cut-out wind speeds, as specified in DLC 1.1 of IEC 61400-1 [9]. Loads occurring during parked conditions, e.g., hurricanes, are not considered in this study.

This section presents the methods that provide the efficient extreme load assessment without relying on computationally intensive MCS or statistical extrapolation.

For the sake of readability, we first summarize SIS methods (more details are available in Ref. [17]). Unlike MCS that samples wind speeds from the original density, f, in Eq. (2), SIS methods sample wind speeds from a biased density, q, so that large load outputs can be observed more frequently.

Two SIS methods, called SIS1 and SIS2, are proposed in Ref. [17]. From a biased density, q, SIS1 independently samples M different wind speeds, x_i, $i=1,…,M$⁠. Then, it runs the simulator N_i times at each sampled x_i, $i=1,…,M$⁠, where N_i is called an allocation size. Figure 1 depicts the overall procedure of SIS1. Let $Yj(i), i=1,…,M, j=1,…,Ni$⁠, denote the load output obtained in the $jth$ run at the input wind speed, x_i. With these outputs, the POE estimate of SIS1 is

where $I(·)$ is an indicator function. The total number of simulation runs is $NT=∑i=1MNi$⁠. Note that at each sampled x_i, the N_i load outputs, $Yj(i), j=1,…,Ni$⁠, are used to estimate the conditional POE, $P (Y>l∣X=xi)$ in Eq. (4). Because the input wind speed is sampled from the biased density, q (rather than the original density, f), SIS1 uses the likelihood ratio, $f(xi)/q(xi)$ in Eqs. (4) and (5), to make the estimator unbiased.

where s(x) is $P(Y>l∣X=x)$ and $Cq1$ is $∫Xff (x)(1/NT)s (x)(1−s (x))+s (x)2 dx$ with $Xf$ denoting the support of f. Note that $qSIS1$ focuses sampling efforts on the important x's that correspond to the large conditional exceedance probability, s(x). For the sampled $xi,i=1,…,M$⁠, SIS1 allocates the larger N_i to those x_i's associated with the relatively smaller $s(xi)$ in order to balance the simulation efforts within the important region (for more details, we refer the reader to Ref. [17]).

where Y_i is an output at x_i, $i=1,…,NT$⁠. Because the allocation size is one at each sampled wind speed, N_T wind speeds are sampled in SIS2, given the total number of simulation runs, N_T. Similar to SIS1, $f(xi)/q(xi)$ in Eq. (8) is used to maintain the unbiasedness of the POE estimator in SIS2.

where $Cq2$ is $∫Xff (x)s (x)dx$⁠.

where SIS collectively denotes SIS1 and SIS2. Below, we explain how to estimate the extreme quantile using the order statistics.

We use the N_T load outputs (i.e., N_T 10 min maximum load responses) from N_T simulation runs. First, we sort the N_T outputs and consider the order statistics of the outputs, $Y(1),…,Y(NT)$⁠, where $Y(k)$ denotes the kth smallest load output. Let $x(k)$ and $N(k)$ be the corresponding wind speed and allocation size, respectively, at which $Y(k)$ is simulated, $k=1,…,NT$⁠. We substitute $Y(k)$ for l in Eq. (10) for $k=1,…,NT$ to obtain the POE estimate at each output.

because no output exceeds the largest output, $Y(NT)$⁠, and therefore, $I(Yj(i)>Y(NT))=0$ for $∀i,∀j$⁠.

The inequality in Eq. (12) is satisfied when l is equal to, or larger than, $Y(NT−1)$⁠. Consequently, $Y(NT−1)$ becomes the extreme load estimate, $l̂PT$⁠, because it is the infimum of l that satisfies the inequality. In fact, $Y(NT−1)$ is the largest extreme load estimate associated with a nonzero P_T, which can be estimated using N_T outputs.

In Eq. (13), we sort the N_T outputs, $Yj(i), i=1,…,M$⁠, and $j=1,…,Ni$⁠, from the smallest to the largest, $Y(a), a=1,…,NT$⁠. $N(a)$ is the allocation size at $x(a)$⁠, where $Y(a)$ is obtained. In Eq. (14), the $(NT−k)$ outputs greater than $Y(k)$⁠, i.e., $Y(k+1),…,Y(NT)$⁠, are selected to compute the POE estimate. As a result, $Y(k)$ becomes the extreme load estimate at the POE level of $PT=(1/M)∑a>kf(x(a))/(N(a)·qSIS1(x(a)))$⁠.

The extreme load estimate with SIS2 can be similarly obtained using SIS2's POE estimate in Eq. (8). Again, by sorting the N_T outputs, we obtain the order statistics of the outputs, $Y(k)$⁠, and the corresponding wind speeds, $x(k)$⁠, for $k=1,…,NT$⁠.

Consequently, $Y(NT−1)$ is the estimated extreme load at the POE level of $f(x(NT))/(NT·q(x(NT)))$⁠.

Consequently, MCS can obtain the largest extreme load estimate at the POE level of $1/NT$⁠. Similarly, the kth smallest load output, $Y(k)$⁠, becomes the estimated extreme load at $PT=(NT−k)/NT$ because $(NT−k)$ outputs exceed $Y(k)$⁠.

Next, we explain how the proposed approach can reduce the computational burden compared to MCS. In MCS, the smallest POE level that MCS can estimate is $1/NT$⁠. On the contrary, in SIS1 and SIS2, the smallest POE level associated with the second largest output, $Y(NT−1)$⁠, is $f(x(NT))/(M·N(NT)·qSIS1(x(NT)))$ and $f(x(NT))/(NT·qSIS2(x(NT)))$ in Eqs. (11) and (15), respectively. The SIS methods control the sampling process to focus on the critical wind condition, where q (i.e., $qSIS1$ or $qSIS2$⁠) is high. When f and q are significantly different, the likelihood ratio, f/q, is mostly smaller than 1 in the sampled speed. This can lead to $f(x(NT))/(M·N(NT)·qSIS1(x(NT)))<1/NT$ and $f(x(NT))/(NT·qSIS2(x(NT)))<1/NT$⁠. In other words, the smallest POE level that SIS1 and SIS2 can estimate is less than that of MCS, i.e., $P̂SIS (Y>Y(NT−1))<P̂MCS (Y>Y(NT−1))$⁠. As such, with the same computational effort, the SIS methods can estimate the larger extreme load associated with a smaller POE than that estimated by MCS.

Moreover, SIS1 and SIS2 reduce uncertainties in extreme load estimation because they minimize the estimator variances. Section 5 presents how much computational efforts and uncertainties can be reduced by employing the proposed approach.

We specify the key details in implementing the proposed approach and summarize the implementation procedure.

where f_R denotes the untruncated Rayleigh density function, $fR(x)=xe−x2/2τ2/τ2$⁠, x > 0 with $τ=2/π·10$ (unit: m/s), and $FR(x)=1−e−x2/2τ2$ is the corresponding cumulative distribution function [14]. We use x_in = 3 m/s and x_out = 25 m/s as the cut-in and cut-out wind speed, respectively. We use the wind regime corresponding to a class I-B in the IEC Standard [9]. Also, the NREL 5 MW baseline turbine model specification [24] is employed in our simulations.

In $qSIS1$⁠, N_i, $i=1,…,M$⁠, and $qSIS2$ in Eqs. (6), (7), and (9), the conditional POE, s(x), is unknown in practice, so it needs to be approximated. As we proposed in Ref. [17], we get a pilot sample consisting of 250 load responses over a range of wind speeds between 3 m/s and 25 m/s. Then, we fit the nonhomogeneous GEV distribution to the load responses to approximate the conditional POE (for the detailed procedure, see Ref. [17]). Note that we use the approximated conditional POE only to guide the sampling efforts for the input, x, and not to estimate the load output.

We also need to prespecify the resistance level, l, in s(x). This resistance level, l, is necessary only for controlling the sampling from q and not for deciding the extreme load level. One desirable observation is that $supp(ql2)⊆supp(ql1)$ for $l2>l1$⁠, where q_l is either $qSIS1$ in Eq. (6) or $qSIS2$ in Eq. (9) with the quantile l, and $supp(q)={x∈ℝ:q(x)>0}$⁠. This implies that the support of SIS density for a higher quantile is included in the support of SIS density for a lower quantile, because the zero-probability failure event region, ${x∈Xf:P(Y>l|X=x)=0}$⁠, increases or remains the same, as l increases. Therefore, when we estimate a higher quantile than l used for obtaining the SIS density, q_SIS, the POE estimations remain unbiased. We use l = 2.33 m and l = 14,600 kN · m for the in-plane tip deflection and flapwise bending moment, respectively, in our implementation.

Figure 2 compares the original input density, f, (i.e., truncated Rayleigh density) and the SIS1 density, $qSIS1$⁠, for the two types of load responses: in-plane tip deflections and flapwise bending moments. The SIS2 densities, $qSIS2$⁠, are similar to $qSIS1$⁠. Since the in-plane tip deflections tend to increase as wind speed increases, in Fig. 2(a), SIS1 reweights the sampling efforts to focus on the input conditions where large in-plane tip deflections are likely to occur. Likewise, in Fig. 2(b), the density of $qSIS1$ is high around 10–15 m/s where large flapwise bending moments are observed.

We summarize the procedure for estimating the extreme load. To quantify the variability of the resulting extreme load estimates in each method, we repeat step 2–4 in the following procedure multiple times:

Step 1: Take a pilot sample and obtain $qSIS1$ in Eq. (6), N_i, $i=1,…,M$⁠, in Eq. (7) and $qSIS2$ in Eq. (9).

Step 2: In SIS1, sample M wind conditions, $x1,…,xM$⁠. At each x_i, run the simulator N_i times to obtain $Yj(i), i=1,…,M, j=1,…,Ni$ (Likewise in SIS2 with N_i = 1, $i=1,…,NT$⁠).

Step 3: Sort the N_T outputs and obtain the order statistics, $Y(1),…,Y(NT)$⁠.

Step 4: Find the POE level at each $Y(k), k=1,…,NT$⁠, using Eqs. (14) and (16) in SIS1 and SIS2, respectively.

In step 1, M should be predetermined in SIS1. Our prior study [17] indicates that the performance of SIS1 is generally insensitive to the ratio of M to N_T. We use M = 500 and $NT=3000$⁠.

For the in-plane tip deflection, we obtain the order statistics, $Y(1),…,Y(NT)$⁠, from $NT=3000$ outputs in each method. Figure 3 compares the results from 16 repetitions of each method. In each repetition, we linearly interpolate the results [16], so Fig. 3 presents the 16 trajectories of each method. We observe some notable features of the proposed approach:

• SIS consistently provides lower standard errors of quantile estimates (i.e., narrower band width). For example, when $PT=1/3000$ in Fig. 3(a), the extreme load estimates of MCS range from 2.27 m to 2.56 m, whereas SIS1 results in a much narrower range, 2.43–2.48 m. The wider range implies high uncertainties, because each repetition generates significantly different results. By comparing the ranges, we can see that compared to MCS, SIS substantially reduces the uncertainty and provides more accurate estimations. Moreover, the range increases much faster in MCS as P_T decreases. Therefore, MCS would produce unreasonably uncertain estimations when POE level is less than $10−4$⁠.

• SIS can estimate much smaller POEs than MCS with the same computational resource, $NT=3000$⁠. Among the 16 trajectories, the median POE level that SIS1 (SIS2) can reach is around $1.1×10−5 (7.7×10−6)$⁠, as seen in Fig. 3(b). This is much smaller than the POE level of $1/3000=3.3×10−4$ which MCS can estimate with the same computational resource. To reach the same POE levels as SIS1 and SIS2 with 3000 simulation runs, MCS requires around 91,000–130,000 runs (from $(1.1×10−5)−1 to(7.7×10−6)−1$⁠).

• The results of SIS2 are similar to those of SIS1, as shown in Fig. 3(b).

For the flapwise bending moment, we conduct more repetitions, namely, 20 repetitions because the flapwise bending moment shows larger variations than the in-plane tip deflection [23]. Figure 4 shows the results. We summarize our observations:

• The proposed approach provides higher accuracy than MCS. When $PT=1/3000$ in Fig. 4(a), the extreme load estimates of MCS range from $1.47×104$ kN · m to $1.56×104$ kN · m, whereas SIS1 results in a narrower range, $1.504×104$ kN · m to $1.55×104$ kN · m. The uncertainty of MCS increases more rapidly than that of SIS1 as the target POE decreases.

• The performance gain in terms of the computational efficiency appears less remarkable compared to the in-plane tip deflection case. The reason is that the SIS density, q_SIS, for the flapwise bending moment is closer to the original density, f, as seen in Fig. 2. Therefore, in the critical zone, the ratio, f/q, for the flapwise bending moment is not as small as that obtained for the in-plane tip deflection. However, overall SIS1 and SIS2 estimate the extreme load levels associated with smaller POEs than MCS, given the same computational effort. To simulate the same POE levels as SIS1 and SIS2 with 3000 runs, MCS requires around 9000–11,000 runs (from $(1.1×10−4)−1 to (9.5×10−5)−1$⁠) (see Fig. 4(b)).

Note that the proposed approach is separately implemented for each load response, whereas sampling from MCS does not depend on the load type. Therefore, it is interesting to see the advantage of the proposed approach for multiple types of extreme load estimations. Recall that we use $NT=3000$ for both in-plane tip deflection and flapwise moment. When MCS uses $NT=6000$ replications, the largest extreme load level that MCS can estimate is associated with $PT=1/6000=1.67×10−4$⁠. This POE is still larger than the POE level that SIS can achieve for both load types. More importantly, considering the large variability of the MCS estimation, the SIS methods are advantageous even though they use a separate implementation for each load type.

Devising the method that can jointly estimate the extreme loads of multiple types will be the subject of our future study. For example, a mixture of the current two SIS densities for in-plane tip deflections and flapwise bending moments will be likely bimodal, focusing sampling efforts on around the rated and cut-out wind speeds. Using such a bimodal SIS density will be more advantageous than using the original input density.

for SIS1 and SIS2, respectively. Here, N_R is the number of repetitions; x_ri and N_ri are the ith sampled input and ith allocation size, respectively, in the rth repetition; and $Yj(ri)$ is the jth output at x_ri for $j=1,…,Nri$⁠. With the POE estimates in Eqs. (17) and (18), we can estimate the extreme loads in the same way as described in Secs. 3.2.1 and 3.2.2 for SIS1 and SIS2, respectively.

For the in-plane tip deflection, Fig. 5 shows that combining 16 repetitions of MCS results with N_T = 3000 leads to the largest extreme load estimate of 2.67 m corresponding to the POE of $1/48,000=$ 2.1 × 10⁻⁵. With the same computational efforts, SIS1 and SIS2 lead to the larger extreme load estimates of 3.01 m and 2.98 m, corresponding to the smaller POEs of 7.6 × 10⁻⁸ (250-year return period) and 3.5 × 10⁻⁷ (54-year return period), respectively. In particular, we obtain the extreme load estimates with 50-year return period as 2.96 m from SIS1 and 2.97 m from SIS2. Considering that MCS needs 2.6 × 10⁶ replications to obtain a 50-year extreme load estimate, this result implies that SIS1 (or SIS2) with 48,000 replications saves the computational effort over MCS by at least 98%.

For the flapwise bending moment, the largest extreme load estimates from MCS, SIS1, and SIS2 with combined repetitions are 15,770 kN · m, 15,830 kN · m, and 15,950 kN · m, corresponding to the POEs of 1.7 × 10⁻⁶ (⁠$1/60,000$⁠), 4.8 × 10⁻⁶, and 6.2 × 10⁻⁶, respectively. The largest extreme load estimates we obtain from SIS1 and SIS2 correspond to a return period of 3–4 years.

It should be noted that these extreme load estimates are subject to sample-to-sample variations. To more accurately estimate the expected extreme load level during a turbine's design life (e.g., 50 years), we need multiple extreme load estimates and compute the average of the estimates. Due to the limitation of our computational resources, we cannot obtain the multiple extreme load levels associated with the 50-year design life in this study. However, our results show the promise of the proposed approach in obtaining the average extreme load levels during a turbine's design life.

We also compare our estimation results with the extreme load levels reported in Ref. [23] that uses much more computational resources with a supercomputer. We, however, note some differences in simulation settings. Barone et al. [23] used the FAST version 7.00.01 a-bjj (released on Nov. 5, 2010) and the inflow turbulence generated on a 20 × 20 square grid with width of 137 m, whereas we use the FAST version v7.01.00 a-bjj (released on Feb. 16, 2012) with 15 × 15 grid.

For the in-plane tip deflection, the extreme load estimates with a return period of 50 years from our approaches, SIS1 and SIS2, are around 3.0 m (see Fig. 5). In Ref. [23], the estimate with the same return period is between 3.1 m and 3.3 m. The difference between our estimates and that in Ref. [23] can be regarded small, considering the differences in the simulation settings and the inherent uncertainties in stochastic simulation.

For the flapwise bending moment, we obtain the extreme load estimates with a return period of 3–4 years, which are around $1.6×104$ kN · m. The result in Ref. [23] suggests that the estimate with the similar return period is $1.7×104$ kN · m. The difference in the order of $0.1×104$ kN · m can be regarded small due to the same reasons mentioned above and the large variability in the flapwise bending moment.

where B is the number of bins. We use the equal bin width, 1 m/s, for each bin.

where the approximation, $P̂(Y>lT∣xb−1<X≤xb)$⁠, is commonly achieved by fitting a parametric distribution. Although the nonparametric estimation of the conditional probability is possible, it requires a large number of simulations for each bin to estimate the extreme tail distribution. Therefore, the parametric approach is more widely employed in practice.

In our implementation, we sample 200 wind speeds at each bin and run the simulators to obtain the corresponding 200 load responses. To obtain $P̂(Y>lT∣xb−1<X≤xb)$⁠, we fit the parametric distribution by the maximum likelihood estimation. From our extensive implementations with different parametric distributions including GEV (note that Gumbel distribution is a special type of GEV distribution), Weibull, Gamma, and Lognormal distributions, the GEV distribution provides the best fit for both in-plane tip deflections and flapwise bending moments.

Figure 6 shows the trajectory from the binning method. The trajectory is obtained by computing $P̂BIN$ in Eq. (20) by incrementally increasing the target quantile, l_T, by 1 throughout the range displayed in Fig. 6. For comparison, Fig. 6 also shows the trajectories of SIS1 and SIS2. Recalling that SIS1 and SIS2 yield unbiased POE estimation, the results suggest that the binning method tends to overestimate the extreme loads for both in-plane tip deflection and flapwise moments. As the target POE decreases, the slope of SIS trajectory tends to decrease, whereas the binning method fails to capture this pattern around the tail and the resulting bias grows larger. In particular, the large bias from the binning method for the flapwise moment indicates that the load extrapolation based on the binning method can be highly unreliable, especially when the method is applied to the data involving the large variability as in the flapwise bending moment.

As a remark, the common form of binning method uses the uniform distribution for sampling in each bin, mainly out of convenience and not for reducing the estimation uncertainty. Such approach using the uniform sampling is in fact similar to using the uniform density function for the SIS density, q(x), in Eq. (5) or (8), if we use the corresponding POE estimator instead of the binning estimator in Eq. (20). However, the resulting estimator of the extreme load cannot be as useful as ours, because the importance of each input condition is not properly reflected by the uniform distribution.

We also note that under the binning framework, bootstrapping has been employed to understand the uncertainty in the design load estimation [25]. However, the bootstrap technique cannot reduce the uncertainty itself as the existing data are resampled, but only quantifies the estimation uncertainty. Therefore, to reduce the variability of the estimator, we need a sophisticated approach to guide the data sampling before the data are observed.

This study proposes a new approach to efficiently assess a wind turbine's extreme loads using stochastic simulations. Randomness of wind turbulence plays one of the most significant roles in causing the variations in the extreme load level, making extrapolations widely inaccurate. Our sampling-based approach using SIS fully considers the uncertainties lying in the wind turbulence by avoiding the load extrapolation. Both MCS and SIS repeat running the simulation multiple times with different random seeds to observe the load variations caused by the stochastic wind turbulence. However, in contrast to MCS, SIS samples important wind speeds more often, helping us observe more extreme values in the load response. Therefore, as we have more information on the tail distribution of the load, the proposed approach is able to estimate the extreme load more accurately.

In the future, we will investigate the optimal way to design the simulation experiment for simultaneously assessing the extreme load responses of multiple types to further reduce the computational effort. We are also interested in investigating the effects of multiple environmental factors on wind turbine loads, in order to help develop reliable designs for turbines operating in various environmental conditions.

We would like to thank the Editor, Associate Editor, and anonymous reviewers for their constructive comments on various aspects of this work. This work was partially supported by the National Science Foundation (Grant No. CMMI-1362513).