FINANCE , STATISTICS, FORECASTING. ACCURACY TESTING AND ERROR MATRICES

lable at ScienceDirect

Renewable Energy 112 (2017) 474e485

Contents lists avai

Renewable Energy

journal homepage: www.elsevier.com/locate/renene

Forecasting ground-level irradiance over short horizons: Time series, meteorological, and time-varying parameter models

Gordon Reikard a, *, Sue Ellen Haupt b, Tara Jensen b

a U.S. Cellular, USA b National Center for Atmospheric Research, USA

a r t i c l e i n f o

Article history: Received 30 December 2016 Received in revised form 28 March 2017 Accepted 5 May 2017 Available online 7 May 2017

Keywords: Solar irradiance Meteorological models Time series models Forecasting

* Corresponding author. E-mail address: [email protected] (G. Reikar

http://dx.doi.org/10.1016/j.renene.2017.05.019 0960-1481/© 2017 Elsevier Ltd. All rights reserved.

a b s t r a c t

One of the key enabling technologies for integrating solar energy into the grid is short-range forecasting. Two issues have emerged in the literature. The first has to do with the relative merits of physics-based versus time series models. The second is how to parameterize short-term variability. One promising approach is time-varying parameter models. Time series models can be updated using moving windows. Meteorological models can be adjusted to match the data more closely. This study evaluates several types of models over forecast horizons ranging from 15 min to 4 h, using data from two locations in the United States. The Weather Research Forecast (WRF) model is a state-of-the art numerical weather prediction system. The Dynamic Integrated Forecast (DICast) system combines meteorological models with statis- tical adjustments. The primary time series model is the ARIMA. Several other techniques are also tested, cloud advection, smart persistence forecasts and regression trees. Each type of model is found to have particular strengths and weaknesses. Among time series models, ARIMAs with time-varying coefficients are superior to fixed coefficient methods. In a direct comparison of meteorological and time series models, the ARIMA is more accurate at short horizons, while the numerical weather prediction models are more accurate as the horizon extends. The convergence point, at which the two methods achieve similar degrees of accuracy, is in the range of 1e3 h. Adjusting meteorological model output using statistical corrections at regular intervals, as in the DICast, consistently outperforms the alternatives at horizons of 2e4 h, and is highly competitive at 1 h.

© 2017 Elsevier Ltd. All rights reserved.

1. Introduction

One of the key enabling technologies for integrating solar en- ergy into the grid is short-term forecasting. In most utilities in North America, balancing reserves, used to buffer imbalances be- tween supply and demand, are calculated at the 1 h horizon. The Federal Energy Regulatory Commission has also mandated 15 min transmission scheduling to assist in integrating variable sources [1]. Some independent system operators are scheduling at horizons of as little as 5 min. Forecasts at somewhat longer horizons are used in operational planning, peak load matching, switching sources and planning.

An extensive literature on forecasting has arisen over the last two decades. Several classes of forecasting models have demon- strated value over particular horizons [2]. At the shortest time

d).

scales, less than 15 min, models based on sky image data can often outperform other methods [3e8]. Time series models have also been found to perform well over short horizons, ranging from a few minutes to several hours. At longer horizons, on the order of 4 h and beyond, meteorological or Numerical Weather Prediction (NWP) models have been found to yield the most accurate predictions [9e13].

Two issues have emerged in the literature. The first has to do with the relative merits of physics-based versus time series methods. In principle, meteorological models are attractive because they can capture the factors influencing ground-level irradiance: cloud cover, precipitation, humidity and aerosols. The corre- sponding disadvantage is that these variables may be difficult to predict accurately. In time series models, the causal factors are captured only implicitly, through lag coefficients. Despite this, they often predict more accurately at short horizons.

The second issue is how to model short-term variability. The atmosphere is known to have multifractal properties: It is charac- terized by high degrees of intermittency and irregular outliers

G. Reikard et al. / Renewable Energy 112 (2017) 474e485 475

[14,15]. One approach is time varying parameter models. Evidence from other fields, primarily econometrics, has demonstrated that stochastic parameter models can often predict more accurately than their fixed-coefficient counterparts [16,17]. Meteorological models can be made time varying by adjusting the forecasts to more closely match the recent data.

This study evaluates meteorological and time series models over horizons ranging from 15 min to 4 h. The models are reviewed in Section 2. The databases and forecasting methodology are reviewed in Section 3. The empirical findings are reported in Section 4. Further analysis is conducted in Sections 5e6. Section 7 concludes.

2. The forecasting models

Table 1 provides a glossary of the forecasting models. The meteorological forecasts were generated using systems developed at the National Center for Atmospheric Research (NCAR). These include two versions of the Weather Research Forecast model and the Dynamic Integrated Forecast system. Two cloud advection models are also tested. The primary time series model tested here is the well-known autoregressive, integrated, moving average (ARIMA) class. The other statistical models include persistence and regression forecasts based on the clearness index, i.e., the ratio of the irradiance reaching the surface of the earth to irradiance impinging on the top of the atmosphere. The DICast and ARIMA incorporate time-varying parameters. The other techniques are fixed coefficient models.

2.1. Meteorological and combined models

The Weather Research and Forecasting (WRF) model uses the Navier Stokes equations of fluid flow [18,19]. It simulates advection in the atmosphere using the initial conditions established by ob- servations and boundary conditions from a global model, with a series of parameterizations for the unresolved and external pro- cesses [20].

A recent version known as WRF-Solar incorporates several new features that enhance its ability to predict irradiance [21]. The solar tracking algorithm accounts for changes in the eccentricity of the Earth's orbit and the obliquity of the Earth's axis. The model output includes diffuse and direct normal as well as global horizontal irradiance. Interpolation algorithms are used to account for

Table 1 Glossary of the forecasting models.

Model Description

Numerical Weather Prediction models WRF-Solar-NOW Weather Research and Forecasting model. Now-casting version

calculated using the National Centers for Environmental Predic WRF-Solar-DA Weather Research and Forecasting model. Day ahead, 3 km reso

NCEP's High Resolution Rapid Refresh system for the first 15 h DICast Dynamic Integrated Forecast system, using statistical adjustme

model forecasts are post-processed to correct the bias using the forecasts are combined using optimized weights.

Cloud Advection Models CIRACast Combines cloud fields identified by satellite observations with MADCast Multisensor Advection Diffusion algorithm. Combined Models Weighted Average Weighted average of the WRF-Solar-Now, CIRACast and MADC Time Series Models ARIMA Autoregressive, integrated, moving average model, with time-v

Coefficient variation estimated using a moving window. Transfer Function ARIMA with the clearness index as an input. The clearness ind Smart Persistence The clearness index is set equal to its previous value, while the s

and day of year. Regression tree The clearness index, forecasted using a regression tree.

irradiance between model runs, and a fast radiative transfer algo- rithm calculates the surface irradiance. A new parameterization is used to capture absorption and scattering of radiation by aerosols. Three dimensional aerosols are allowed to interact with the cloud microphysics. The model also accounts for the feedbacks in short- wave irradiance from smaller clouds.

The forecasts are from two different versions of WRF-Solar. The first is the now-casting version of WRF-Solar (hereafter WRF-Solar- Now), which was run hourly at a 9 km resolution over the contig- uous United States. The initial and boundary conditions are calcu- lated using the National Centers for Environmental Prediction (NCEP) Rapid Refresh model. The second version, WRF-Solar Day Ahead (hereafter WRF-Solar-DA), uses a higher resolution of 3 km. It is run once per day, targeting day-ahead decision points. This version is initialized using NCEP's High Resolution Rapid Refresh system for the first 15 h.

Prior studies have found that models combining atmospheric physics and statistical adjustment can predict more accurately than meteorological models alone, over horizons of 1e5 h [22e26]. The forecasts are from the Dynamic Integrated forecast (DICast®) sys- tem, which adjusts NWP model output based on recent data [27]. Forecasts have been created for several measures of irradiance at selected locations where observational data was provided by pri- vate utilities [13,28]. DICast uses a two-step process. The NWP model forecasts are post-processed using a Model Output Statistics (MOS) approach [29]: the model output is adjusted upward or downward to match the most recent actual values. Then, several of these adjusted forecasts are combined using optimized weights. In effect, DICast is an artificial intelligence system that continually updates the forecasts so as to more closely approximate the observed data.

2.2. Cloud advection models

Two cloud advection models are also tested. Analyses of satellite-derived advection schemes have found that these fill a gap between short-term sky-imaging methods and methods based on meteorological models alone, as well as capturing ramps in irra- diance [30,31]. Empirical tests of satellite-derived advection fore- casts have concluded that these are effective primarily within very short time frames [32e42].

