To harvest more energy from wind, wind turbine size has rapidly increased, entailing the serious concern on the reliability of the wind turbine. Accordingly, the international standard requires turbine designers to estimate the extreme load that could be imposed on a turbine during normal operations. At the design stage, physics-based load simulations can be used for this purpose. However, simulating the extreme load associated with a small load exceedance probability is computationally prohibitive. In this study, we propose using importance sampling combined with order statistics to reduce the computational burden significantly while achieving much better estimation accuracy than existing methods.

## Introduction

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 [1013]. 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 [1013]. 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.

## Problem Description

Let Y denote the maximum load response on a turbine subsystem during a specific time interval. Since the load distribution depends on the weather condition, POE can be expressed as follows:
$PT=P(Y>lT)$
(1)

$=∫xinxoutP(Y>lT|X=x)f(x)dx$
(2)

where lT 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, xin m/s, to the cut-out speed, xout 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].

Assessing the extreme load implies determining the resistance level, lT, that satisfies Eqs. (1) and (2), given the target POE level, PT [5]. Mathematically, the extreme load refers to the extreme quantile, lT, given by
$lT=inf{l∈ℝ:P(Y>l)≤PT}$
(3)

where “$inf$” implies the infimum.

Note that the extreme load level depends on the target POE, PT. When PT is small, a large extreme quantile, lT, needs to be estimated, and vice versa. When a wind turbine is designed to resist the resulting load level, lT, the anticipated load exceedance probability is smaller than, or equal to, the target level, PT; this is why lT 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.

## Methodology

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

### SIS Methods.

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, xi, $i=1,…,M$. Then, it runs the simulator Ni times at each sampled xi, $i=1,…,M$, where Ni 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, xi. With these outputs, the POE estimate of SIS1 is

$P̂SIS1(Y>l)=1M∑i=1MP̂ (Y>l|X=xi)f(xi)q(xi)$
(4)

$=1M∑i=1M(1Ni∑j=1NiI(Yj(i)>l))f(xi)q(xi)$
(5)

where $I(·)$ is an indicator function. The total number of simulation runs is $NT=∑i=1MNi$. Note that at each sampled xi, the Ni 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.

According to Ref. [17], the variance of SIS1 POE estimator can be minimized when the following optimal SIS density, $qSIS1$, and allocation size, Ni, $i=1,…,M$, are used:
$qSIS1 (x)=1Cq1f (x)1NTs (x)(1−s (x))+s (x)2$
(6)

$Ni=NTNT(1−s(xi))1+(NT−1)s(xi)∑j=1MNT(1−s(xj))1+(NT−1)s(xj), i=1,…,M$
(7)

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 Ni to those xi'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]).

The SIS2 method is similar to SIS1, except that it restricts the allocation size, Ni, to one. SIS2 samples wind speeds from a biased density, q, and runs the simulator once at each sampled wind speed. The resulting POE estimate of SIS2 is
$P̂SIS2(Y>l)=1NT∑i=1NTI(Yi>l)f(xi)q(xi)$
(8)

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

The variance of SIS2 POE estimator is minimized with the following optimal SIS density:
$qSIS2 (x)=1Cq2f (x)s (x)$
(9)

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

Recall that the extreme load corresponds to the extreme quantile, lT, as shown in Eq. (3). Given the target POE level, PT, the extreme quantile can be estimated by the following quantile estimator [21]:
$l̂PT=inf{l∈ℝ:P̂SIS (Y>l)≤PT}$
(10)

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

#### Extreme Load Assessment With SIS1.

We use the NT load outputs (i.e., NT 10 min maximum load responses) from NT simulation runs. First, we sort the NT 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.

With the largest load output, $Y(NT)$, the POE estimate is
$P̂SIS1 (Y>Y(NT))=1M∑i=1M(1Ni∑j=1NiI(Yj(i)>Y(NT)))f(xi)qSIS1(xi)=0$

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

