Maximum Instantaneous Wind Speed Forecasting and Performance Evaluation by Using Numerical Weather Prediction and On-Site Measurement

Abstract

A maximum instantaneous wind speed forecast methodology based on the autoregressive with exogenous inputs (ARX) model is proposed, in which numerical weather prediction and on-site measurement are used as inputs and the model parameters are estimated using non-parametric regression with forgetting factors. The accuracy of prediction using a proposed dynamic model is then evaluated and compared to the conventional static model output statistics (MOS) model. It is found that the prediction accuracy is improved by utilizing numerical weather prediction with a higher horizontal resolution. Finally, the predictability of the maximum instantaneous wind speed higher than 15m/s is evaluated using the receiver operating characteristic (ROC) curve and the area under the curve (AUC). The optimal quantile level of the maximum instantaneous wind speed is derived using a cost function.

1. Introduction

Many forms of infrastructure such as railways and roads are constructed in complex terrain in Japan. The operation and maintenance of these infrastructures are influenced by strong gusts of wind. For example, train operation has to be suspended when a certain value of instantaneous wind speed is recorded. This requires the forecasting of strong gusts of wind, including information on when the strong gust of wind happens and the intensity of the gust, with a forecast horizon of up to 24 h. Since decisions concerning operation are usually based on the instantaneous wind speed, forecasting of instantaneous wind speed is needed, as shown in Misu and Ishihara [1].

Kobayashi and Shimamura [2] proposed a method of predicting the instantaneous wind speed up to 15 min ahead based on the on-site measurement data and the Kalman filter for the application of train operation regulation. Hoppmann et al. [3] proposed a method for the prediction of instantaneous wind speed with a confidence interval up to two minutes based on past measurement data using a linear trend model for railway operations in Germany. Nevertheless, these methods based on previously measured data do not have the ability to perform forecasts up to 24 h ahead, and for this it is necessary to use forecast methods which can consider weather conditions, such as numerical weather prediction data.

Accordingly, much research has been conducted on wind power forecasting of up to 24 h, which is needed for power grid operation and power trading [4,5]. The most important part of wind power forecasting is the forecasting of the wind speed, which can be done using a combination of a physical model, such as numerical weather prediction (NWP) and a statistical model based on on-site measurement. Landberg [6,7] predicted the wind speed in a wind farm based on numerical weather prediction and measured wind speed using a statistical method called model output statistics (MOS). In MOS, the accuracy of the forecast is improved by using a linear regression model based on the available measurement data for each wind direction. Although this method reduces the bias of the forecast, long-term measurement data are required and the method cannot be used for cases in which the regression model changes depending on the season or in which the linear regression model is not applicable. Advances in wind power forecasting have been summarized by Sideratos and Hatziargyriou [8], Soman et al. [9], Giebel and Kariniotakis [10], Dittner and Vasel [11] and Dhiman et al. [12].

The autoregressive with exogenous inputs (ARX) model was proposed by Joensen et al. [13], in which the real-time measurement of the power generation output and numerical weather prediction are used as the input, and parametric regression with forgetting factors is used for the estimation of the model parameters and functions. They performed the forecast of the 1-h average wind power output up to 24 h ahead and verified the results at a wind farm in Denmark. The properties required for non-parametric probabilistic forecasts of wind power were pointed out by Pinson et al., along with an evaluation of their use [14]. Although the average power output can be forecasted with high accuracy by this model, no model has been proposed for the forecasting of the maximum instantaneous wind speed. Moreover, its applicability in complex terrain has not been clarified.

In this study, a forecast model for the maximum instantaneous wind speed using the ARX model is proposed in Section 2. Its applicability to maximum instantaneous wind speed forecasting in a mountainous area in Japan is then clarified and the effects of horizontal resolution on the numerical weather forecast data, used as inputs on the forecast accuracy, are investigated in Section 3. Furthermore, the predictability of strong gusts of wind with wind speeds higher than 15 m/s is quantitatively evaluated using receiver operating characteristic (ROC) curves and the area under the curve (AUC). The optimal quantile level of the maximum instantaneous wind speed is derived using a cost function.

2. Forecast Model

In this study, a model to forecast the maximum instantaneous wind speed up to 24 h ahead is developed based on the numerical weather prediction and on-site measurement of wind speed. Figure 1 shows an overview of the model. The first step of this model consists of a mean wind speed forecast model, a fluctuating wind speed forecast model and a peak factor forecast model. Each model performs a forecast based on numerical weather prediction and on-site measurement. Then, a maximum instantaneous wind speed forecast model calculates the expected valued of the maximum instantaneous wind speed by adding the fluctuating wind speed multiplied by the peak factor to the mean wind speed. Finally, the maximum instantaneous wind speed with a quantile level is evaluated by considering the forecast error of the expected value.