The CIRACast algorithm identifies cloud fields from satellite

Temporal Resolution

, 9 km resolution. Initial and boundary conditions tion (NCEP) Rapid Refresh model.

15 min, 1 h

lution. Initial and boundary conditions calculated using , transitioning to the Rapid Refresh at longer horizons.

1 h

nt and blending of multiple NWP model outputs. The recent forecast error. Then, several of the adjusted NWP

1 h

wind forecasts from the NOAA Global Forecast System. 15 min, 1 h 15 min, 1 h

ast. 15 min, 1 h

arying coefficients. Specification: ARIMA (1,0,0)(1,1,0). 15 min, 1 h

ex is forecasted using a regression on proximate lags. 15 min, 1 h olar angle changes In relation to the specific time of day 15 min, 1 h

15 min, 1 h

G. Reikard et al. / Renewable Energy 112 (2017) 474e485476

observations [43e45]. The Colorado State University's Cooperative Institute for Research in the Atmosphere (CIRA) team has devel- oped an operational satellite-derived forecast for predicting irra- diance over intervals of 15 min to 3 h. The CIRA method utilizes real-time cloud imagery obtained from the Pathfinder Atmospheres Extended (PTMOS-x) retrieval suite, which is based on geosta- tionary satellite observations [46]. The cloud motion is then pre- dicted using wind forecasts derived from the Global Forecast System (GFS). Ground level irradiance is computed using a radiative transfer model.

The second forecast uses the Multisensor Advection Diffusion (MADCast) algorithm of Aulign�e [47,48]. The method is in three stages. The first is to retrieve observations and express them as cloud fractions. Irradiance is calculated using a radiative transfer model, under the assumption of clear skies. Departures between the irradiance measured by the satellites and irradiance estimated by the model are computed. The cloud fraction is based on the difference between the two. The cloud fraction profiles are then interpolated to the model grid points. The interpolations are enhanced using the information from multiple satellite platforms, to obtain better estimates of their horizontal and vertical resolu- tions. The second stage is to forecast the cloud fraction, using the dynamic core of the WRF model. The third stage is to convert these forecasts into ground-level irradiance.

2.3. Time series models

Time series models can be implemented using commercially available software, and the programming is straightforward. The ARIMA class of models is well-established [49]. Using standard notation, let f(L) be the autoregressive operator, represented as a polynomial in the backshift operator: f(L) ¼ 1 e f1L - … - fpLp. Let F(L) be the cyclical autoregressive operator, defined the same way. Let q(L) be the moving average operator: q(L) ¼ 1 þ q1L þ … þ qqLq, and Q(L) be the cyclical moving average operator. Let the super- script x denote the order of differencing, and the superscript z denote the order of cyclical differencing. Let the superscript f denote the cyclical frequency, for the hourly data, 24 h. Let Yt denote ground-level irradiance and the t-subscript denote time variation. The model is then of the form:

(1-L)x(1-Lf)zYt ¼ [qt(L) Qt(L) / ft(L) Ft(L)] εt (1)

where ε is the residual and the coefficients are stochastic. The success of the ARIMA and related models in forecasting

irradiance traces back to their ability to reproduce the diurnal cycle. An early review of time series techniques found that ARIMA models were able to outperform most of the alternatives, including neural networks and unobserved components models, which use trigo- nometric terms to capture cyclical behavior [50]. This has been confirmed in more recent studies [51].

The ARIMA can also include causal inputs, in which case it is referred to as a transfer function. Inputs that have been considered include a range of meteorological variables such as temperature, precipitation, cloud cover, and processed satellite images [52e56]. Let Mt denote the input, l and d denote the moving average and autoregressive polynomials for the input. The transfer function is of the form:

(1-L)x(1-Lf)zYt ¼ [qt (L)t Q t (L) / ft (L) Ft (L)] εt þ [lt (L) / dt (L) ] Mt (2)

Other popular techniques include neural networks [57e65]. Training neural networks directly on irradiance has been found to forecast effectively only at very short horizons. Instead, the approach used by many of these studies has been to train the net on

the clearness index. The advantage of this method is that it iden- tifies the effect of aerosols and cloud cover, and takes the solar angle into account; this is computed explicitly. The corresponding disadvantage is that cloud cover can be highly intermittent, making it difficult to predict.

Two forecasts using the clearness index are also tested. These were run at NCAR as part of the Sun4Cast solar forecasting system [13], and are reproduced here. The first is simply a “smart” persistence forecast, which is often used as the baseline for the other forecasts to beat: The clearness index is assumed to be equal to its previous value, while the solar angle changes as computed for the specific time of day and day of year. The second is a regression tree forecast based on the open source Cubist software [53,66]. These forecasts were run hourly to predict 15 min interval forecasts out to 3 h.

3. The databases and the forecasting methodology

3.1. The data

Seven data sets at two locations in the United States were used in the forecasting tests. Table 2 shows the datasets, the interval spanned, and the number of usable observations, i.e., observations for which there are meteorological model forecasts, less missing and nighttime values. The data are available at two resolutions, 15- min and hourly. All the databases include global horizontal irradi- ance (GHI) and the clearness index. Irradiance is denominated in watts per meter squared (W/m2). The clearness index ranges be- tween values of 0 and 1, where 1 indicates completely clear skies.

Figs. 1e2 show the locations on maps. Four sites are from the Sacramento Municipal Utility District (SMUD) in California (latitude 38.33 N, longitude 121.28 W). The district spans about 1400 square kilometers. Site 67 lies east of Sacramento. Sites 68e69 are north- east of the city, near the towns of Folsom and Roseville respectively. Site 70 is further south, and west of the Interstate 5 highway. At the 15-min resolution, the data run from January 2, 2015 through April 28, 2016, and consist of 17,701 to 19,636 usable values. At the hourly resolution, the data run from January 17, 2015 through June 6, 2016. The number of usable observations ranges from 3056 to 5232.

Three sites are at Brookhaven National Laboratory, in Upton, New York (latitude 40.52 N, longitude 72.53 W). The Long Island Solar Farm (LISF) is a 32-MW solar photovoltaic plant built as a joint venture between the Department of Energy and the Long Island Power Authority. Since the power plant spans only about 200 acres, these sites are closer together. The data run from February 6, 2015 through April 12, 2016. The number of usable values ranges from 3659 to 3815 for the hourly data, and 15,170 to 15,808 for the 15- min data.

Figs. 3e4 show ground level irradiance at the 24-h resolution, at Sacramento site 69, and Brookhaven site 13. Irradiance is domi- nated by the diurnal cycle, but also shows high degrees of nonlinear variability. The data at Brookhaven is much more volatile, due to higher precipitation and cloud cover.

In all the data sets, there were large numbers of missing ob- servations. Missing values are not a problem for the meteorological models, since all of the forecasts correspond with actual observa- tions. However, they are an issue for the ARIMAs, which need to be estimated over continuous data streams. Several interpolation methods were tried. The optimal method for Sacramento was to estimate an ARIMA over the actual observations, and use the fitted values to fill in the missing data points. At Brookhaven, this was impossible, since with more missing observations, the likelihood function could not be maximized. Instead, the procedure was to use data from Brookhaven's solar bay station, where there were no missing values. The bay station is located less than 1 km from the

Table 2 The data.

Site Start date End date Usable Observations

Sacramento 15-min resolution Site 67 January 2, 2015, 02:00 h April 28, 2016, 23:00 h 19441 Site 68 January 2, 2015, 02:00 h April 28, 2016, 23:00 h 19523 Site 69 January 2, 2015, 02:00 h April 28, 2016, 23:00 h 19536 Site 70 January 2, 2015, 02:00 h April 28, 2016, 23:00 h 17701 Hourly resolution Site 67 January 17, 2015, 20:00 h June 6, 2016, 0:100 h 5232 Site 68 January 17, 2015, 20:00 h June 6, 2016, 0:100 h 5305 Site 69 January 17, 2015, 20:00 h June 6, 2016, 0:100 h 5215 Site 70 January 17, 2015, 20:00 h April 14, 2016, 10:00 h 3056 Brookhaven 15-min resolution Site 13 February 7, 2015, 01:00 h April 12, 2016, 14:00 h 15714 Site 18 February 7, 2015, 01:00 h April 12, 2016, 14:00 h 15170 Site 24 February 7, 2015, 01:00 h April 12, 2016, 14:00 h 15808 Hourly resolution Site 13 February 6, 2015, 20:00 h April 12, 2016, 14:00 h 3797 Site 18 February 6, 2015, 20:00 h April 12, 2016, 14:00 h 3659 Site 24 February 6, 2015, 20:00 h April 12, 2016, 14:00 h 3815