With the second largest load output, $Y(NT−1)$, we obtain the POE estimate
$P̂SIS1 (Y>Y(N−1))=1M∑i=1M(1Ni∑j=1NiI(Yj(i)>Y(NT−1)))f(xi)qSIS1(xi)=f(x(NT))M·N(NT)·qSIS1(x(NT))$
(11)
where $N(NT)$ in Eq. (11) is the allocation size at $x(NT)$, where the largest load event, $Y(NT)$, occurs. Equation (11) holds because only the largest load output, $Y(NT)$, is greater than the second largest output, $Y(NT−1)$. Then, we plug the POE level in Eq. (11) into PT in Eq. (10) to obtain
$l̂PT=inf{l∈ℝ:P̂SIS (Y>l)≤f(x(NT))M·N(NT)·qSIS1(x(NT))}$
(12)

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 PT, which can be estimated using NT outputs.

In general, the POE estimate with the kth smallest (or equivalently, $(NT−k+1)th$ largest) load output, $Y(k), k=1,…,NT$, is
$P̂SIS1 (Y>Y(k))=1M∑i=1M(1Ni∑j=1NiI(Yj(i)>Y(k)))f(xi)qSIS1(xi)=1M∑i=1M∑j=1Ni(1NiI(Yj(i)>Y(k))f(xi)qSIS1(xi))=1M∑a=1NT(1N(a)I(Y(a)>Y(k))f(x(a))qSIS1(x(a)))$
(13)

$=1M∑a>k(1N(a)f(x(a))qSIS1(x(a)))$
(14)

In Eq. (13), we sort the NT 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)))$.

#### Extreme Load Assessment With SIS2.

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

With the largest load output, $Y(NT)$, the POE estimate becomes zero as in SIS1. With the second largest load output, $Y(NT−1)$, the POE estimate is given by
$P̂SIS2 (Y>Y(NT−1))=1NT∑i=1NTI(Yi>Y(NT−1))f(xi)qSIS2(xi)=f(x(NT))NT·qSIS2(x(NT))$
(15)

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

In general, the kth smallest (or equivalently, $(NT−k+1)th$ largest) load output, $Y(k), k=1,…,NT$, becomes the extreme load estimate associated with the POE level
$P̂SIS2 (Y>Y(k))=1NT∑i=1NTI(Yi>Y(k))f(xi)qSIS2(xi)=1NT∑i>kf(x(i))qSIS2(x(i))$
(16)

### Implications.

The standard MCS method provides an unbiased estimation of the POE as follows [22]:
$P̂MCS(Y>l)=1NT∑i=1NTI(Yi>l)$
where Yi is the load response from the ith simulation run at the wind speed independently sampled from the original density, f. After sorting the NT outputs, with the second largest output, $Y(NT−1)$, from MCS, the POE estimate becomes
$P̂MCS (Y>Y(NT−1))=1NT$

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))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.

## Implementation Details

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

### Wind Turbine Simulator.

We use the NREL wind turbine simulators, TurbSim [7] and FAST [8]. We use a 10 min simulation as one simulation run and take the maximum in-plane tip deflection or flapwise bending moment as a response variable, Y. Note that in Refs. [14] and [17], the flapwise bending moments are computed using the outputs from FAST. However, in our implementation, we use the flapwise bending moments directly generated from FAST as in Refs. [15] and [23]. For the original input density, f, we use the following truncated Rayleigh density:
$f(x)=fR(x)FR(xout)−FR(xin)$

where fR 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 xin = 3 m/s and xout = 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.

### Obtaining SIS Densities.

In $qSIS1$, Ni, $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 ql 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, qSIS, 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.

### Procedure Summary.

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), Ni, $i=1,…,M$, in Eq. (7) and $qSIS2$ in Eq. (9).

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

Step 3: Sort the NT 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 NT. We use M = 500 and $NT=3000$.

## Results

### Comparison With MCS.

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 PT 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, qSIS, 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.

### Estimating Extreme Loads Associated With Smaller POEs.