The numerical weather prediction and the on-site measurement, which are used as the input to the model, are first presented in Section 2.1. The mean wind speed forecast model, fluctuating wind speed forecast model, peak factor forecast model and the maximum instantaneous wind speed forecast model are then described in Section 2.2. Finally, non-parametric regression with forgetting factors is presented in Section 2.3.

2.1. Input Data

The Japan Meteorological Agency provides several numerical weather prediction datasets with different resolutions and forecast horizons, as shown in

Table 1. Summary of numerical weather prediction data used in this study.

Model	GPV-MSM	GPV-GSM (Japan Area)
Forecast horizon	39 h	84 h
Delivery time	About 3 h after initial time	About 3 h after initial time
Temporal resolution	1 h (surface), 3 h (barometric level)
Initial time (JST)	3:00, 9:00, 15:00, 21:00	3:00, 9:00, 15:00, 21:00
Forecast variables	Sea surface rehabilitation pressure (ground surface), altitude (pressure surface), horizontal wind, updraft, temperature, relative humidity, accumulated precipitation
Number of vertical layers	60 layers
Horizontal discretization scheme and resolution	Grid point model resolution: 20 km	Spectrum model cut-off wave number: 519
Output horizontal grid resolution (surface)	0.05 degrees north-south × 0.0625 degrees east-west	0.2 degrees north-south × 0.25 degrees east-west
Output Domain	22.4° N–47.6° N, 120° E–150° E	20° N–50° N, 120° E–150° E

. In this study, the GPV-MSM dataset, which has the finest horizontal resolution of about 5 km at the surface, and the GPV-GSM dataset (Japan area), which has the second-finest horizontal resolution of about 20 km at the surface, are used in order to investigate the effect of the resolution of the NWP. The forecast interval of the GPV-MSM is three hours (i.e., eight forecasts per day), though only four forecasts per day are used for the sake of the comparison with GPV-GSM (Japan area). Among the forecast variables available in these NWPs, only the east–west component and the north–south components of the surface wind velocity are used as input data for the forecast model in this study.

In the proposed model, real-time on-site measurement of average wind speed, average wind direction, standard deviation of the wind speed and maximum instantaneous wind speed are used as input data. As described in Figure 1, the average wind speed and wind direction are the inputs to the average wind speed forecast model, the standard deviation of the wind speed is the input to the fluctuating wind speed forecast model, and the average wind speed and the standard deviation of the wind speed, as well as the maximum instantaneous wind speed, are the inputs to the peak factor forecast model. The averaging time of the wind speed can be set between 10 min and 1 h. In this study the averaging time is set to 30 min because the maximum instantaneous wind speed is forecasted every 30 min based on the 30-min average wind speed and the standard deviation of the wind speed for every 30 min. The forecast is also conducted every 30 min.

Figure 2 shows the forecast schedule when GPV-MSM is used as the input NWP. Forecast data is delivered online from the Meteorological Service Support Center approximately three hours after the initial time but may be delayed for up to three hours. For this reason, each forecast dataset is considered to become available six hours after the initial time and can be used for the wind speed forecast. For example, if the question of whether a strong wind event will occur during the daytime has to be determined at 6:00 am, the NWP data that can be used is based on the initial time at 21:00. This data is the forecast value from the initial value (21:00 the previous day) until 12:00 the next day, which is 39 h after the initial time.

2.2. Forecast Model

The NWP does not contain the effects of high-resolution terrain and surface roughness because of the limitation of the spatial resolution. In this study, the local forecast is used in which the local wind speed, affected by the high-resolution topography and surface roughness, is expressed as a function of wind speed and direction of NWP and is identified by non-parametric regression using the measurements and the NPW data. After $\hat{k}$ steps ($k$ hours) at time $t$, the local wind speed forecast value | $\left|{\overline{u}}_{t+k|t}^{\mathrm{local}}\right|$ is modeled as a function of the NWP wind speed $\left|{\overline{u}}_{t+k|t}^{\mathrm{nwp}}\right|$ and NWP wind direction ${\mathsf{\theta }}_{t+k|t}^{\mathrm{nwp}}$ as shown in Equation (1), where the nearest neighbor of the output grid is used to extract exact-location wind forecasts from the NWP and is expressed as

$$\left|{\overline{u}}_{t+k|t}^{\mathrm{local}}\right|=f\left(\left|{\overline{u}}_{t+k|t}^{\mathrm{nwp}}\right|,{\theta }_{t+k|t}^{\mathrm{nwp}}\right)$$

(1)

Note that $k=\hat{k}\Delta t$. $\left|{\overline{u}}_{t+k|t}^{\mathrm{nwp}}\right|$ and ${\theta }_{t+k|t}^{\mathrm{nwp}}$ are the absolute value of the mean wind speed and wind direction of NWP, which can be calculated from the east–west component $u_{t+k|t}^{\mathrm{nwp}}$ and the northsouth component $v_{t+k|t}^{\mathrm{nwp}}$ of the wind speed, as shown in Equation (2).