Fig. 1. Sites 67e70, Sacramento Municipal Utility District (SMUD).

G. Reikard et al. / Renewable Energy 112 (2017) 474e485 477

power plant. The observations for irradiance were similar to those at the three sites, making these values suitable to use as

interpolations. In the tests of forecast accuracy, however, all the interpolated

Fig. 2. The long Island solar Farm, at Brookhaven National Laboratory, New York.

Fig. 3. Global Horizontal Irradiance at Sacramento site 69. Left scale: watts per meter squared. Resolution: Hourly. Time span: May 1 to May 31, 2015. Source: Sacramento Municipal District and National Center for Atmospheric Research.

Fig. 4. Global Horizontal Irradiance at Brookhaven site 13. Left scale: watts per meter squared. Resolution: Hourly. Time span: May 1 to May 31, 2015. Source: Brookhaven National Laboratory.

G. Reikard et al. / Renewable Energy 112 (2017) 474e485478

values were omitted. To insure comparability, the time series forecasts are evaluated only for the same intervals as for the meteorological models.

3.2. Design of the tests

The forecasting experiments for Sacramento were set up as follows. For the 15-min data, the following models were used: the

smart persistence and regression tree, the WRF-Solar-Now, the ARIMA, the transfer function, the CIRACast and MADCast, and a weighted average of the WRF-Solar-Now and cloud advection models. For the hourly data, the models also include the WRF- Solar-DA, and the DICast. Two measures of forecast accuracy are used, the mean absolute error (MAE), in W/m2, and the root mean squared error (RMSE). The RMSE assigns a stronger penalty to large errors.

Table 3 Model Accuracy, 15 min resolution, Sacramento.

Site Forecast Horizon

15 min 30 min 45 min

Part 1: The mean absolute error (watts per meter squared)

Sacramento 67 Smart Persistence 86.5 99.8 106.2 Regression tree 78.1 85.7 94.2 NWP-Solar-NOW 67.6 67.8 66.5 CIRACast 70.7 71.4 70.9 MADCast 75.6 81.9 84.5 Weighted Average 77.7 86.6 88.6 ARIMA 42.4 54.8 61.5 Transfer Function 42.9 57.2 65.2 Sacramento 68 Smart Persistence 85.9 98.5 111.8 Regression tree 73.6 80.7 88.3 NWP-Solar-NOW 83.4 80.6 80.4 CIRACast 84.8 85.6 84.6 MADCast 70.8 74.5 80.8 Weighted Average 71.1 84.7 91.2 ARIMA 41.1 54.3 60.7 Transfer Function 42.7 57.6 65.6 Sacramento 69 Smart Persistence 81.6 94.8 105.7 Regression tree 68.5 75.2 83.7 NWP-Solar-NOW 68.1 68.6 68.5 CIRACast 71.9 72.8 73.2 MADCast 80.9 83.1 75.5 Weighted Average 72.5 80.1 85.7 ARIMA 41.4 54.2 60.9 Transfer Function 42.7 56.9 64.8 Sacramento 70 Smart Persistence 79.1 92.2 104.9 Regression tree 67.6 72.3 81.6 NWP-Solar-NOW 65.1 64.3 64.5 CIRACast 74.3 75.5 77.3 MADCast 60.6 65.6 69.5 Weighted Average 70.8 79.4 86.9 ARIMA 42.2 53.1 60.7 Transfer Function 43.4 56.9 65.1 Average, Four sites Smart Persistence 83.3 96.3 107.2 Regression tree 72.0 78.5 87.0 NWP-Solar-NOW 71.1 70.3 70.0 CIRACast 75.4 76.3 76.5 MADCast 72.0 76.3 77.6 Weighted Average 73.0 82.7 88.1 ARIMA 41.8 54.1 61.0 Transfer Function 42.9 57.2 65.2

Part 2: The root mean squared error

Sacramento 67 Smart Persistence 126.5 152.9 157.5 Regression tree 146.7 158.4 170.2 NWP-Solar-NOW 159.9 155.5 148.7 CIRACast 130.1 155.2 145.8 MADCast 124.6 123.8 124.8 Weighted Average 117.4 136.8 130.1 ARIMA 72.4 92.7 103.2 Transfer Function 72.9 91.5 100.9 Sacramento 68 Smart Persistence 128.1 142.8 156.8 Regression tree 133.4 139.8 150.7 NWP-Solar-NOW 169.4 172.4 169.1 CIRACast 135.2 140.1 166.1 MADCast 115.7 116.2 115.4 Weighted Average 119.3 125.1 126.7 ARIMA 71.8 90.5 99.6 Transfer Function 72.4 92.1 102.7 Sacramento 69 Smart Persistence 121.1 143.7 158.3 Regression tree 137.1 150.4 265.4 NWP-Solar-NOW 140.1 149.8 161.1 CIRACast 146.6 156.7 168.9

(continued on next page)

G. Reikard et al. / Renewable Energy 112 (2017) 474e485 479

The ARIMA models were specified as ARIMA (1,0,0)(1,1,0), i.e., the model is differenced at the cyclical horizon; it includes one proximate lag and one lag corresponding to the diurnal cycle. For the hourly data, the interval of differencing is of course 24 h. For the 15-min data, the interval of differencing is 96 periods. Various specifications were essayed, including longer lags. However, the simpler specification produced the lowest forecast errors.

In the transfer function, the clearness index is used as an input. The issues involved in estimating the transfer function were com- plex. First, missing values in the clearness index had to be inter- polated. A battery of interpolation algorithms was run; the best results were obtained using a regression with stochastic co- efficients. Spurious interpolations such as negative values or values in excess of unity were constrained to lie between 0 and 1. Second, the clearness index itself needs to be forecasted. Several methods were tried, but none were able to predict effectively beyond very short horizons. The model used here was a regression on lags. A neural network produced very similar results.

In the time series models for the 15-min data, the first 2000 observations were used as a training sample. The irradiance and clearness index series were then forecasted iteratively, over hori- zons of 15, 30 and 45 min. In the tests for the hourly data, the first 500 observations were used as the training sample, and the fore- casts were run over horizons of 1e4 h. In each instance, the models were estimated over prior values, forecasted, then re-estimated over the most recent value and forecasted again, until the end of the time series. All the predictions are true out-of-sample forecasts, in that they only use data prior to the start of the horizon. The forecasts for horizons beyond one observation are for this interval only, and skip the intervening values.

Time-varying parameter regressions can be estimated either using a Kalman filter [67] or a moving window. With an unre- stricted Kalman filter, the coefficients behave as a random walk, reducing predictive accuracy, so the moving window was used instead. Narrower widths allow high degrees of coefficient varia- tion, while wider widths make the coefficients more inertial [68]. Preliminary tests were run over a range of moving windows. For the 15-min data, the lowest errors were found over window width of 1400e1700 observations. The width used in the tests was 1580 observations (365 h or roughly 15 days). For the hourly data, the smallest errors were found for widths in the range of 400e600 h. In the tests, a window width of 480 h (20 days) was used.

4. The forecasting tests

Because of the number of models tested, it is useful to report the salient findings up front. At horizons of 15e45 min, the ARIMA is generally superior. At horizons of 2e4 h, the DICast achieves the most accurate forecasts. At 1 h, the contest between the ARIMA and the DICast is close at the California sites, but the DICast is superior at Brookhaven. The ARIMA is generally more accurate than WRF- Solar models at short horizons, but as the horizon extends, the meteorological models achieve greater accuracy.

4.1. Sacramento

Table 3 shows the results for the 15-min data at Sacramento. Parts 1 and 2 report the MAE and RMSE respectively. At this reso- lution, the ARIMA easily achieves the most accurate forecasts in terms of the MAE. The average error for the ARIMA is 41.1 W/m2 at the 15-min horizon, increasing to 53.6 W/m2 at 30 min and 60.2 W/ m2 at 45 min. The error for the transfer function is slightly higher. At first sight, this might appear counterintuitive: including the clearness index should make the model more sensitive to cloud cover. However, including more terms on the right hand side can

Table 3 (continued )

Site Forecast Horizon

15 min 30 min 45 min