To illustrate how the proposed approach can estimate the extreme loads associated with smaller POEs (e.g., the POE of 3.8 × 10−7 corresponding to the 50-year return period [9]) without extrapolation, we combine the multiple repetitions to use all the simulation results. The POE estimate with combined repetitions is
$P̂SIS1 (Y>l)=1NR∑r=1NR1M∑i=1M(1Nri∑j=1NriI(Yj(ri)>l))f(xri)q(xri)$
(17)
and
$P̂SIS2 (Y>l)=1NR∑r=1NR1NT∑i=1NTI(Yri>l)f(xri)q(xri)$
(18)

for SIS1 and SIS2, respectively. Here, NR is the number of repetitions; xri and Nri are the ith sampled input and ith allocation size, respectively, in the rth repetition; and $Yj(ri)$ is the jth output at xri 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 NT = 3000 leads to the largest extreme load estimate of 2.67 m corresponding to the POE of $1/48,000=$ 2.1 × 10−5. 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−8 (250-year return period) and 3.5 × 10−7 (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 × 106 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−6 ($1/60,000$), 4.8 × 10−6, and 6.2 × 10−6, 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.

### Comparison With the Binning-Based Load Extrapolation Method.

As one of the most common load extrapolation methods in practice, we mentioned the binning method in Sec. 1. Mathematically, this method is based on the fact that POE in Eq. (2) can be expressed as
$P(Y>lT)=∫xinxoutP(Y>lT|X=x)f(x)dx=∑b=1BP(Y>lT|xb−1
(19)

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

We can obtain the probability, $P(xb−1, in Eq. (19), because the distribution of X, f(x), is known. However, because the conditional probability, $P(Y>lT∣xb−1, is unknown, the binning method approximates it, using the estimator
$P̂BIN=∑b=1BP̂(Y>lT|xb−1
(20)

where the approximation, $P̂(Y>lT∣xb−1, 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, 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, lT, 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.

## Summary

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.

## Acknowledgment

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).

## References

References
1.
Fossum
,
P. K.
,
Frøyd
,
L.
, and
Dahlhaug
,
O. G.
,
2013
, “
Design and Fatigue Performance of Large Utility-Scale Wind Turbine Blades
,”
ASME J. Sol. Energy Eng.
,
135
(
3
), p.
031019
.
2.
Byon
,
E.
,
2013
, “
Wind Turbine Operations and Maintenance: A Tractable Approximation of Dynamic Decision-Making
,”
IIE Trans.
,
45
(
11
), pp.
1188
1201
.
3.
Asmus
,
P.
,
2010
, “
The Wind Energy Operations and Maintenance
,” Wind Energy Update, London, UK.
4.
Byon
,
E.
,
Ntaimo
,
L.
, and
Ding
,
Y.
,
2010
, “
Optimal Maintenance Strategies for Wind Power Systems Under Stochastic Weather Conditions
,”
IEEE Trans. Reliab.
,
59
(
2
), pp.
393
404
.
5.
Lee
,
G.
,
Byon
,
E.
,
Ntaimo
,
L.
, and
Ding
,
Y.
,
2013
, “
Bayesian Spline Method for Assessing Extreme Loads on Wind Turbines
,”
Ann. Appl. Stat.
,
7
(
4
), pp.
2034
2061
.
6.
Agarwal
,
P.
, and
Manuel
,
L.
,
2008
, “
Extreme Loads for an Offshore Wind Turbine Using Statistical Extrapolation From Limited Field Data
,”
Wind Energy
,
11
(
6
), pp.
673
684
.
7.
Jonkman
,
B. J.
,
2009
, “
TurbSim User's Guide: Version 1.50
,” National Renewable Energy Laboratory, Golden, CO, Technical Report No.
NREL
/TP-500-46198.
8.
Jonkman
,
J. M.
, and
Buhl
,
M. L.
, Jr.
,
2005
, “
FAST User's Guide
,” National Renewable Energy Laboratory, Golden, CO, Technical Report No.
NREL
/EL-500-38230.
9.
IEC
,
2005
, “
Wind Turbines—Part 1: Design Requirements
,” International Electrotechnical Commission, Geneva, Switzerland, IEC/TC88, 61400–1 ed. 3.
10.
Toft
,
H. S.
,
Sørenson
,
J. D.
, and
Veldkamp
,
D.
,
2011
, “
Assessment of Load Extrapolation Methods for Wind Turbines
,”
ASME J. Sol. Energy Eng.
,
133
(
2
), p.
021001
.
11.
Regan
,
P.
, and
Manuel
,
L.
,
2008
, “
Statistical Extrapolation Methods for Estimating Wind Turbine Extreme Loads
,”
ASME J. Sol. Energy Eng.
,
130
(
3
), p.
031011
.
12.
Fogle
,
J.
,
Agarwal
,
P.
, and
Manuel
,
L.
,
2008
, “
Towards an Improved Understanding of Statistical Extrapolation for Wind Turbine Extreme Loads
,”
Wind Energy
,
11
(
6
), pp.
613
635
.
13.
Natarajan
,
A.
, and
Holley
,
W.
,
2008
, “
,”
ASME J. Sol. Energy Eng.
,
130
(
3
), p.
031017
.
14.
Moriarty
,
P.
,
2008
, “
Database for Validation of Design Load Extrapolation Techniques
,”
Wind Energy
,
11
(
6
), pp.
559
576
.
15.
Manuel
,
L.
,
Nguyen
,
H. H.
, and
Barone
,
M. F.
,
2013
, “
On the Use of a Large Database of Simulated Wind Turbine Loads to Aid in Assessing Design Standard Provisions
,”
AIAA
Paper No. 2013-0197.
16.
Cannamela
,
C.
,
Garnier
,
J.
, and
Iooss
,
B.
,
2008
, “
Controlled Stratification for Quantile Estimation
,”
Ann. Appl. Stat.
,
2
(
4
), pp.
1554
1580
.
17.
Choe
,
Y.
,
Byon
,
E.
, and
Chen
,
N.
,
2015
, “
Importance Sampling for Reliability Evaluation With Stochastic Simulation Models
,”
Technometrics
,
57
(
3
), pp.
351
361
.
18.
Yampikulsakul
,
N.
,
Byon
,
E.
,
Huang
,
S.
,
Sheng
,
S.
, and
You
,
M.
,
2014
, “
Condition Monitoring of Wind Power System With Nonparametric Regression Analysis
,”
IEEE Trans. Energy Convers.
,
29
(
2
), pp.
288
299
.
19.
Lee
,
G.
,
Ding
,
Y.
,
Xie
,
L.
, and
Genton
,
M. G.
,
2015
, “
A Kernel Plus Method for Quantifying Wind Turbine Performance Upgrades
,”
Wind Energy
,
17
(
7
), pp.
1207
1219
.
20.
Saranyasoontorn
,
K.
,
Manuel
,
L.
, and
Veers
,
P. S.
,
2004
, “
A Comparison of Standard Coherence Models for Inflow Turbulence With Estimates From Field Measurements
,”
ASME J. Sol. Energy Eng.
,
126
(
4
), pp.
1069
1082
.
21.
Glynn
,
P. W.
,
1996
, “
Importance Sampling for Monte Carlo Estimation of Quantiles
,”
Second International Workshop on Mathematical Methods in Stochastic Simulation and Experimental Design
, St. Petersburg, Russia, June 18–22, pp.
180
185
.
22.
Kroese
,
D. P.
,
Taimre
,
T.
, and
Botev
,
Z. I.
,
2011
,
Handbook of Monte Carlo Methods
,
Wiley
,
New York
.
23.
Barone
,
M.
,
Paquette
,
J.
,
Resor
,
B.
, and
Manuel
,
L.
,
2012
, “
,”
AIAA
Paper No. 2012–1288.
24.
Jonkman
,
J. M.
,
Butterfield
,
S.
,
Musial
,
W.
, and
Scott
,
G.
,
2009
, “
Definition of a 5-MW Reference Wind Turbine for Offshore System Development
,” National Renewable Energy Laboratory, Golden, CO, Technical Report No.
NREL
/TP-500-38060.
25.
Agarwal
,
P.
, and
Manuel
,
L.
,
2008
, “
The Influence of the Joint Wind-Wave Environment on Offshore Wind Turbine Support Structure Loads
,”
ASME J. Sol. Energy Eng.
,
130
(
3
), p.
031010
.