$$\left|u_{t+k|t}^{\mathrm{nwp}}\right|=\sqrt{{\left(v_{t+k|t}^{\mathrm{nwp}}\right)}^2+{\left(u_{t+k|t}^{\mathrm{nwp}}\right)}^2},{\theta }_{t+k|t}^{\mathrm{nwp}}={\tan }^{-1}\left(v_{t+k|t}^{\mathrm{nwp}},u_{t+k|t}^{\mathrm{nwp}}\right)$$

(2)

The function $f\left(\left|{\overline{u}}_{t+k|t}^{\mathrm{nwp}}\right|,{\theta }_{t+k|t}^{\mathrm{nwp}}\right)$ is estimated using non-parametric regression with a forgetting factor, which is set to the value of 0.999, as proposed by Joensen et al. [13].

The local wind speed obtained using Equation (1) is considered to contain no bias component of error, but the phase error, which is originally included in NWP, is still included. To reduce this error for short-term forecasts, a model in which the weighted sum of the persistent forecast, which is a forecast that the future weather condition will be the same as the present condition, and the local forecast obtained by Equation (1) are used as the final forecast value, as

$$\left|{\overline{u}}_{t+k|t}^{\mathrm{pred}}\right|=a\left(k,{\theta }_{t+k|t}^{\mathrm{nwp}}\right)\left|{\overline{u}}_t^{\mathrm{meas}}\right|+b\left(k,{\theta }_{t+k|t}^{\mathrm{nwp}}\right)\left|{\overline{u}}_{t+k|t}^{\mathrm{local}}\right|,$$

(3)

where $a$ and $b$ are the smooth functions of the forecast horizon $k$ and the wind direction of NWP, and are estimated by non-parametric regression with the forgetting factor of 0.999 [13].

The fluctuating component of the wind speed ${\sigma }_{t+k|t}^{\mathrm{local}}$ is also assumed to be a function of the mean wind speed and wind direction of NWP.

$${\sigma }_{t+k|t}^{\mathrm{local}}=f_{\sigma }\left(\left|u_{t+k|t}^{\mathrm{nwp}}\right|,{\theta }_{t+k|t}^{\mathrm{nwp}}\right)$$

(4)

Similar to the average wind speed, the forecast value of the fluctuating wind speed ${\sigma }_{t+k|t}^{\mathrm{pred}}$ is modeled as the weighted sum of the persistent forecast and the local forecast so that the phase error of the forecast can be reduced by Equation (5).

$${\sigma }_{t+k|t}^{\mathrm{pred}}=a_{\sigma }\left(k,{\theta }_{t+k|t}^{\mathrm{nwp}}\right){\sigma }_t^{\mathrm{meas}}+b_{\sigma }\left(k,{\theta }_{t+k|t}^{\mathrm{nwp}}\right){\sigma }_{t+k|t}^{\mathrm{local}},$$

(5)

where $a_{\sigma }$ and $b_{\sigma }$ are the smooth functions of the forecast horizon $k$ and the wind direction of the NWP, and are estimated by non-parametric regression with the forgetting factor of 0.999 [13].

Based on the average wind speed $\left|{\overline{u}}_{t+k|t}^{\mathrm{pred}}\right|$ and the fluctuating component of the wind speed ${\sigma }_{t+k|t}^{\mathrm{pred}}$, the expected value of the maximum instantaneous wind speed $u_{t+k|t}^{\max ,\mathrm{pred}}$ is obtained using Equation (6).

$$u_{t+k|t}^{\max ,\mathrm{pred}}=\left|{\overline{u}}_{t+k|t}^{\mathrm{pred}}\right|+p_t{\sigma }_{t+k|t}^{\mathrm{pred}}$$

(6)

where $p_t$ is a peak factor and is assumed to be a value that depends only on the initial forecast time, not on the wind speed, wind direction or forecast period. It can be estimated as a constant value by using the regression with the forgetting factor, in which the error shown in Equation (7) is minimized. Since the peak factor fluctuates greatly depending on the weather conditions, the forgetting factor corresponds to a relatively short time scale, and 0.917 is used.

$$ϵ=p_t-\frac{u_t^{\max ,\mathrm{meas}}-\left|{\overline{u}}_t^{\mathrm{meas}}\right|}{{\sigma }_t^{\mathrm{meas}}}$$

(7)

where $u_t^{\max ,\mathrm{meas}}$, ${\sigma }_t^{\mathrm{meas}}$ and $\left|{\overline{u}}_t^{\mathrm{meas}}\right|$ are the maximum instantaneous, standard deviation and mean of the measured wind speed.