MADCast 119.2 123.3 124 Weighted Average 111.2 124.8 129.7 ARIMA 70.2 88.3 98.1 Transfer Function 70.8 89.9 100.7 Sacramento 70 Smart Persistence 123.8 145.5 167.4 Regression tree 127.7 134.6 145.2 NWP-Solar-NOW 170.1 167.2 165.7 CIRACast 129.2 140.4 163.1 MADCast 129.5 123.2 130.4 Weighted Average 114.9 126.9 136.2 ARIMA 72.0 89.3 99.1 Transfer Function 72.5 90.8 110.7 Average, Four sites Smart Persistence 124.9 146.2 160.0 Regression tree 136.2 145.8 182.9 NWP-Solar-NOW 159.9 161.2 161.2 CIRACast 135.3 148.1 161.0 MADCast 122.3 121.6 123.7 Weighted Average 115.7 128.4 130.7 ARIMA 71.6 90.2 100.0 Transfer Function 72.1 91.1 103.8

Table 4 Model Accuracy, hourly resolution, Sacramento.

Site Forecast Horizon

1 h 2 h 3 h 4 h

Part 1: Mean absolute error (watts per meter squared)

Sacramento 67 Smart Persistence 120.4 153.9 182.7 195.4 Regression tree 103.9 130.5 e e CIRACast 88.5 98.5 e e MADCast 72.1 73.8 75.6 78.1 Weighted average 89.9 91.7 99.1 90.3 NWP-Solar-NOW 67.3 65.6 64.8 64.3 NWP-Solar-DA e e e e DICast 38.9 43.1 44.5 44.6 ARIMA 46.6 61.1 65.5 68.3 Transfer Function 52.4 70.6 75.5 80.2 Sacramento 68 Smart Persistence 125.1 158.8 186.4 198.6 Regression tree 100.1 128.2 e e CIRACast 85.0 91.2 e e MADCast 87.9 89.2 90.1 93.3 Weighted average 95.1 104.8 112.7 106.7 NWP-Solar-NOW 66.1 81.7 80.8 81.7 NWP-Solar-DA 66.0 65.6 66.6 67.1 DICast 49.1 66.2 66.1 66.8 ARIMA 47.4 62.1 68.1 71.8 Transfer Function 53.4 70.6 79.4 84.4 Sacramento 69 Smart Persistence 130.5 157.4 187.2 201.9 Regression tree 93.3 122.3 e e CIRACast 90.1 95.8 e e MADCast 71.7 73.2 75.4 77.9 Weighted average 87.4 92.4 100.4 93.2 NWP-Solar-NOW 67.9 67.1 68.1 66.1 NWP-Solar-DA 49.8 49.1 49.9 50.1 DICast 36.2 43.2 44.1 49.3 ARIMA 48.5 62.8 69.1 72.1 Transfer Function 54.1 70.1 78.3 83.4 Sacramento70 Smart Persistence 117.8 156.3 186.8 202.1 Regression tree 89.4 113.3 e e CIRACast 74.9 77.6 e e MADCast 75.2 76.1 77.1 77.6 Weighted average 90.7 95.4 103.6 92.7 NWP-Solar-NOW 67.3 63.3 63.1 63.2 NWP-Solar-DA 63.5 61.5 67.5 67.1 DICast 44.5 58.9 58.6 58.9 ARIMA 43.3 55.9 60.6 62.7 Transfer Function 49.9 66.9 67.4 68.8 Average, Four sites Smart Persistence 123.5 156.6 185.8 199.5 Regression tree 96.7 123.6 e e CIRACast 84.6 90.8 e e MADCast 76.7 78.1 79.6 81.7 Weighted average 90.8 96.1 104.0 95.7 NWP-Solar-NOW 67.2 69.4 69.2 68.8 NWP-Solar-DA 59.8 58.7 61.3 61.4 DICast 42.2 52.9 53.3 54.9 ARIMA 46.4 60.5 65.8 68.7 Transfer Function 52.5 69.6 75.2 79.2

Part 2: Root mean squared error

Sacramento 67 Smart Persistence 163.2 199.1 230.6 247.4 Regression tree 145.2 174.4 e e CIRACast 128.1 144.9 e e MADCast 114.1 116.6 119.2 119.9 Weighted average 136.2 139.9 151.8 152.9 NWP-Solar-NOW 114.3 111.6 116.4 119.7 NWP-Solar-DA e e e e DICast 73.8 74.9 75.4 76.1 ARIMA 76.8 96.7 106.1 112.5 Transfer Function 90.7 113.1 123.8 132.5 Sacramento 68 Smart Persistence 177.8 211.4 241.2 256.3 Regression tree 147.7 177.1 e e

G. Reikard et al. / Renewable Energy 112 (2017) 474e485480

cause the model to become “over-parameterized”. In effect, the clearness index adds a second term for behavior that is already captured by the lags. As a result, the model terms are not inde- pendent and may interfere with each other, reducing predictive accuracy.

At the 15-min horizon, the results from several of the other models are similar. The MAEs for the regression tree, WRF-Solar and CIRACast all lie in a range of 74.4e77.3 W/m2. However, at 30 and 45 min, the WRF-Solar-Now model is clearly superior. The MADCast model generates lower MAEs at 15 min, although at 45 min the two cloud advection models achieve similar degrees of accuracy. At all horizons, these models are consistently able to beat the smart persistence forecast. The regression tree is somewhat better, but the errors remain prohibitively high.

The findings for the RMSE confirm that the ARIMA and transfer function are more accurate at short horizons. When this measure is used, the transfer function is more competitive. Taking the average of all four sites, the transfer function error is only negligibly higher than the ARIMA at 15e30 min, although at 45 min the ARIMA is clearly better. At Site 67, the transfer function is actually slightly better at 30e45 min, but this is an anomalous result. The more typical finding is that the ARIMA and transfer function are close, with the ARIMA slightly more accurate.

At 15 min, the weighted average of three models runs in third place, while the MADCast and persistence forecasts run fourth and fifth. The MADCast model continues to do reasonably well at 30e45 min, when it is substantially better than the CIRACast. The accuracy of the persistence forecast falls away very quickly.

Table 4 shows the results for the hourly data. Using the MAE, at the 1 h horizon, the contest between the DICast and the ARIMA is fairly close, with DICast winning in two cases, while the ARIMA is more accurate in the other two. Averaging the four sites, the mean absolute error for the DICast is 42.8 W/m2, while the mean absolute error for the ARIMA is 46.5 W/m2. The other models all show much higher errors. The WRF-Solar-Now shows an MAE of 70.5 W/m2. The WRF-Solar-DA is more accurate: the MAE is 67.5 W/m2. The cloud advection models achieve MAEs in the 83e86 W/m2 range; the MADCast is slightly better than the CIRACast. The transfer function MAE averages 51.6 W/m2, higher than the ARIMA. Again, the smart persistence and regression tree do poorly, consistently showing the highest errors.

Table 4 (continued )

Site Forecast Horizon

1 h 2 h 3 h 4 h

CIRACast 129.8 140.6 e e MADCast 145.5 146.4 148.2 149.6 Weighted average 147.7 162.3 174.2 181.2 NWP-Solar-NOW 131.2 140.5 148.8 148.7 NWP-Solar-DA 131.1 132.3 135.6 138.3 DICast 90.2 129.1 129.7 137.8 ARIMA 77.1 101.4 113.8 120.2 Transfer Function 87.5 114.8 129.4 137.6 Sacramento 69 Smart Persistence 130.2 203.6 236.1 255.2 Regression tree 136.2 169.1 e e CIRACast 129.5 142.7 e e MADCast 112.8 114.6 116.9 117.8 Weighted average 121.9 139.3 153.2 155.1 NWP-Solar-NOW 91.2 110.3 108.1 105.3 NWP-Solar-DA 80.7 83.9 88.5 91.5 DICast 62.9 63.8 73.9 79.3 ARIMA 77.3 100.7 113.8 118.0 Transfer Function 90.8 117.4 129.4 139.3 Sacramento70 Smart Persistence 180.9 221.3 251.4 267.3 Regression tree 159.3 187.2 e e CIRACast 131.9 140.1 e e MADCast 161.2 160.7 160.4 161.5 Weighted average 152.1 174.9 185.9 185.9 NWP-Solar-NOW 143.3 153.7 150.8 146.3 NWP-Solar-DA 156.1 143.6 143.7 143.7 DICast 120.8 121.6 123.6 124.5 ARIMA 74.6 96.5 104.5 113.2 Transfer Function 82.2 105.2 118.2 126.3 Average, Four sites Smart Persistence 163.0 208.9 239.8 256.6 Regression tree 147.1 177.0 e e CIRACast 129.8 142.1 e e MADCast 133.4 134.6 136.2 137.2 Weighted average 139.5 154.1 166.3 168.8 NWP-Solar-NOW 120.0 129.0 131.0 130.0 NWP-Solar-DA 122.6 119.9 122.6 124.5 DICast 86.9 97.4 100.6 104.4 ARIMA 76.5 98.8 109.6 116.0 Transfer Function 87.8 112.6 125.2 133.9

Table 5 Model Accuracy, 15-min resolution, Brookhaven.

Site Forecast Horizon

15 min 30 min 45 min

Part 1: Mean absolute error (watts per meter squared)

Brookhaven 13 Smart Persistence 87.9 100.1 109.1 NWP-Solar-NOW 84.2 85.2 86.8 CIRACast e 82.6 93.3 MADCast 95.1 97.7 96.8 Weighted average 81.8 89.6 93.2 ARIMA 66.4 89.1 101.8 Brookhaven 18 Smart Persistence 81.6 93.9 101.1 NWP-Solar-NOW 81.5 83.4 84.3 CIRACast e 83.5 93.8 MADCast 94.2 96.1 95.2 Weighted average 79.1 86.4 89.4 ARIMA 65.1 88.4 100.9 Brookhaven 24 Smart Persistence 87.2 93.5 101.1 NWP-Solar-NOW 81.3 82.1 82.9 CIRACast e 84.4 95.1 MADCast 90.8 91.7 90.9 Weighted average 78.7 85.1 87.7 ARIMA 64.9 88.2 100.5 Average of all sites Smart Persistence 85.6 95.8 103.8 NWP-Solar-NOW 82.3 83.6 84.7 CIRACast e 83.5 94.1 MADCast 93.4 95.2 94.3 Weighted average 79.9 87.0 90.1 ARIMA 65.5 88.6 101.1

Part 2: Root mean squared error

Brookhaven 13 Smart Persistence 132.9 147.7 156.5 NWP-Solar-NOW 126.5 135.1 134.7 CIRACast e 124.8 149.3 MADCast 155.4 155.8 155.7 Weighted average 126.5 134.1 135.7 ARIMA 103.7 136.1 139.6 Brookhaven 18 Smart Persistence 125.2 142.6 149.5 NWP-Solar-NOW 138.1 142.2 143.1 CIRACast e 125.3 149.8 MADCast 154.6 154.9 155.1 Weighted average 125.5 132.2 132.9 ARIMA 102.6 134.1 138.5 Brookhaven 24 Smart Persistence 126.2 140.9 148.4 NWP-Solar-NOW 137.6 139.6 139.9 CIRACast e 127.3 152.9 MADCast 150.7 150.2 150.8 Weighted average 124.6 129.6 130.1 ARIMA 102.3 133.1 138.3 Average of all sites Smart Persistence 128.1 143.7 151.5 NWP-Solar-NOW 134.1 139.0 139.2 CIRACast e 125.8 150.7 MADCast 153.6 153.6 153.9 Weighted average 125.5 132.0 132.9 ARIMA 102.9 134.4 138.8

G. Reikard et al. / Renewable Energy 112 (2017) 474e485 481

At 2 h, the DICast comes in first, with an error of 55 W/m2. The ARIMA runs second, with an MAE of 60.3 W/m2. The WRF-Solar-DA runs third, and the WRF-Solar-Now runs fourth.

At 3e4 h, the DICast wins unambiguously, with MAEs in the range of 55e57 W/m2, while the ARIMA MAE increases to 65.8 W/ m2 and 68.7 W/m2. At the 3 h horizon, the WRF-Solar-DA and the ARIMA are tied almost exactly. At 4 h however, the WRF-Solar-DA is slightly better.

There is of course some variation at the individual sites. For instance, at Site 68, the results for the WRF-Solar-DA, ARIMA and DICast are very similar at 2e3 h. By comparison, at Site 69, the WRF-Solar-DA is better at all horizons beyond 1 h, while at Site 70, the ARIMA is better for the first 2 h.

The finding for the RMSE, reported in Part II of Table 4, differ in certain respects, but on the whole produce similar findings. The contest is again primarily between the DICast and the ARIMA. Using an average of the four Sacramento sites, the ARIMA achieves the lowest error at the 1 h horizon. The convergence point, at which the two methods achieve comparable degrees of accuracy, is 2 h. The DICast is better at 3e4 h. The WRF-Solar-DA and WRF-Solar-NOW achieve comparable degrees of accuracy at 1 h, but at all other horizons, the WRF-Solar-DA is better. The MADCast is slightly better than CIRACast at 1 h, but at these horizons the cloud advection models are not particularly competitive.

4.2. Brookhaven

Tables 5 and 6 show the results for Brookhaven. Not all the models were available for this site. Instead, the methods used here are the smart persistence, the NWP-Solar-Now, CIRAcast, MADcast, the weighted average and the ARIMA. The transfer function also could not be tested, due to extended gaps in the clearness index.

At the 15 min resolution, the ARIMA achieves the highest degree of accuracy, but the errors here are considerably higher than at the

Table 6 Model Accuracy, hourly resolution, Brookhaven.

Site Forecast Horizon

1 h 2 h 3 h 4 h

Part 1: Mean absolute error (watts per meter squared)

Brookhaven 13 Smart Persistence 115.1 153.9 171.6 187.4 MADCast 98.6 99.9 102.1 104.7 NWP-Solar-NOW 88.4 87.1 89.1 88.5 NWP-Solar-DA 76.7 77.1 77.7 79.1 DICast 51.1 63.2 63.8 64.2 ARIMA 71.4 114.8 145.1 151.4 Brookhaven 18 Smart Persistence 109.6 142.6 165.1 180.1 MADCast 96.5 98.1 100.5 102.8 NWP-Solar-NOW 86.4 85.7 85.5 86.4 NWP-Solar-DA 76.1 76.3 76.9 78.2 DICast 49.6 60.3 60.7 61.4 ARIMA 70.1 112.2 124.4 149.1 Brookhaven 24 Smart Persistence 110.6 142.7 165.6 171.1 MADCast 92.1 94.4 96.9 99.1 NWP-Solar-NOW 84.6 83.9 85.8 85.9 NWP-Solar-DA 73.5 73.7 74.3 75.3 DICast 48.2 58.8 59.3 59.8 ARIMA 70.2 113.6 126.2 150.3 Average of all sites Smart Persistence 111.8 146.4 167.4 179.5 MADCast 95.7 97.5 99.8 102.2 NWP-Solar-NOW 86.5 85.6 86.8 86.9 NWP-Solar-DA 75.4 75.7 76.3 77.5 DICast 49.6 60.8 61.3 61.8 ARIMA 70.6 113.5 131.9 150.3

Part 2: Root mean squared error

Brookhaven 13 Smart Persistence 163.4 200.6 226.8 248.7 MADCast 155.1 155.5 158.1 159.4 NWP-Solar-NOW 143.4 140.8 144.1 144.3 NWP-Solar-DA 124.7 125.1 125.6 127.6 DICast 82.7 101.8 102.9 103.8 ARIMA 122.7 179.3 198.8 219.8 Brookhaven 18 Smart Persistence 159.6 196.2 222.1 243.1 MADCast 155.7 156.6 159.1 160.6 NWP-Solar-NOW 146.6 145.3 146.7 146.6 NWP-Solar-DA 124.9 125.2 125.7 127.4 DICast 81.3 99.2 99.4 100.6 ARIMA 119.1 172.4 186.5 215.7 Brookhaven 24 Smart Persistence 160.4 195.1 220.2 240.5 MADCast 151.8 152.9 156.8 157.3 NWP-Solar-NOW 142.2 139.8 144.3 144.9 NWP-Solar-DA 122.5 122.8 123.1 124.8 DICast 79.6 96.8 97.3 98.1 ARIMA 119.6 169.2 185.7 217.6 Average of all sites Smart Persistence 161.1 197.3 223.0 244.1 MADCast 154.2 155.0 158.0 159.1 NWP-Solar-NOW 144.1 142.0 145.0 145.3 NWP-Solar-DA 124.0 124.4 124.8 126.6 DICast 81.2 99.3 99.9 100.8 ARIMA 120.5 173.6 190.3 217.7

G. Reikard et al. / Renewable Energy 112 (2017) 474e485482

Sacramento sites. The weighted average achieves the second smallest error, followed by WRF-Solar-Now. As the horizon in- creases, the accuracy of all the models falls away very quickly. The ARIMA does poorly, running fourth at 30 min. At 45 min, the WRF- Solar-Now does better than the alternatives. The accuracy of CIR- ACast falls away very quickly between 30 and 45 min. The MADCast achieves similar degrees of accuracy at all three horizons.