Since the proposed model is based on the average relationship between the predicted and measured wind speeds, the final output will be the expected value of the maximum instantaneous wind speed. During an actual strong wind event, the maximum instantaneous wind speed may fluctuate around this value. In this study, the standard deviation of this fluctuation can be modeled by using the mean square error of the maximum instantaneous wind speed forecast of the past ${ϵ}_t\left(k\right)$ and only the function of the forecast horizon $k$, as shown in Equation (8). When the maximum instantaneous wind speed forecast is carried out, Equation (9) is used to estimate the quantile level of the maximum instantaneous wind speed, which is important for engineering applications. The value of $\gamma$ should be determined depending on the application and will be discussed in Section 3.

$${ϵ}_t\left(k\right)=\sqrt{\frac{1}{N}{\displaystyle \sum }_t^{}{\left(u_{t+k}^{\max ,\mathrm{meas}}-u_{t|t+k}^{\max ,\mathrm{pred}}\right)}^2}$$

(8)

$$u_{t+k|t}^{\max }=u_{t+k|t}^{\max ,\mathrm{pred}}+\gamma {ϵ}_t\left(k\right)$$

(9)

2.3. Non-Parametric Regression with Forgetting Factor

In this study, non-parametric regression with a forgetting factor is used to estimate the parameters or functions of each forecast model described in Section 2.2. This method can be formulated in the form of Equation (10).

$$y_s=z_s{\mathsf{\varphi }}^T\left(q_s\right)+{ϵ}_s,$$

(10)

where $z_s={\left(\begin{array}{cccc}{z}_{s\left(1\right)}& z_{s\left(2\right)}& \cdots & z_{s\left(M\right)}\end{array}\right)}^T$ and $q_s=\left(\begin{array}{cccc}{q}_{s\left(1\right)}& q_{s\left(2\right)}& \cdots & q_{s\left(N\right)}\end{array}\right)$ are the explanatory variables, $y_s$ is the objective variable, ${\mathsf{\varphi }}^T\left(q_s\right)=\left(\begin{array}{cccc}{\varphi }_{\left(1\right)}\left(q\right)& {\varphi }_{\left(2\right)}\left(q\right)& \cdots & {\varphi }_{\left(M\right)}\left(q\right)\end{array}\right)$ is a smooth function of $q_s$ estimated by non-parametric regression and ${ϵ}_s$ is an model error modeled as white noise. The subscript s is the index of time series data. The models proposed in Section 2.2 can be expressed in the form of Equation (10). For example, Equation (1) can be shown as Equation (10) with the parameters defined in Equation (11):