At the hourly resolution, the models are the smart persistence, MADCast, WRF-Solar-Now, WRF-Solar-DA, DICast and ARIMA. The

DICast achieves the highest degree of accuracy at all horizons. The forecast error from the DICast increases sharply between 1 and 2 h, but then levels off. The ARIMA places second at 1 h, but again the accuracy falls away very rapidly. The WRF-Solar-DA runs in third place at 1 h and second at 2e4 h. The WRF-Solar-Now runs in fourth place at 1 h and third thereafter.

5. Findings from the forecasting tests

Among the statistical models, the ARIMA is superior to the two models based on the clearness index, and perhaps surprisingly to the transfer function as well. This outcome is at variance with other studies, so it requires some explanation. Several other studies comparing neural networks with ARIMAs have used long training periods, extended holdout periods for the forecasts, and fixed co- efficients [70]. This procedure, however, massively understates the power of the ARIMA. Allowing regression coefficients to vary over time can capture a great deal of nonlinear variability. The main caveat is that the ARIMA works better at Sacramento. At Broo- khaven, the ARIMA is much less effective, except at the shortest horizons. The higher errors from the regression tree are probably attributable to the use of fixed coefficients.

The tests also provide further evidence on how well statistical methods compare with meteorological models. At Sacramento, the convergence point between the WRF-Solar models and the ARIMA is in the range of 3 h, using the MAE. At horizons beyond 3 h, the WRF-Solar models are more accurate. At horizons of 15 min to 1 h, ARIMA models are clearly superior. At 2 h, the contest is much closer, with the ARIMA slightly better for the California data sets. Conversely, at Brookhaven, the convergence point appears to be about an hour. At all horizons beyond 1 h, the WRF-solar models are superior.

The contest between the ARIMA and the DIcast is fairly close at short horizons. At the Sacramento sites, the ARIMA is competitive at 1 h using the MAE at Sacramento, and over somewhat longer horizons using the RMSE. The DICast, however, is more accurate at all horizons at Brookhaven.

6. Regressions for the physics models

As a further gauge of the relative strengths of the meteorological and cloud advection models, the actual values were regressed on the model forecasts at 1 h, in natural logs, using a grid-search correction for serial correlation. Tables 7 and 8 present these re- sults for Sacramento and Brookhaven respectively. The coefficients are elasticities: they express the percent change in the data in response to the percent change in the model forecast. Ideally, the elasticity should equal to 1. Rho (r) is the coefficient of fractional differencing, estimated by the serial correlation correction. Lower values of rho are consistent with less structure in the residuals.

The cloud advection models do not fit the data particularly well. the CIRAcast elasticities range from 0.70 to 0.77. The adjusted R- square ranges from 0.65 to 0.71, indicating that a substantial share of the variance remains unexplained. The MADcast is only some- what better: the elasticities range from 0.71 to 0.83 at Sacramento, and 0.73 to 0.76 at Brookhaven. The weighted average also shows poor results. The elasticities are in the range of 0.68e0.75.

The WRF models yield rather good results at Sacramento. The elasticities are close to unity at three of the sites, although site 67 shows a lower value. At Brookhaven, however, the WRF elasticities are considerably lower. Similarly, the R-squares are high at three of the Sacramento sites, but lower at Brookhaven.

The DICast shows the best results. At Sacramento, the elasticity for the DICast ranges from 0.95 to 1.06, averaging out to roughly unity. At Brookhaven, the elasticities are just short of 1. The R-

Table 7 Regression coefficients, Sacramento sites.

Site CIRAcast MADcast Weighted Average NWP-Solar-NOW NWP-Solar-DA DICast

Sacramento 67 Constant 0.31 0.84 0.82 1.21 e �0.39 Elasticity 0.89 0.82 0.74 0.75 e 1.05 R-bar-square 0.65 0.86 0.57 0.81 e 0.90 Rho 0.81 0.75 0.86 0.71 e 0.11 Sacramento 68 Constant 1.45 1.04 0.89 �0.41 �0.41 0.16 Elasticity 0.71 0.78 0.75 1.03 1.02 0.95 R-bar-square 0.69 0.84 0.59 0.87 0.87 0.89 Rho 0.91 0.91 0.89 0.89 0.81 0.77 Sacramento 69 Constant 1.05 0.71 0.86 �0.14 �0.43 �0.43 Elasticity 0.77 0.83 0.72 0.96 1.05 1.06 R-bar-square 0.70 0.85 0.58 0.89 0.90 0.91 Rho 0.91 0.84 0.88 0.78 0.31 0.31 Sacramento 70 Constant 1.66 1.58 1.43 0.71 1.96 0.11 Elasticity 0.70 0.71 0.68 0.95 0.63 0.95 R-bar-square 0.66 0.81 0.58 0.87 0.74 0.88 Rho 0.92 0.95 0.91 0.91 0.93 0.61

Statistics are regression coefficients, unless otherwise indicated. Rho is the coefficient of fractional differencing.

Table 8 Regression coefficients, Brookhaven sites.

Site MADcast Weighted Average NWP-Solar-NOW DICast

Brookhaven 13 Constant 1.04 1.27 1.27 �0.05 Elasticity 0.74 0.71 0.71 0.96 R-bar-square 0.84 0.79 0.79 0.95 Rho 0.91 0.84 0.84 0.44 Brookhaven 18 Constant 0.92 1.27 1.30 �0.12 Elasticity 0.76 0.71 0.71 0.96 R-bar-square 0.84 0.79 0.80 0.96 Rho 0.91 0.83 0.83 0.31 Brookhaven 24 Constant 1.12 1.33 1.33 �0.06 Elasticity 0.73 0.71 0.71 0.97 R-bar-square 0.84 0.79 0.81 0.96 Rho 0.91 0.82 0.83 0.29

Statistics are regression coefficients, unless otherwise indicated. Rho is the coefficient of fractional differencing.

G. Reikard et al. / Renewable Energy 112 (2017) 474e485 483

squares are in the range of 0.88e0.91 at Sacramento, and 0.96 at Brookhaven. The constant terms and coefficients of fractional dif- ferencing are also lower, pointing to less serial correlation in the residuals.

7. Conclusions

Each model has particular strengths and weaknesses. Time se- ries models are able to predict more accurately at the shortest horizons primarily because the data is dominated by dependence between proximate time points. As this dependence falls away, over longer horizons, the time series models become less effective. Expressed another way, the atmosphere has moved beyond the Lagrangian integral time scale at this point, i.e., the time span over which the atmosphere “remembers” its prior state. Among time series models, ARIMAs with stochastic parameters are superior to fixed coefficient models. While the persistence and regression tree models showed much larger errors, this does not necessarily invalidate the idea of converting irradiance to an index. In this respect, there is an emerging literature on clear sky models, which quantify the impact of factors ranging from the air mass, to cloud

cover and atmospheric turbidity [71]. A detailed comparison of clear sky models with ARIMAs, both using stochastic parameters, will be the topic of a future study.

Cloud advection models are not found to fit the data closely, and their predictive accuracy is generally lower than large-scale meteorological models. The reasons most likely have to do with the larger number of causal factors taken into account in numerical weather prediction systems.

Large-scale meteorological models do better for time scales beyond about 2 h. The reason is their ability to predict the synoptic and mesoscale changes in cloud cover and weather patterns. Their main limitation is that they do not capture all the short-term variability in the data. Statistical adjustment enables the models to track the data more effectively. The DICast combines the physics embodied in meteorological models with the ability to adjust to changing atmospheric conditions.

In a direct comparison of time series and physics-based models, the ARIMA is more accurate at short horizons, while the meteoro- logical models are more accurate as the horizon extends. The convergence points lie in a range of only 1e3 h. Comparisons for wind and wave energy have found longer convergence points, in the range of 5e6 h [72e74]. The expected integral time scale for wave energy is longer than for solar. The shorter convergence points for solar may also reflect the fact that the WRF-Solar-Now model was run in nowcasting mode, updating every hour, which is seldom available for other NWP models.

When NWP models are statistically adjusted, they are able to outperform all the alternatives at horizons of 2e4 h, while they are competitive at the 1-h horizon. The more general implication is that over horizons where the data exhibits high degrees of variability, enabling models to adapt to changing conditions raises forecast accuracy.

Acknowledgements

The authors thank Sacramento Municipal Utility District and Brookhaven National Laboratory for use of their solar irradiance measurements. They also thank the modelers who produced the forecasts used as a comparison for this work and Jared Lee for the maps used in Figs. 1e2. Finally, the authors thank the reviewers, whose comments have made this a better paper.

G. Reikard et al. / Renewable Energy 112 (2017) 474e485484

References

[1] Federal Energy Regulatory Commission (FERC), Order No. 764, Federal Energy Regulatory Center, Washington, DC, 2012.

[2] A. Tuohy, J. Zack, S.E. Haupt, J. Sharp, M. Ahlstrom, S. Dise, E. Grimit, C. M€ohlren, M. Lange, M.G. Casado, J. Black, M. Marquis, Craig Collier, Solar forecasting: methods, challenges and performance, IEEE Power Energy Mag. 13 (2015) 50e59.

[3] C.W. Chow, N. Urquhart, M. Lave, A. Dominquez, J. Kleissl, J. Shields, B. Washom, Intra-hour forecasting with a total sky imager at the UC San Diego solar energy testbed, Sol. Energy 85 (2011) 2881e2893.

[4] R. Marquez, C.F.M. Coimbra, Intra-hour DNI forecasting based on cloud tracking image analysis, Sol. Energy 91 (2013) 327e336.

[5] R. Marquez, H.T.C. Pedro, C.F.M. Coimbra, Hybrid solar forecasting method uses satellite imaging and ground telemetry as inputs to ANNs, Sol. Energy 92 (2013) 176e188.

[6] S. Quesada-Ruiz, Y. Chu, J. Tovar-Pescador, H.T.C. Pedro, C.F.M. Coimbra, Cloud-tracking methodology for intra-hour DNI forecasting, Sol. Energy 102 (2014) 267e275.

[7] Y. Chu, H.T.C. Pedro, C.F.M. Coimbra, Hybrid intra-hour DNI forecasts with sky image processing enhanced by stochastic learning, Sol. Energy 98 (2013) 592e603.

[8] Y. Chu, H.T.C. Pedro, M. Li, C.F.M. Coimbra, Real-time forecasting of solar irradiance ramps with smart image processing, Sol. Energy 114 (2015) 91e104.

[9] E. Lorenz, J. Kuhnert, D. Heinemann, Overview on irradiance and photovoltaic power prediction, Weather Matters Energy (2012) 429e454.

[10] J. Kleissl, Solar Energy Forecast and Resource Assessment, Academic Press, 2013.

[11] P. Mathiesen, J. Kleissl, Evolution of numerical weather prediction for intra- day solar forecasting in the continental United States, Sol. Energy 85 (2011) 967e997.

[12] J. Zhang, B.M. Hodge, S. Lu, H.F. Hamann, B. Lehman, J. Simmons, E. Campos, V. Banunarayanan, J. Black, J. Tedesco, Baseline and target values for regional and point PV power forecasts: toward improved solar forecasting, Sol. Energy 122 (2015) 804e819.

[13] S.E. Haupt, B. Kosovic, T. Jensen, J. Lee, P. Jimenez, J. Lazo, J. Cowie, T. McCandless, J. Pearson, G. Weiner, S. Alessandrini, L. Delle Monache, D. Yu, Z. Peng, D. Huang, J. Heiser, S. Yoo, P. Kalb, S. Miller, M. Rogers, L. Hinkleman, The SunCast Solar Power Forecasting System: the Results of the Public-private- academic Partnership to Advance Solar Power Forecasting, NCAR Technical- Report TN-526þSTR, 2016, http://dx.doi.org/10.5065/D6N58JR2, 307 pp.

[14] D. Schertzer, S. Lovejoy, F. Schmitt, Y. Chigirinskaya, D. Marsan, Multifractal cascade dynamics and turbulent intermittency, Fractals 5 (1997) 427e471.

[15] S. Lovejoy, D. Schertzer, The Weather and Climate: Emergent Laws and Multifractal Cascades, Cambridge University Press, Cambridge, 2013.

[16] C.W.J. Granger, Non-linear Models: where Do We Go Next e Time Varying Parameter Models? Studies in Nonlinear Dynamics and Econometrics 12: Article 1, 2008. http://www.bepress.cxom/snde/vol12/iss3/art1.

[17] D.W. Bunn, Modelling Prices in Competitive Electricity Markets, Wiley, New York, 2004.

[18] J.A. Dutton, The Ceaseless Wind: an Introduction to the Theory of Atmospheric Motion, McGraw-Hill, New York, 1976.

[19] G.K. Vallis, Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-scale Circulation, Cambridge University Press, Cambridge, UK, 2006.

[20] W.C. Skamarock, J.B. Klemp, J. Dudhia, D.O. Gill, D.M. Barker, M.G. Duda, X.- Y. Huang, W. Wang, J.G. Powers, A Description of the Advanced Research WRF Version 3. NCAR Tech. Note NCAR/TN-475þSTR, 2008, http://dx.doi.org/ 10.5065/D68S4MVH, 113 pp.

[21] P. Jimenez, J. Hacker, J. Dudhia, S.E. Haupt, J. Ruiz-Arias, C. Gueymard, G. Thompson, T. Eidhammer, A. Deng, WRF-Solar: an augmented NWP model for solar power prediction. Model description and clear sky assessment, Bull. Am. Meteorological Soc. 97 (7) (2015) 1249e1264, http://dx.doi.org/10.1175/ BAMS-D-14-00279.1.

[22] M. Bouzerdoum, A. Mellit, A. Massi Pavan, A hybrid model (SARIMA-SVM) for short-term power forecasting of a small-scale grid-connected photovoltaic plant, Sol. Energy 98 (2013) 226e235.

[23] C. Voyant, M. Muselli, C. Paoli, M.L. Nivet, Numerical Weather Prediction (NWP) and hybrid ARMA/ANN to predict global radiation, Energy 39 (1) (2014) 341e355.

[24] C. Voyant, M. Muselli, C. Paoli, M.L. Nivet, Hybrid methodology for hourly global radiation forecasting in Mediterranean area, Renew. Energy 53 (2013) 1e11.

[25] W. Myers, G. Wiener, S. Linden, S.E. Haupt, A consensus forecasting approach for improved turbine hub height wind speed predictions, in: Proceedings, Wind Power, 2011. Anaheim, CA, May 24, 2011, 2011.

[26] L. Delle Monache, F.A. Eckel, D.L. Rife, B. Nagarajan, K. Searight, Probabilistic weather predictions with an analog ensemble, Mon. Weather Rev. 141 (2013) 3489e3516, http://dx.doi.org/10.1175/MWR-D-12-00281.1.

[27] S.E. Haupt, W.P. Mahoney, Wind power forecasting, in: IEEE Spectrum, Nov 2015, 2015, pp. 46e52.

[28] S.E. Haupt, B. Kosovic, Variable generation power forecasting as a big data problem, IEEE Trans. Sustain. Energy 8 (2) (2017) 725e732.

[29] H.R. Glahn, D.A. Lowry, The use of model output statistics (MOS) in objective

weather forecasting, J. Appl. Meteorology 11 (1972) 1203e1211. [30] Z. Peng, S. Yoo, D. Yu, D. Huang, Solar irradiance forecast system based on

geostationary satellites, in: 2013 IEEE International Conference on Smart Grid Communications, 2013, pp. 708e713.

[31] Z. Peng, S. Yoo, D. Yu, D. Huang, P. Kalb, J. Heiser, 3D cloud detection and tracking for solar forecast using multiple sky imagers, in: Proceedings of the 29th Annual ACM Symposium on Applied Computing, Association for Computing Machinery 2014, 2014, pp. 512e517.

[32] R. Perez, S. Kivalov, J. Schlemmer, K. Hemker, D. Renne, T.E. Hoff, Validation of short and medium-term operational solar radiation forecasts for the U.S, Sol. Energy 84 (2010) 2161e2172.

[33] R. Perez, E. Lorenz, S. Pelland, M. Beautharnois, G. Van Knowe, K. Hemker, D. Heinemann, J. Remund, S.C. Muller, W. Traumiller, et al., Comparison of numerical weather prediction solar irradiance forecasts in the U.S., Canada and Europe, Sol. Energy 94 (2013) 305e326.

[34] H.G. Beyer, C. Costanzo, D. Heinemann, Modifications of the Heliosat pro- cedure for irradiance estimates from satellite data, Sol. Energy 56 (1996) 121e207.

[35] I. Bilionis, E.M. Constantinescu, M. Anitescu, Data-driven model for solar irradiation based on satellite observations, Sol. Energy 110 (2014) 22e38.