$$\begin{gathered} M=1, \\ N=2, \\ {\varphi }_{\left(1\right)}\left(q\right)=f\left(q_{\left(1\right)},q_{\left(2\right)}\right), \\ \{\begin{array}{c}{q}_{\left(1\right)}=u_{t+k|t}\\ q_{\left(2\right)}={\theta }_{t+k|t}\end{array}, \\ {z}_{\left(1\right)}=1, \end{gathered}$$

(11)

Equation (3) can be shown as Equation (10) with the parameters defined in Equation (12) as

$$\begin{gathered} M=2, \\ N=2, \\ \{\begin{array}{c}{\varphi }_{\left(1\right)}\left(q\right)=a\left(q_{\left(1\right)},q_{\left(2\right)}\right)\\ {\varphi }_{\left(2\right)}\left(q\right)=b\left(q_{\left(1\right)},q_{\left(2\right)}\right)\end{array} \\ \{\begin{array}{c}{q}_{\left(1\right)}=k \\ q_{\left(2\right)}={\theta }_{t+k|t}\end{array} \\ \{\begin{array}{c}{z}_{\left(1\right)}=u_{t+k|t}^{\mathrm{local}}\\ z_{\left(2\right)}=u_{t+k|t}^{\mathrm{meas}}\end{array} \end{gathered}$$

(12)

In non-parametric regression, a function $\mathsf{\varphi }\left(q_s\right)$ is estimated by using the pair of explanatory variables and training data obtained in the past locally around the described points $q_p={\left(\begin{array}{cccc}{q}_{p\left(1\right)}& q_{p\left(2\right)}& \cdots & q_{p\left(N\right)}\end{array}\right)}^{T }$. When the training data obtained at the time s is defined as $y_s$, $z_s$ and $q_s={\left(\begin{array}{cccc}{q}_{s\left(1\right)}& q_{s\left(2\right)}& \cdots & q_{s\left(N\right)}\end{array}\right)}^T$, the function ${\hat{\mathsf{\varphi }}}_{p,t}\left(q\right)$ near $q_p$ at time $t$ is obtained by minimizing the weighted error, as shown in Equation (13).

$$ϵ={\displaystyle \sum }_{s=1}^t{\lambda }^{t-s}w\left(q_s,q_p\right){\left(y_s-z_s^T{\hat{\mathsf{\varphi }}}_{p,t}\left(q\right)\right)}^2$$

(13)

where the weight ${\lambda }^{t-s}w\left(q_s,q_p\right)$ can be divided into two parts as follows.

$\lambda$ is a model parameter called a forgetting factor that takes a value in the range of $0<\lambda \le 1$. Since $t$ is the current time and $s$ is the time when the measurement data were obtained, ${\lambda }^{t-s}$ is small for the old measurement data and large for the new data. By using this forgetting factor, a dynamic model which forgets past measurement data and emphasizes the weight of new measurement data is realized. When a larger value of $\lambda$ is used, longer period data is remembered, and when a smaller value of $\lambda$ is used, the training focuses more on the latest data. An index $N_{eff}$ has been proposed [12], as shown in Equation (14), to show the effective number of data used for the training when different $\lambda$ values are used.

$$N_{eff}=\frac{\lambda }{1-\lambda }$$

(14)

For example, if $\lambda =1$, the effective number of data used will be infinite, which is equivalent to using all the past measurement data. In the field of wind power forecasting, 0.999 is used as a forgetting factor [12].

$w\left(q_s,q_p\right)$ is a weight function that takes a large value when the distance is small between the training data $q_s$ at time $s$ and the point $q_p$, in the vicinity of which the function is estimated as small, and takes a small value when the distance is large. With this weight function, it is possible to increase the weight of the training data in the vicinity of the estimation point of the function. The multidimensional weight function $w\left(q_s,q_p\right)$ is expressed as the product of the one-dimensional weight function $W\left(x\right),$ as shown in Equation (15).

$$w\left(q_s,q_p\right)={\displaystyle \prod }_{j=1}^{N}W\left(\frac{\left|q_{s\left(j\right)}-q_{p\left(j\right)}\right|}{{\hslash }_j}\right)$$

(15)

where ${\hslash }_j$ is a quantity called bandwidth, which is a parameter that determines the smoothness of the estimated function. A function shown in Equation (16) is used as the one-dimensional weight function $W\left(x\right)$.

$$W\left(x\right)=\{\begin{array}{cc}{\left(1-x^3\right)}^3& 0\le x<1\\ 0& 1\le x\end{array}$$

(16)

$z_s$ and $\lambda$, which minimize the error shown in Equation (13), can be calculated using the recurrence formula. The details are shown in Appendix A.

3. Forecast Results and Error Evaluation

First, the strong wind events and analysis conditions used in this study are explained in Section 3.1. Next, examples of forecasts performed in the study and models for comparison are introduced and the effect of differences in numerical weather forecast data used as inputs on forecast accuracy is clarified in Section 3.2. The performance of the forecast of the strong wind event with a maximum gust wind speed of 15 m/s or more is then quantitatively evaluated using the ROC curve and AUC in Section 3.3. Finally, the optimal quantile level of the maximum instantaneous wind speed is discussed in Section 3.4.

3.1. Analysis Conditions

In this study, we used the wind speed measured by an anemometer at a height of 10 m above the ground level in the mountainous region in East Japan during three winters, i.e., from 6 January, 2014 to 31 March, 2014; from 6 December, 2014 to 31 March, 2015; and from 6 December, 2015 to 31 March, 2016. The anemometer used was a N-262LVS from Nippon Electric Instruments, set at 138.769167 degrees East and 36.944089 degrees North. The sampling frequency was 4 Hz and a 3-s averaged value was used as the instantaneous wind speed.

Table 2. List of targeted strong wind events and corresponding weather conditions

Date and Time Weather Conditions	Maximum Instantaneous Wind Speed (m/s)	Date and Time Weather Conditions	Maximum Instantaneous Wind Speed (m/s)
2014/01/30 16:00 South-coast cyclone	15.3	2015/01/31 07:30 Siberian High and Aleutian Low	15.5
2014/01/31 13:00 Siberian High and Aleutian Low	15.6	2015/02/01 10:00 Siberian High and Aleutian Low	15.0
2014/02/05 10:30 Siberian High and Aleutian Low	15.1	2015/02/13 15:00 Siberian High and Aleutian Low	15.9
2014/02/16 12:00 Siberian High and Aleutian Low	21.2	2015/02/15 14:30 Siberian High and Aleutian Low	16.0
2014/03/06 10:30 Siberian High and Aleutian Low	15.3	2015/02/26 07:30 South-coast cyclone	15.7
2014/03/10 12:30 Siberian High and Aleutian Low	16.2	2015/02/27 12:30 Siberian High and Aleutian Low	15.7
2014/03/20 06:30 South-coast cyclone	15.7	2015/03/01 06:00 South-coast cyclone	15.2
2014/03/21 06:30 Siberian High and Aleutian Low	16.2	2015/03/02 07:00 Siberian High and Aleutian Low	20.1
2014/03/30 10:00 South-coast cyclone	18.8	2015/03/03 16:30 South-coast cyclone	15.5
2014/03/31 06:00 Siberian High and Aleutian Low	19.4	2015/03/09 06:00 South-coast cyclone	15.4
2014/12/16 06:30 South-coast cyclone	17.8	2015/03/10 17:00 Siberian High and Aleutian Low	15.5
2014/12/18 15:30 Siberian High and Aleutian Low	17.6	2015/12/11 14:30 Siberian High and Aleutian Low	17.4
2014/12/20 07:00 South-coast cyclone	16.3	2016/01/04 18:00 Siberian High and Aleutian Low	15.9
2015/01/06 15:30 Siberian High and Aleutian Low	17.8	2016/01/18 06:00 South-coast cyclone	21.8
2015/01/07 07:00 Siberian High and Aleutian Low	16.2	2016/01/20 07:30 Siberian High and Aleutian Low	17.9
2015/01/17 13:00 Siberian High and Aleutian Low	15.1	2016/02/09 16:00 Siberian High and Aleutian Low	17.1
2015/01/22 10:30 South-coast cyclone	15.8	2016/02/10 06:00 Siberian High and Aleutian Low	16.7
2015/01/23 11:30 Siberian High and Aleutian Low	15.1	2016/03/01 06:00 Siberian High and Aleutian Low	16.8

shows a list of the 36 strong wind events with a maximum instantaneous wind speed of 15 m/s or more during the measurement campaign. All the 36 events are targeted for the maximum instantaneous wind speed forecast. The weather conditions for 36 strong wind events and corresponding weather conditions are also shown in Table 2. Figure 3 shows an aerial photograph and elevation contours of the area. The center of Figure 3 is the location of the anemometer. The measurement data required for correction of the local average wind speed and the local fluctuation wind speed were updated every 30 min using the latest data obtained at the time of the forecast. The initial values of the forecast model were trained by using the measurement data from December 2013 prior to the forecast evaluation.

The non-parametric regressions for each model were carried out under the conditions listed in

Table 3. Summary of model parameters.

Parameter		Value
Forgetting factor		0.917 (peak factor estimation) 0.999 (other than that)
Bandwidth	Wind speed	4.0 (m/s)
	Wind direction	11.25 (deg)
	Forecast horizon	0.5 (h)

. The forgetting factor $\lambda$ was set to 0.999 in the models other than the peak factor forecast model. The effective number of data in Equation (14) corresponds to approximately 20 days, i.e., the models are always trained by using the past 20 days of data, approximately. On the other hand, 0.917 was used for estimating the peak factor, which corresponds to about 5.5 h. When performing non-parametric regression, it is necessary to set the initial values of the local wind speed conversion and the wind speed correction functions. In this study, a forecast was performed in advance from 1 December to 31 December 2013, prior to the analysis period, and the initial values of each model were obtained by training.

3.2. Forecast Examples and Comparative Models

In this section, model functions are identified by non-parametric regression and an example of the identified functions is discussed. Figure 4 shows the function $f\left(\left|{\overline{u}}_{t+k|t}^{\mathrm{nwp}}\right|,{\theta }_{t+k|t}^{\mathrm{nwp}}\right)$ identified by the mean wind speed prediction model for the case of ${\theta }_{t+k|t}^{\mathrm{nwp}}=0°$ and ${\theta }_{t+k|t}^{\mathrm{nwp}}=180°$. It is found that the local wind speed generally does not vary linearly with the wind speed forecasted by the numerical weather prediction and also depends on the wind direction.

In the mean wind speed forecast model, the functions $a\left(k,{\theta }_{t+k|t}^{\mathrm{nwp}}\right)$ and $b\left(k,{\theta }_{t+k|t}^{\mathrm{nwp}}\right),$ which are used to take into account the on-site measurement, are also trained. Figure 5 shows the training results of these functions for ${\theta }_{t+k|t}^{\mathrm{nwp}}=180°$. When the forecast horizon $k$ is short, the weight of the latest measurement data is large, and if the forecast horizon $k$ time becomes longer, the weight of the forecast value based on the numerical weather forecast increases.

The concept of using the wind speed ratio for each wind direction has been widely used to evaluate the influence of local topography on wind conditions [1]. In addition, in the wind power output forecast model MOS, the wind speed is multiplied by a constant for each wind direction. Based on this idea, it is possible to perform the forecast by multiplying the coefficient, which is a function of each wind direction, to the numerical weather forecast data released by the Japan Meteorological Agency. In this study, MOS is a function of the NWP with GPV-GSM, and the maximum instantaneous wind speed forecast model based on this concept is constructed as shown in Equations (17) to (19) and is used as the forecast model for comparison. Average wind speed $\left|{\overline{u}}_{t+k|t}^{\mathrm{MOS}}\right|$ is expressed as

$$\left|{\overline{u}}_{t+k|t}^{\mathrm{MOS}}\right|=C\left({\theta }_{t+k|t}^{\mathrm{nwp}}\right)\left|{\overline{u}}_{t+k|t}^{\mathrm{nwp}}\right|$$

(17)

where $C\left({\theta }_{t+k|t}^{\mathrm{nwp}}\right)$ is a function of the wind direction. $S\left({\theta }_{t+k|t}^{\mathrm{nwp}}\right)$ in Equation (18), which is used to forecast the fluctuating wind speed ${\sigma }_{t+k|t}^{\mathrm{MOS}}$, is also a function of the wind direction and is identified basd on the past measurement data.

$${\sigma }_{t+k|t}^{\mathrm{MOS}}=S\left({\theta }_{t+k|t}^{\mathrm{nwp}}\right)\left|{\overline{u}}_{t+k|t}^{\mathrm{nwp}}\right|$$

(18)

The peak factor $p^{\mathrm{MOS}}$ is a constant and was determined based on the research by Ishizaki [15]. The maximum instantaneous wind speed $u_{t+k|t}^{\max ,\mathrm{MOS}}$ based on the static MOS model can be calculated by

$$u_{t+k|t}^{\max ,\mathrm{MOS}}=\left|{\overline{u}}_{t+k|t}^{\mathrm{MOS}}\right|+p^{\mathrm{MOS}}{\sigma }_{t+k|t}^{\mathrm{MOS}}$$

(19)

Figure 6 shows a comparison between the forecasted maximum instantaneous wind speed at 6:00 a.m. and measurement on the days when a typical strong wind event occurred. It is found that the prediction accuracy is improved by using the dynamic model proposed in this study compared to the conventional static MOS model, in which there is a large difference between the forecast value based on the numerical forecast and the observed value at 6:00 in the morning. There is also a significant difference between the MOS model and the dynamic model immediately after the start of the forecast. In addition, the prediction accuracy is better when the GPV-MSM data is used as the input NWP.

It is noted that the forecasted maximum instantaneous wind speeds in Figure 6 express the expected values at $\gamma =0$ as shown in Equation (9) and are lower than the optimal values. In this study, $\gamma =0.9$ was used as an optimal quantile level. This implies that the optimal forecasted maximum instantaneous wind speeds can match the measured ones, as shown in Figure 6. However, all model significantly underestimated the maximum instantaneous wind speeds on 30 March, 2014. This strong wind event was caused by atmospheric stratification and the mountain, as pointed out Saito and Ikawa [16]. The operational forecast models such as JMA-GSM or JMA-MSM utilize data assimilation of the global forecast data at all levels to achieve a high forecast score. However, this may be not the best solution for all cases. The use of an ensemble of models may improve the overall performance, as mentioned by Meng and Zhang [17], Mohrlen and Jorgensen [18].

3.3. Evaluation of Forecast Results

In this study, the accuracy of the forecast is evaluated using the root mean square error (RMSE) of the maximum instantaneous wind speed as a function of forecast horizon and the ROC curve of the strong wind event. The root mean square error of the maximum instantaneous wind speed is calculated for the entire forecast period of 11 months as a function of the forecast horizon $k$. Figure 7 shows the RMSE of maximum instantaneous wind speed predicted by the conventional static MOS model and the proposed dynamic method as functions of the forecast horizon. The static MOS model is almost constant regardless of the forecast horizon, whereas the proposed dynamic model uses the most up-to-date measurement data, and the forecast error within 6 h is greatly reduced. In addition, the use of the GPV-MSM forecast value as the input improves the forecast error for all forecast horizons compared to the GPV-GSM forecast value. This is due to the high horizontal resolution of the GPV-MSM model, which reflects the effects of mountainous terrain on the prediction.

The model for estimating the quantile level of the maximum instantaneous wind speed described in Section 2.2 is important for the capability of forecasting strong wind events, as shown in Table 2. Since the proposed dynamic model predicts the average value of the maximum instantaneous wind speed, it is necessary to estimate the quantile level of the maximum instantaneous wind speed in order to predict whether a strong wind event will occur. When the larger value of the parameter $\gamma$ is used, the high quantile level of the maximum instantaneous wind speed is predicted, and the number of predicted strong wind events increases. On the other hand, when the smaller value of the parameter $\gamma$ is used, the low quantile level of the maximum instantaneous wind speed decreases, and the number of predicted strong wind events is reduced.

It is common to use ROC curves [19,20] to evaluate forecasts including such model parameters. The ROC curve shows the false positive rate on the horizontal axis and the true positive rate (also called sensitivity) on the vertical axis and shows the change in the false positive rate and the true positive rate when the model parameters are changed.

Table 4. Event occurrence table.

		Measurement
		Yes	None
Forecast	Yes	a (True positive)	b (False positive)
	None	c (False negative)	d (True negative)

shows the definitions of the true positive $a$, false positive $b$, false negative $c$ and true negative $d$. The false positive rate is defined as $b/\left(b+d\right)$, which means $b$ cases were predicted with a strong wind event, in a total of $\left(b+d\right)$ cases in which no strong wind event occurred. On the other hand, the true positive rate can be defined as $a/(a+c)$, which means $a$ cases of strong wind events were predicted among $\left(a+c\right)$ cases in which the strong wind event occurred. With a perfect forecast, the true positive rate is 1 and the false positive rate is 0, which means the ROC curve passes through the upper left corner of the figure and the value of the AUC is equal to one. For the random model, the true positive rate and the false positive rate are expected to be equal and the ROC curve is a straight line which connects the lower left corner and the upper right corner. The ROC curve of the actual forecast model is located between these two curves. With a better forecast model, the ROC curve approaches the upper left corner. The forecast model can be evaluated quantitatively using the AUC. A larger AUC value indicates a better model [21].

In this study, the definition of the occurrence and forecast of a strong wind event is as follows. If a maximum instantaneous wind speed of 15 m/s or more is recorded between 6:00 a.m. and 18:00 p.m., it is defined as an occurrence of a strong wind event. Similarly, if a maximum instantaneous wind speed of 15 m/s or more is forecasted, it is defined as a forecast of a strong wind event. Figure 8 shows the ROC curves for the static MOS model and the proposed dynamic models, with GPV-GSM and GPV-MSM as the input NWPs. Regardless of the model used, the true positive rate can be increased by increasing the model parameter $\gamma$, but this also causes an increase in the false positive rate. When GPV-MSM is used as the input NWP, the ROC curve is the closest to the upper left corner, indicating that the true positive rate shows the highest value, whereas false positive rate has been kept low.

Table 5. Area under the curve (AUC) of proposed model and conventional model.

Model	AUC
Static MOS model	0.84
Dynamic model with GPV-GSM	0.92
Dynamic model with GPV-MSM	0.94

shows the AUC for each model. The AUC values were improved from 0.84 in the static MOS model to 0.94 in the proposed dynamic model with GPV-MSM, since NWP with GPV-MSM increases the resolution of land use and topography, resulting in the improvement of NWP performance, as mentioned by Kikuchi et al. [22]. Moreover, it should be noted that NWPs with very high resolution may not improve the performance of predictions, as shown in Goit et al. [23], in which the resolution of the NWP was increased from 2 km to 333 m and an accuracy improvement was not observed.

3.4. Optimal Quantile Level of the Maximum Instantaneous Wind Speed

In order to find an optimal quantile level of the maximum instantaneous wind speed for engineering applications, the cost function of Pinson et al. [24] was used. The total loss caused by the false negative $c$ and the false positive $b$ is considered as the cost function. Assuming that the loss caused by one false negative event is $l$ and the loss caused by one false positive is $\alpha l$, the total loss $l_{\mathrm{total}}$ over the forecast period can be written as Equation (20). In this study, by changing the value of $\alpha$, the cost function as function of $\gamma$ is investigated.

$$l_{\mathrm{total}}=cl+b\alpha l$$

(20)

Figure 9 shows the dimensionless loss $l_{\mathrm{total}}/l$ as a function of $\gamma$ when $\alpha$ is 0.5 and 1. There is an optimum value at $\gamma =0.9$ when $\alpha =1$ and at $\gamma =0.7$ when $\alpha =0.5$. Since the number of occurrences of strong wind events is limited and the number of false negatives is limited, whereas the number of false positives is not limited, the graph is asymmetrical. It is obvious that the optimal quantile level of the maximum instantaneous wind speed depends on the losses by misses and false alarms. When $\alpha =1$, i.e., the losses by misses are equal to those by false alarms, the cost decreases when using a higher quantile level of the maximum instantaneous wind speed. However, when $\alpha =0.5$, i.e., the losses by false alarms are lower than those by misses, the cost decreases when using a lower quantile level of the maximum instantaneous wind speed.

4. Conclusions

In this study, a maximum instantaneous wind speed forecast model has been proposed, based on numerical weather prediction data and on-site measurement data, and validated using data measured in a mountainous area. The conclusions obtained are as follows:

• The maximum instantaneous wind speed forecast model is constructed using the ARX model. A non-parametric regression with multiple time scale forgetting factors is proposed to identify the model parameters. The maximum instantaneous wind speeds calculated using the proposed dynamic model show favorable agreement with the measurements, whereas the conventional static MOS model underestimates them.
• The prediction accuracy of the maximum instantaneous wind speed forecast is improved when high-resolution forecast data are used as the input NWP and when multiple time scale forgetting factors are adopted for the proposed dynamic model.
• The predictability of a strong wind event with a maximum instantaneous wind speed of 15 m/s or more has been evaluated using the ROC curve and the AUC. The proposed dynamic model increases the true positive rate and the AUC increases from 0.84 in the static MOS model to 0.94 in the proposed dynamic model.