[36] S. Cros, O. Liandrat, N. S�ebastien, N. Schmutz, Extracting cloud motion vectors from satellite images for solar power forecasting, in: Geoscience and Remote Sensing Symposium (IGARSS), 2014 IEEE International, 2014, pp. 4123e4126.

[37] D.R. Stauffer, N.L. Seaman, Use of four-dimensional data assimilation in a limited area mesoscale model. Part I: experiments with synoptic-scale data, Mon. Weather Rev. 181 (1990) 1250e1277.

[38] D.R. Stauffer, N.L. Seaman, Multiscale four-dimensional data assimilation, J. Appl. Meteorology 33 (1994) 416e434.

[39] R. Tapakis, A.G. Charalambides, Equipment and methodologies for cloud detection and classification: a review, Sol. Energy 95 (2013) 392e430.

[40] A. Zagouras, A. Kazantzidis, E. Nikitidou, A.A. Argiriou, Determination of measuring sites for solar irradiance, based on cluster analysis of satellite- derived cloud estimations, Sol. Energy 97 (2013) 1e11.

[41] A. Hammer, D. Heinemann, E. Lorenz, B. Lückehe, Short-term forecasting of solar radiation: a statistical approach using satellite data, Sol. Energy 67 (1999) 139e150.

[42] E. Lorenz, A. Hammer, D. Heinemann, Short term forecasting of solar radiation based on satellite data. EUROSUN2004, ISES Eur. Sol. Congr. (2004) 841e848.

[43] S.D. Miller, M.A. Rogers, A.K. Heidinger, I. Laszlo, M. Sengupta, . Cloud advection schemes for short-term satellite-based insolation forecasts, in: Proceedings of the World Renewable Energy Forum, 17 May 2012, Denver, CO, 2012, pp. 1963e1967.

[44] S.D. Miller, J.M. Forsythe, P.T. Partain, J.M. Haynes, R.L. Bankert, M. Sengupta, C. Mitrescu, J.D. Hawkins, T.H. Vonder Haar, Estimating three-dimensional cloud structure via statistically blended satellite observations, J. Appl. Mete- orology Climatol. 53 (2014) 437e455.

[45] M.A. Rogers, S.D. Miller, J.M. Haynes, A.K. Heidinger, S.E. Haupt, M. Sengupta, Improvements in satellite-derived short-term insolation forecasting: statisti- cal comparisons, challenges for advection-based forecasts, and new tech- niques, in: Sixth Conference on Weather, Climate, and the New Energy Economy at the 95th American Meteorological Society Annual Meeting, Phoenix, AZ, 6 Jan 2015, 2015.

[46] A.K. Heidinger, M.J. Foster, A. Walther, X. Zhao, The pathfinder atmospheres extended AVHRR climate data set, Bull. Am. Meteorological Soc. (2013), http:// dx.doi.org/10.1175/BAMS-D-12-00246.1.

[47] T. Aulign�e, Multivariate minimum residual method for cloud retrieval. Part I: theoretical aspects and simulated observations experiments, Mon. Weather Rev. 142 (2014a) 4383e4398, http://dx.doi.org/10.1175/mWR-D-13-00172.1.

[48] T. Aulign�e, Multivariate minimum residual method for cloud retrieval. Part II: real observations experiments, Mon. Weather Rev. 142 (2014b) 4399e4415, http://dx.doi.org/10.1175/mWR-D-13-00173.1.

[49] G.E.P. Box, G.M. Jenkins, Time Series Analysis: Forecasting and Control, Holden-Day, San Francisco, 1976.

[50] G. Reikard, Predicting solar radiation at high resolutions: a comparison of time series forecasts, Sol. Energy 83 (2009) 342e349.

[51] A.M. Nobre, C.A. Severiano, S. Karthik, M. Marek Kubis, L. Zhao, F.R. Martins, E.B. Pereira, R. Rüther, T. Reindl, PV power conversion and short-term fore- casting in a tropical, densely-built environment in Singapore, Renew. Energy 94 (2016) 496e509.

[52] G. Lopez, F.J. Batlles, J. Tovar-Pescador, Selection of input parameters to model direct solar irradiance by using artificial neural networks, Energy 30 (2005) 1675e1684.

[53] T.C. McCandless, S.E. Haupt, G.S. Young, Short Term Solar Radiation Forecasts Using Weather Regime-dependent Artificial Intelligence Techniques, in: 12th Conference on Artificial Intelligence: Applications of Artificial Intelligence Methods for Energy, Atlanta, GA, American Meteorological Society, J3.5, 2014.

[54] T.C. McCandless, S.E. Haupt, G.S. Young, A regime-dependent artificial neural network technique for short-range solar irradiance forecasting, Renew. En- ergy 89 (2016) 351e359.

[55] T.C. McCandless, G.S. Young, S.E. Haupt, L.M. Hinkelman, Regime-dependent short-range solar irradiance forecasting, J. Appl. Meteorology Climatol. 55 (2016) 1599e1613.

[56] C.L. Fu, H.Y. Cheng, Predicting solar irradiance with all-sky image features via regression, Sol. Energy 97 (2013) 537e550.

[57] A. Mellit, Artificial intelligence technique for modeling and forecasting of solar

G. Reikard et al. / Renewable Energy 112 (2017) 474e485 485

radiation data: a review, Int. J. Artif. Intell. Soft Comput. 1 (2008) 52e76. [58] A. Mellit, A. Massi Pavan, V. Lughi, Short-term forecasting of power production

in a large-scale photovoltaic plant, Sol. Energy 105 (2014) 401e413. [59] G. Notton, C. Paoli, S. Vasileva, M.L. Nivet, J.-L. Canaletti, C. Cristofari, Esti-

mation of hourly global solar irradiation on tilted planes from horizontal one using artificial neural networks, Energy 39 (2012) 166e179.

[60] S. Bhardwaj, V. Sharma, S. Srivastava, O.S. Sastry, B. Bandyopadhyay, S.S. Chandel, J.R.P. Gupta, Estimation of solar radiation using a combination of Hidden Markov Model and generalized Fuzzy model, Sol. Energy 93 (2013) 43e54.

[61] M. Diagne, M. David, P. Lauret, J. Boland, N. Schmutz, Review of solar irradi- ance forecasting methods and a proposition for small-scale insular grids, Renew. Sustain. Energy Rev. 27 (2013) 65e76.

[62] R.H. Inman, H.T.C. Pedro, C.F.M. Coimbra, Solar forecasting methods for renewable energy integration, Prog. Energy Combust. Sci. 39 (2013) 535e576.

[63] E. Fernandez, F. Almonacid, N. Sarmah, P. Rodrigo, T.K. Mallick, P. Perez- Higueras, A model based on artificial neuronal network for the prediction of the maximum power of a low concentration photovoltaic module for building integration, Sol. Energy 100 (2014) 148e158.

[64] F. Almonacid, P.J. P�erez-Higueras, E.F. Fern�andez, L. Hontoria, A methodology based on dynamic artificial neural network for short-term forecasting of the power output of a PV generator, Energy Convers. Manag. 85 (2014) 389e398.

[65] H. Morf, Sunshine and cloud cover prediction based on Markov processes, Sol. Energy 110 (2014) 615e626.

[66] J.R. Quinlan, Improved use of continuous attributes in C4.5, J. Artif. Intell. Res. 4 (1996) 77e90.

[67] R.E. Kalman, A new approach to linear filtering and prediction problems. Transactions of the American Society of Mechanical Engineers,, J. Basic Eng. 83D (1960) 35e45.

[68] B. Rossi, A. Inoue, Out-of-sample forecast tests robust to the choice of window size, J. Bus. Econ. Statistics 30 (2012) 432e453.

[70] H.T.C. Pedro, C.F.M. Coimbra, Assessment of forecasting techniques for solar power prediction with no exogenous inputs, Sol. Energy 86 (2012) 2017e2028.

[71] M.J. Reno, C.W. Hansen, J.S. Stein, Global Horizontal Clear Sky Models: Implementation and Analysis, Sandia National Laboratory report SAND, 2012, pp. 2012e2389.

[72] G. Reikard, P. Pinson, J.R. Bidlor, Forecasting Ocean Wave Energy: the ECMWF wave model and time series methods, Ocean. Eng. 38 (2011) 1089e1099.

[73] A.M. Foley, P.G. Leahy, A. Marvuglia, E.J. McKeough, Current methods and advances in forecasting of wind power generation, Renew. Energy 37 (2012) 1e8.

[74] G. Reikard, B. Robertson, J.R. Bidlot, Combining wave energy with wind and solar: short-term forecasting, Renew. Energy 81 (2015) 442e456.