Toward scalable stochastic unit commitment. Part 1: load scenario generation

Abstract

Unit commitment decisions made in the day-ahead market and during subsequent reliability assessments are critically based on forecasts of load. Traditional, deterministic unit commitment is based on point or expectation-based load forecasts. In contrast, stochastic unit commitment relies on multiple load scenarios, with associated probabilities, that in aggregate capture the range of likely load time-series. The shift from point-based to scenario-based forecasting necessitates a shift in forecasting technologies, to provide accurate inputs to stochastic unit commitment. In this paper, we discuss a novel scenario generation methodology for load forecasting in stochastic unit commitment, with application to real data associated with the Independent System Operator of New England (ISO-NE). The accuracy of the expected load scenario generated using our methodology is consistent with that of point forecasting methods. The resulting sets of realistic scenarios serve as input to rigorously test the scalability of stochastic unit commitment solvers, as described in the companion paper. The scenarios generated by our method are available as an online supplement to this paper, as part of a novel, publicly available large-scale stochastic unit commitment benchmark.

Full text

Translate

Turn on search term navigation

http://crossmark.crossref.org/dialog/?doi=10.1007/s12667-015-0146-8&domain=pdf

Web End = http://crossmark.crossref.org/dialog/?doi=10.1007/s12667-015-0146-8&domain=pdf

Web End = Energy Syst (2015) 6:309329

DOI 10.1007/s12667-015-0146-8

ORIGINAL PAPER

Yonghan Feng1 Ignacio Rios2 Sarah M. Ryan1

Kai Sprkel3 Jean-Paul Watson4

Roger J.-B. Wets5 David L. Woodruff 6

Received: 30 April 2014 / Accepted: 16 March 2015 / Published online: 9 April 2015 Springer-Verlag Berlin Heidelberg 2015

Abstract Unit commitment decisions made in the day-ahead market and during subsequent reliability assessments are critically based on forecasts of load. Traditional, deterministic unit commitment is based on point or expectation-based load forecasts. In contrast, stochastic unit commitment relies on multiple load scenarios, with asso-

B Jean-Paul Watson

[email protected]

Yonghan Feng [email protected]

Ignacio Rios [email protected]

Sarah M. Ryan [email protected]

Kai Sprkel [email protected]

Roger J.-B. Wets [email protected]

David L. Woodruff [email protected]

1 Department of Industrial and Manufacturing Systems Engineering, Iowa State University,

Ames, IA, USA2 University of Chile, Santiago, Chile3 University of Duisburg-Essen, Essen, Germany4 Analytics Department, Sandia National Laboratories, Albuquerque, NM, USA5 Department of Mathematics, University of California Davis, Davis, CA, USA6 Graduate School of Management, University of California Davis, Davis, CA, USA

http://crossmark.crossref.org/dialog/?doi=10.1007/s12667-015-0146-8&domain=pdf

Web End = http://crossmark.crossref.org/dialog/?doi=10.1007/s12667-015-0146-8&domain=pdf

Web End = Toward scalable stochastic unit commitment. Part 1: load scenario generation

123

310 Y. Feng et al.

ciated probabilities, that in aggregate capture the range of likely load time-series. The shift from point-based to scenario-based forecasting necessitates a shift in forecasting technologies, to provide accurate inputs to stochastic unit commitment. In this paper, we discuss a novel scenario generation methodology for load forecasting in stochastic unit commitment, with application to real data associated with the Independent System Operator of New England (ISO-NE). The accuracy of the expected load scenario generated using our methodology is consistent with that of point forecasting methods. The resulting sets of realistic scenarios serve as input to rigorously test the scalability of stochastic unit commitment solvers, as described in the companion paper. The scenarios generated by our method are available as an online supplement to this paper, as part of a novel, publicly available large-scale stochastic unit commitment benchmark.

1 Introduction

Because operational constraints on thermal generation units require them to be committed well in advance of when they may be needed, system operators regularly solve unit commitment optimization problems for the day ahead. Typically, scheduling decisions for a day d are made on day d 1, based on forecasts of uncertain quantities such

as hourly load and renewables output; these quantities are generally aggregated across buses, for each of a systems load zones. In the context of traditional deterministic unit commitment procedures, such forecasts take the form of point or expected-value quantitiesrepresenting a single time series for each forecasted quantity. Uncertainty associated with such forecasts is addressed by maintaining a non-trivial level of generation reserves, which compensate for deviations from the predicted quantities as day d operations proceed.

In contrast, stochastic unit commitment procedures [16,24,25] assume the availability of a number of forecast scenarios, each representing a distinct potential time series of the forecasted quantities. Throughout, we use the term scenario in a narrow sense, representing a full specication of all random data required to instantiate a unit commitment problem, with associated probability of occurrence. In aggregate, the set of scenarios should represent the range of possible behaviors on day d. By explicitly representing forecast uncertainty through sets of scenarios, it should be possible to signicantly decrease generation reserve margins and consequently reduce overall system operation costs [21]. However, the need for multiple scenarios imposes fundamentally novel requirements on forecasting technologies, which have yet to be adequately addressed.

Our goal in this paper is to present approaches and data sources for generating quantiably accurate and realistic load scenarios for use in stochastic unit commitment. We focus on load, as opposed to renewables production, to tractably scope our study. The procedures described below extend to wind and solar plant production, which we will address in a future contribution. While stochastic bidding or trading [23] and market clearing formulations [2,15] for the day ahead wholesale electricity market have also been proposed, we consider the reliability unit commitment (also called resource adequacy assessment) process conducted by system operators. Consequently, we generate scenarios specifying zonal or system-wide load rather than demand bids

123

Toward scalable stochastic unit commitment 311

by load-serving entities. We demonstrate the ability of Regional Transmission Operators (RTOs) and other operator entities to generate accurate load scenarios for use with stochastic unit commitment procedures, using data that are readily available. Our experiments proceed in the context of publicly available data from the Independent System Operator of New England (ISO-NE). The resulting scenarios are then used in the companion paper [4] to rigorously test the scalability of a stochastic unit commitment solver. These scenarios and the resulting test cases are publicly available, lling a critical need for researchers investigating the scalability of stochastic unit commitment solvers. The value of stochastic unit commitment compared to deterministic or other approaches, in terms of cost savings or reliability improvement, is not addressed in either paper but is a subject for ongoing rigorous testing.

The remainder of this paper is organized as follows. Our methodology for data transformation and tting of load regression models is detailed in Sect. 2, which also includes a discussion of related literature and approaches. Preliminary results of the proposed methodology for modeling summer data in a single ISO-NE load zone were reported in [7]. This paper describes enhancements to our load modeling method by accounting for the non-weather-dependent component of demand (Sect. 2.3), transforming data to aggregate across load zones as well as days of the week (Sect. 2.4), and improving the process for segmenting days by weather forecast (Sect. 2.5). We then describe our procedures for generating load scenarios in Sect. 3. Comprehensive experimental results for all zones and all seasons of a year are presented for data associated with ISO-NE in Sect. 4. We conclude with a summary of our results in Sect. 5.

2 Load forecasting methodology

When forecasting load, the information available to operators on day d 1 includes

weather forecasts for day d, historical records of previous weather forecasts, and historical actual system loads. Historical system load data exhibit temporal patterns that vary according to season of the year, day of the week, and hour of the day. While some temporal load patterns are predictable based on knowledge of business hours and diurnal light patterns, the portion of load derived from heating and cooling (both industrial and residential) depends strongly on weather. And while numerical weather prediction models have become increasingly accurate over the past several decades, there remains signicant uncertainty associated with day-ahead weather forecasts. The challenge for system operators is to form an accurate and comprehensive picture of the day-ahead load, which not only includes point forecasts of the load in each hour, but also acknowledges the precision (or lack thereof) associated with those forecasts. To address this challenge, we introduce a novel optimization-based method to develop a stochastic model for the load on day d based on the weather forecast available on day d 1. Optimization models are solved to t functional regression models

for expected temporal load patterns (Sect. 2.3), obtain transformations between load patterns for different day types (Sect. 2.4), and approximate forecast error density functions (Sect. 2.5). We note that our approach does not assume normality or other parametric forms for error density functions, as are often used in experimental analyses of stochastic unit commitment (e.g., see [12]).

123

312 Y. Feng et al.

2.1 Background

Common methods for short-term load forecasting can be categorized as either based on articial intelligence or statistical techniques [10]. Methods from articial intelligence, such as neural networks, are widely used but do not provide probabilistic information that could be used to generate multiple probability-weighted scenarios. Among statistical approaches, which can provide the required probabilistic information, the most prevalent methods are time series and regression models. Due to limited space, we do not provide a complete review but refer the reader to recent surveys such as [6]. Instead, we now highlight samples of statistical approaches for load forecasting from the recent literature, focusing on achieved accuracy and limitations for purposes of scenario generation. Similarly, forecasting models for wholesale market prices have been assessed for their suitability in generating scenarios [13].

The weather variable most commonly used to predict load is temperature, due to its inuence on heating and cooling requirements. Other variables considered include humidity and cloud cover, although their impacts on load are much smaller than that of temperature. Humidity increases load in the summer, again due to cooling requirements. In contrast, cloud cover increases load in the winter (due to increased lighting requirements) and reduces load in the summer.

Liu et al. [14] analyze the nonlinear relationship between temperature and load, using estimates derived from a nonparametric regression method. They t a time series model to the residuals of the load-temperature regression and consider lags of 1, 24, and 168 h in their day-ahead forecasting model. Using actual historical temperature and load data obtained from a US utility, they demonstrate an out-of-sample mean absolute percent error (MAPE) of 1.2% for their 24-hour-ahead forecasts. Hong et al. [10] develop a multiple linear regression model of load that considers temperature, hour, type of day, and month as independent variables; the model additionally contains a linear term trend, and terms to capture interactions among the independent variables. Using actual weather data to predict hourly loads for a US utility over a one-year time period, they obtain an out-of-sample MAPE of 4.6%. Black [1] also uses a multiple linear regression model to examine the inuence of weather on load, but instead focuses on summer weekdays in the region served by ISO-NE. Time-of-day effects are captured through a separate regression model for each hour of the day, considering temperature, humidity, and irradiance as independent variables. The out-of-sample MAPEs yielded by these models average 23% for the whole New England region and 34% for individual subregions such as Connecticut and Southeast Massachusetts.

While hind-casting studies of the type described above (which use actual historical weather data as input) are useful for identifying factors that inuence hourly loads, they do not reect the accuracy or precision of load forecasts available in practice which necessarily rely on day-ahead weather forecasts, as opposed to actual quantities. Further, while load forecasting remains a topic of active investigation, considerable emphasis persists on identifying a single most accurate forecast trajectory [11]. For investigations into stochastic unit commitment, several approaches for generating multiple trajectories have been used. An ad hoc procedure was used in an early study to create scenarios with large increasing or decreasing load ramps [25]. A common approach is to perturb a single forecast with some stochastic error terms. In stud-

123

Toward scalable stochastic unit commitment 313

ies aimed at exploring the impacts of uncertainty on unit commitment, these errors have been sampled from a Gaussian process [3], densities derived from a jackknife procedure applied to historical forecast errors [26], a uniform distribution [28], and an autoregressive process [20]. Typically, the same process or distribution is used to generate perturbations in each hour of each day in the study, whereas the hourly load forecast error distributions estimated in this paper vary considerably. Perturbations that are independent across hours may result in trajectories that do not well reect the underlying temporal patterns in the load in particular, the ramps that pose the biggest challenge for operators may not be depicted accurately. Time series methods have been used to directly generate multiple trajectories for load [5]. Because of weekly and daily patterns in the load as well as its persistence, the tted processes contain autoregressive terms with multiple lags. The propagation of uncertainty through the lagged terms exposes a drawback of this approach when forecasting more than a few hours in advance. Scenarios for hourly load on day d generated in the early afternoon of day d 1 depend on loads in the evening of day d 1. Even if the forecasts a few

hours ahead are highly reliable, the process of stepping through the forecast hours on day d accumulates increasing numbers of independent random shocks, so that the variance of the generated loads for the later hours of day d grows unrealistically large. See, for example, Fig. 3 in [5], where the dispersion about the mean in the generated set of scenarios appears to increase geometrically with the length of time ahead. In reality, while weather forecasts diminish somewhat in precision as they extend over a day, the ability to forecast load does not drop off as drastically as the ensemble of time series trajectories would suggest.

Because statistical scenario generation approaches require large sample sizes to obtain reliable solutions, scalability of solution procedures may be compromised. To mitigate this difculty, a general method of generating scenarios based on approximating rather than sampling from error distributions was proposed and motivated by load scenario generation in [17]. Our contributions in this paper are (1) the development of methods to incorporate critical data such as seasons, weather patterns, day-of-the-week, and load zones; (2) description of the details necessary to actually generate scenarios; (3) in combination with the companion paper, optimization models that use them and tests on the efcacy of those scenarios; and (4) the scenarios themselves as a resource for the research community. Incorporating numerical weather prediction models in stochastic optimization methods [22,29] may be a promising alternative to generating scenarios from weather forecasts in advance of scheduling. However, becauseto the best of our knowledgesystem operators do not yet have access to such models but rather rely on purchased weather forecasts, we believe our approach (starting with a single weather forecast for the day ahead) is more relevant to current and near-term practice. Another method for short-term load forecasting is to identify similar days within a historical database, where similarity is based on weather, day of the week, and time of year. For example, ISO-NE identies up to ve similar days drawn from the same season with the same day-type according to similarity of their actual temperatures to the forecast temperature of the given day as well as similarity of forecast loads in the last hour of the previous day [8]. Our method has some commonality with this approach, in that we create segments of days that are similar in

123

314 Y. Feng et al.

some sense. Then, within each segment we employ a functional regression method to approximate the probability distribution of load in each hour of the day ahead.

2.2 Methodology overview

We use a multi-step procedure to control for season and type of day, and then approximate the relationships between weather forecast data and the distribution of hourly load sequences within segments of similar days. Starting from a historical database of day-ahead hourly weather forecasts and corresponding actual loads, our load forecasting methodology proceeds as follows:

1. Identify date ranges, or seasons, in which the relationship between weather and loaddisregarding day-of-week effectsis likely to be similar. This ad hoc characterization qualitatively accounts for diurnal light patterns, heating vs. air conditioning impacts, and sociological factors such as whether school is in session.

2. Within each date range, transform the data to account for day-of-the-week and zonal differences within the system. Then, segment the data into bands based on forecast temperature. Data segmentation can in principle proceed using multiple forecast quantities (e.g., temperature and humidity). However, our experiments (summarized very briey in Table 4) indicate that these additional factors do not signicantly improve load forecast accuracy.

3. Within each segment, approximate the relationship between weather and load via a regression function. Additionally, approximate the distribution of residuals associated with the resulting regression function.

Following completion of the segmentation and approximation steps, the procedure for generating load scenarios for a given day d (with associated weather forecast generated on day d 1) is given as follows:1. Identify the date range DR to which day d belongs. Within DR, identify the temperature segment to which the weather forecast for day d belongs.

2. Apply the regression function associated with the identied segment to the weather forecast for day d.

3. Generate forecasted load scenarios for day d using distributions of the forecast errors.

4. As necessary, perform inverse transformations of the load sequences to match the day of the week and the zone.

The segmentation and approximation steps are fully described in the remainder of this section. Details of the load scenario generation process are described subsequently in Sect. 3.

2.3 Estimating regression functions

The main idea to build the regression function is to consider the weather forecast from day d 1 to day d and with this information build a regression function. For this we use

123

Toward scalable stochastic unit commitment 315

second order epi-splines [27] that minimize the deviations from the observed load at the day d. These functions, which have piecewise constant second derivatives, provide great exibility to t varying temporal patterns. For each day d in a given historical segment it is assumed that the following information is available: the hourly load, ldh, and the weather prediction for day d made on day d 1. We use the temperature for

all months and in summer, we also add dew point.We split the 24 h of day d in NR sub-intervals (hk1, hk] of length = 24/NR.

This process determines the total number of coefcients that need to be estimated,i.e., 2 (NR + 2). In the summer, our regression curves will be built by relying on two

epi-splines of order 2, one associated with temperature, and a second one associated with dew point. Each epi-spline is specied by NR constant values for the second derivative in sub-intervals and two integration constants. The construction implicitly assumes that the curves we are tting are twice differentiable, but not necessarily C2. In addition to the parameter NR there is a maximum curvature parameter, . The impact of the parameters on the accuracy of the t is explored in Sect. 4.2.

Further, we assume that the load can be represented as the sum of two components:

A non-weather component, i.e., a component which does not depend on the weather forecast and which is related to the normal behavior of the consumers under conditions specied by each segment.

A weather component, which depends on the weather forecast.

It is natural to think that the non-weather component depends on the segment considered: independently of the weather, people in winter use a different amount of energy than in the summer. Consequently, for each segment we estimate the baseline load, which is the average load for each hour in the segment.

Algorithmic details of regression using epi-splines are provided in [17]. That paper provides an example illustrating how to t a regression function from the weather for day d to R24, using data from some set of days D as

r(d; D).

In the present paper we include details of transformations needed to make the method work effectively in our application. For these descriptions, we use the notations d and do when respectively using forecast and observed weather to estimate transformation functions. A more formal statement would include the regression control parameters, but we omit these in the interest of clarity.

2.4 Data transformation

Within a date range, we use transformations to combine data from disparate day types and zones, prior to constructing regression and error distribution models. Inverse transformations are performed to create load forecasts for particular day types and zones. Without such transformations, there is typically insufcient historical data to yield accurate forecasts. For example, dividing a year into 4 seasonal date ranges yields approximately 12 samples for each week day type. Such small sample sizes can

123

316 Y. Feng et al.

lead to model over-tting and associated out-of-sample prediction inaccuracies. This problem is made more acute if clustering or other classication schemes are imposed on the data within a date range.

Suppose we are given observed load proles ld = (ld1, . . . , ld24)

R24 for a range

of dates d D. A portion of the the load is dependent upon weather factors, but load

proles also depend on the type of the day, e.g., load patterns differ between weekends and weekdays. Based on discernible differences in the average daily load patterns [7], we consider six day types: one for each weekday and one representing weekend days and holidays. We denote the set of day types as T . The set of all dates belonging to a

day type T T is denoted by DT . Clearly, each date maps to a unique day type, such

that [uniontext]T T DT = D and DT DT = if T = T .

While a regression model could be developed for each day type, this would decrease the amount of data available for tting signicantly. Instead, we compute a transformation for each day type to a standard reference day, which we somewhat arbitrarily select as Wednesday. Midweek days have the highest average load and have load proles that are very similar to each other, fairly similar to Monday, and similar to the start of Friday.

The transformation is easily inverted, so that observed loads can be transformed to Wednesday and forecast loads can be transformed back to the original day type. In our analysis, we consider linear transformations to allow for simple forward and inverse computation. We use observed weather to nd a linear transformation from each day type to Wednesday, but the transformation itself is not based on observed weather so load forecasts and scenario generation can be based entirely on data that are available in advance.

For each d DT , let wd = (wd1, . . . , wd24) denote the expected loads for our

reference day type (Wednesday) corresponding to the native day type loads ld =

(ld1, . . . , ld24). We assume that for T T , DT = . For each d DT of day type T ,

wd is computed as a regression ro(do; Dwed), where do denotes the observed weather

for day d.Our goal is then to nd a 24 24 matrix AT such that

AT ld wd d DT . (1)

We formally characterize Eq. 1 as an optimization problem in which the coefcients of AT appear as variables:

min

AT [summationdisplay]

dDT

||AT ld wd||. (2)

Note that this (very small) optimization problem might be linear or non-linear, depending on the choice of norm. In our experiments, we used the l2 norm and employed the open-source non-linear solver Ipopt to solve the resulting optimization problems (see https://projects.coin-or.org/Ipopt

Web End =https://projects.coin-or.org/Ipopt ).

Depending on the relative size of the available dataset {ld|d DT } and the num

ber of coefcients in AT , it may be necessary to introduce additional constraints and regularization terms to formulation (2). Another potential issue is that the resulting

123

Toward scalable stochastic unit commitment 317

AT may be ill conditioned, causing difculties in the calculation A1T. Finally, non-singularity needs to be enforced in (2), which can be easily achieved by requiring that all coefcients above and below the diagonal of AT equal zero; we use this simple approach in the experiments reported below. The linear transformation preserves pairwise correlations among hourly loads, which appear to be fairly consistent across day types [7]. In an application where these correlations differed substantially, the use of a nonlinear transformation could be explored instead.

In practice, a balancing area is typically divided into zones, for which loads are forecasted and reported. In order to increase the data available for regression and error distribution estimation, we additionally combine data from disparate zones in a fashion analogous to that described above for converting data associated with different day types to a reference day type. Data are converted to correspond to data from a reference zone in the same way that Wednesday is used as a reference day.

2.5 Segmentation

For each date range, we partition the weather data for the composite dates into distinct segments. The idea is to limit regression and error distribution inaccuracy by only considering data with similar response characteristics. We segment based on the forecast temperature, as inclusion of additional weather variables (e.g., via k-means clustering) failed to improve prediction accuracy in our experiments. For each date range, we form a temperature distribution over which we introduce cutting points that dene the segments.

Let td be a scalar representation of the hourly forecast temperatures for day d D,

where D is the set of days in the date range. In our experiments, we dene td = td12,

although alternative metrics such as average hourly temperature can be substituted. We require scalar representations of hourly temperature vectors to prevent data for a given day to be mapped to multiple error categories (see Sect. 3).

We estimate the probability density function ft() of the temperature scalar td by

tting an exponential epi-spline [18,19]. We denote the corresponding cumulative distribution function by Ft(). To obtain NS segments of equal probability, we introduce

the break points {b1, . . . , bNS+1} = {0, 1/NS, . . . , (NS 1)/NS, 1} for

Fth() and

then calculate the limit temperatures for each segment Si as

(ti, ti) = (

F1t(bi),

F1t(bi+1)), i = 1, . . . , NS.

Finally, considering the limit temperatures of each segment Si, i = 1, . . . , NS, we

group the days in the date range according to the rule:

d Si td [ti, ti).

3 Scenario generation

To capture the notion that both the mean load response and associated error distributions vary during the day, we split the day into parts and then categorize each portion

123

318 Y. Feng et al.

Fig. 1 Illustrative load regression model and error distributions with hour partition H = {1, 12, 24}

according to the relative error. A regression model is constructed for each of the resulting day parts and associated error categories. Scenarios are then constructed by constructing paths whose trajectories are determined by selecting a specic error category for each day part. This process is a specic instantiation of the general scenario generation methodology detailed in [17].

Let H be the set of hours that dene a partition of the hours in a day, specied as follows:

H = {Hi}|H|i=1, H1 = 1, H|H| = 24, Hi < Hi+1.

The elements Hi represent the partition end-points, e.g., the i-th part of the day is given by the set of hours {Hi, . . . , Hi+1}. For each partition boundary Hi, we compute the

observed regression error di for each day d D as:

di = ldHi rHi (d; D).

We estimate the distribution of these errors by tting an exponential epi-spline [18,19]. This process is graphically illustrated in Fig. 1.

For each partition hour i, we denote the corresponding error probability density function by f i (). Categories within f i () can then be dened through identication of

break points of the associated cumulative distribution function F i (), as was performed

for temperature segmentation. Specically, to generate NK equally sized categories, we select the break points {k1, . . . , kNK+1} = {0, 1/NK , 2/NK , . . . , (NK 1)/NK , 1}

in order to obtain equally weighted categories for each partition of the day. The resulting categories Cki are then dened as Cki = [

F1 i(ki), F1 i(ki+1)), where

123

Toward scalable stochastic unit commitment 319

i {1, . . . , |H| 1} denotes the day partition and k {1, . . . , NK } denotes the

category.Given an hour partition H and associated error categories Cki, we sub-segment the

days d D according to the observed regression model error at the corresponding

partition i. Specically, let Dki denote the set of days in the segment for hour i and

error category k, dened as follows:

d Dki di Cki.

For each sub-segment Dki, we t a regression model r(d; Dki), from which a vector

of predicted hourly loads rd,k is extracted for each day d Dki. Recall that each sub

segment has an associated regression function. Consequently, for boundaries in the middle of the day, there can be two points (i.e., two predicted loads, one from the left day part and one from the right, which might be a different error category.) We avoid discontinuities deriving a single regression curve per category, rd,k, by combining the

predicted points at the limit hours. In particular, for h [Hi, Hi+1],

rd,kh =

[parenleftbigg] Hi+1 h Hi+1 Hi

This process is illustrated in Fig. 2.

Given regression models for each category Cki, the observed load forecast errors are

computed as

kh,d = ldh rd,kh.

We obtain a corresponding probability density function fk

i () by tting an exponential

epi-spline. An illustrative example of this step is shown in Fig. 3, where NK = 2.

The error densities fk

i () serve as the primary input to the scenario generation

process. The rst step in scenario generation involves the identication of a set of distribution cut points C = {cz}|C|z=1 given as parameters, subject to c1 = 0.0 and

c|C| = 1.0. For each partition i and category k, we then calculate the conditional

expected value of the error in each interval dened by a pair of adjacent cutting points:

E[bracketleftBig]F1k

i (cz), F1ki(cz+1)

= k,zi

where k,zi denotes the expected error in category k at cutting point z for hour i. The conditional probability associated with this error point is the width of the interval and the unconditional probability is obtained by multiplying by the width of the error category (which is 1/NK ). The number of cutting points can vary per hour, as shown

[parenrightbigg] rd,kh,i +

[parenrightbigg] rd,kh,i+1.

[parenleftbigg] h Hi

Hi+1 Hi

[bracketrightBig]

[integraldisplay] cz+1

x f

ki (x)dx

i (x)dx

123

320 Y. Feng et al.

Fig. 2 Regression of the error category construction

in Fig. 4. In this illustrative example, C = {0.0, 1.0} for i = 1, but C = {0.0, 0.5, 1.0}

for i = 12, 24. The cutting points also need not be symmetric, although they are in

this example.Given regression models (rk,dh,i) for each hour partition boundary i and category k,

we compute loads at the partition boundaries via:

ld,k,zHi = rd,kHi + k,zi.

For each hour Hi and each category Cki, this step yields |Ci| forecast load samples. The

nal step in our scenario generation process is to connect these samples in order to construct a set of paths that approximates the stochastic process representing load for the full day. This is simply done by calculating the scenario loads at time h [Hi, Hi+1)

by assuming that the deviation from the forecast varies between the deviation at hour Hi and hour Hi+1. This process is illustrated in Fig. 5. Under this methodology, the

number of paths (i.e., scenarios) generated depends on the parameters. For example, if the number of cutting points is the same value, |C|, for every hour partition boundary,

then the number of scenarios would be:

123

Toward scalable stochastic unit commitment 321

Fig. 3 Illustrative error category regressions and associated error distributions with hour partition H = {1, 12, 24} and NK = 2 categories

Fig. 4 Distribution cutting points for scenario generation when H = {1, 12, 24} and NK = 2. For i = 1, C = {0.0, 1.0}. For i {12, 24}, C = {0.0, 0.5, 1.0}

123

322 Y. Feng et al.

Fig. 5 Illustrative scenario paths corresponding to the cutting points shown in Fig. 4

NK (|C| 1)|H|1.

The probability associated with each scenario must be computed as the product of the unconditional probabilities of the points use to construct it. Note that the generation process is deterministic, given a xed set of historical input data.

4 Experimental results

Although substantially better MAPEs can be achieved using standard leave-one-out validation, we have generated our forecasts and scenarios by simulating a rolling horizon as would be seen by a real-world system operator. In our experiments, we consider data for ISO-NE. We begin by tting the models using only data from 2009 and 2010, and consider operations during 2011. As we simulate the progression through the year, data from 2011 are added to the t process as it becomes historical. Specically, we execute the complete methodology outline in Sect. 2.2 for each simulated day. The entire procedure takes minutes of wall clock time to complete, and is therefore feasible in practice. Because our interest is in demonstrating load scenario generation methods, and not actually providing ISO-NE with load forecasts, we begin our simulation on

123

Toward scalable stochastic unit commitment 323

Table 1 Load zones for

ISO-NE, with weather stations and corresponding weights

January 2, 2011 and end on November 20, 2011. Excellent methods exist for dealing with the holiday season in the US [9], but their use is beyond the present scope. Further, we ignore August 2830 of 2011, due to a hurricane event in the region. We partition the days of the year into the following date ranges: Winter (January 2March31), Spring (April 1May 14), Summer (May 15September 14), and Fall (September 15November 20).

As reported in Table 1, ISO-NE is divided into 8 load zones. Weather data for each zone are taken from one or two weather stations. For zones with two stations, an aggregate weather forecast is computed by weighting the composite station data appropriately. Historical load data were obtained by ISO-NE through their public web site (http://www.iso-ne.com/markets/hstdata/znl_info/hourly/index.html

Web End =http://www.iso-ne.com/markets/hstdata/znl_info/hourly/index. http://www.iso-ne.com/markets/hstdata/znl_info/hourly/index.html

Web End =html ). Hourly day-ahead temperature and dew point forecasts were provided directly by ISO-NE.

As indicated in Sect. 2.4, we aggregate data from disparate zones to make more data available to the t process. Zone aggregation proceeds as follows, by employing two reference zones: CT and NEMASSBOST. The RI, SEMASS, and WCMASS zones are aggregated with the CT zone, while the ME, NH, and VT zones are aggregated with the NEMASSBOST zone. This aggregation corresponds to a partition of ISONE in approximately northern and southern regions, which in turn share similar load characteristics. Within each partition, we select the zone with the greatest load as the reference zone, minimizing the total load transformed.

4.1 Forecast MAPEs

We quantify load forecast accuracy as the Mean Absolute Percent Error (MAPE), denoted by MAPE(NR, ), as

[summationtext]zZ [summationtext]dD [summationtext]hH |(lz,dh

E(lz,dh)(NR, ))|/lz,dh |Z| |D| |H|

Load zone Weather stations Weights

ME PWM 1.0

NH CON 1.0

VT BTV 1.0

CT BDL, BDR 0.13, 0.87

RI PVD 1.0

SEMASS PVD 1.0

WCMASS BDL, ORH 0.5, 0.5

NEMASSBOST BOS 1.0

100 (3)

where Z, D, and H respectively denote the sets of load zones, dates, and hours under

consideration. NR and denote regression t parameters. The aggregated MAPEs obtained for each date range in 2011 are reported in Table 2, considering NR = 24,

123

324 Y. Feng et al.

Table 2 Aggregated ISO-NE

MAPEs for 2011 Season Segments (NS)

1 3 5 7

Fall 5.45 4.66 4.20 3.99

Spring 3.10 2.88 2.67 2.73

Summer 10.25 4.82 4.14 4.19

Winter 5.25 3.32 3.29 3.47

= 500, and a variable number of temperature segments NS. The corresponding dis-

aggregated (by zone) MAPEs are reported in Table 3. We observe that segmentation of the data by temperature does improve load forecast accuracy, although the benets either stagnate or decrease once NS 7. Further, the data in Table 2 demonstrate

load forecast accuracies that are consistent with those obtained by ISO-NE in practice, e.g,. see http://www.iso-ne.com/support/training/courses/wem101/10_forecast_scheduling_callan

Web End =http://www.iso-ne.com/support/training/courses/wem101/10_forecast_ http://www.iso-ne.com/support/training/courses/wem101/10_forecast_scheduling_callan

Web End =scheduling_callan . However, our method additionally provides estimates of forecast variability, represented through a collection of load scenarios.

We now briey consider inclusion of dew point in addition to temperature in the segmentation process. Intuitively, dew point can inuence load by impacting cooling requirementsparticularly in the summer. To address this case, we segment forecast data by forming a scalar weather quantity based on a linear combination of temperature and dew point, using observations at h = 12. The resulting MAPEs, aggregated across

all zones, are shown in Table 4, and indicate that inclusion of dew point quantities can marginally improve load forecast accuracy in specic contexts (e.g., in summer).

Finally, we provide exemplars of forecast load scenarios for the CT zone in ISONE, for early summer days in 2011. We have selected low and high variance load 50-scenario cases, respectively shown in Figs. 6 and 7. We show only the rst 24 h of the time series, which repeat (for the purpose of stochastic unit commitment) for an additional 24 h. With the exception of low-load hours 7 and 8, the scenarios shown in Fig. 6 envelop the actual load. Further, during peak load for this day, the actual load closely mirrors the regression load. The situation is reversed in Fig. 7, where the regression load mirrors the actual load in early hours of the day, but diverges after hour 18. These two particular cases serve as an integral component of test cases for assessing our stochastic unit commitment solver, described in the companion paper. All load scenario sets generated for ISO-NE for 2011, either at the zonal or aggregate level, are publicly available at https://prescient.sandia.gov

Web End =https://prescient.sandia.gov for unrestricted use. We note that specic causes of low versus high variance days are very difcult to identify. Rather, we simply observe that the differences are due to differences in the historical forecast errors associated with days of similar types and under similar forecast weather conditions.

4.2 Parameter sensitivity for model tting via epi-splines

To construct our regression models via epi-splines, it is necessary to specify values for the following key parameters:

123

Toward scalable stochastic unit commitment 325

Table 3 Disaggregated zonal

ISO-NE MAPEs for 2011 Season Zone Segments (NS)

1 3 5 7

Fall NH 4.44 4.22 4.41 4.14

VT 3.52 2.91 2.95 3.04

ME 4.13 4.14 4.02 4.15

CT 8.31 7.12 6.77 6.51

RI 5.93 5.30 4.49 4.00

SEMASS 5.75 4.99 4.40 3.84

WCMASS 7.20 6.60 6.34 6.66

NEMASSBOST 4.82 4.63 4.19 4.07

Spring NH 3.14 3.46 3.69 3.93

VT 3.23 3.18 3.04 3.22

ME 4.25 3.90 4.07 4.00

CT 3.46 3.64 3.42 3.89

RI 3.24 2.97 3.27 3.32

SEMASS 3.05 2.93 3.17 3.24

WCMASS 3.74 3.50 3.80 3.98

NEMASSBOST 3.24 3.41 3.35 3.48

Summer NH 9.29 5.65 4.95 5.05

VT 5.64 3.75 3.41 3.40

ME 7.54 4.65 4.83 4.54

CT 11.22 6.42 5.98 5.84

RI 12.90 7.15 5.69 5.86

SEMASS 12.72 6.76 5.60 5.74

WCMASS 9.47 5.51 4.57 4.63

NEMASSBOST 11.34 6.17 5.78 5.32

Winter NH 4.99 3.88 3.78 3.82

VT 4.28 3.77 4.02 4.17

ME 4.10 3.87 3.65 4.05

CT 6.17 4.35 4.21 4.32

RI 5.23 3.81 3.56 3.84

SEMASS 5.34 4.05 3.76 4.11

WCMASS 5.46 3.69 3.84 3.95

NEMASSBOST 5.42 3.54 3.62 3.67

NR: The number of intervals into which the hours of a day are sub-divided. : The maximum curvature of the regression model.

While it is necessary to specify a particular norm (e.g., L1 or L2) when solving the embedded optimization problems, in practice the choice of norm has almost no impact on the resulting load forecast MAPEs. In contrast, accuracy is more sensitive to the choice of NR and , as we now demonstrate.

123

326 Y. Feng et al.

Table 4 Aggregated ISO-NE

MAPE for Summer 2011, using segmentation based on temperature both with or without dew point

First, we consider the impact of NR on load forecast accuracy. We compute the MAPEs for NR {6, 12, 18, 24, 32, 48}, xing NS = 5 and = 500. The results

indicate that the MAPE changes from about 4.7 to 4.5% over the range, with values of 24 and above yielding nearly equal MAPEs.

Next, we consider the impact of on load forecast accuracy. Fixing NR = 24 and

NS = 5, we compute MAPEs for {20, 40, 60, 80, 100, 150, 200, 500}. Extremely

small values of do adversely impact the MAPEs by severely restricting the curvature. However, over the range of we tested, there is essentially no sensitivity.

5 Conclusion

In this paper, we introduced novel methods for obtaining distributions of forecast load in each hour of day d, based on weather forecasts available on day d 1. Our

goal is to estimate trends and error distributions from these data, in order to generate probabilistic scenarios for day-ahead load. We make use of transformations of day types and load zones to provide data for general, functional regression based on epi-

Condition Segments (NS)

1 3 5 7

Without DewPt 10.25 4.82 4.14 4.19

With DewPt 10.55 4.73 4.06 4.06

Fig. 6 Scenarios for low variance load forecast day in 2011 for zone CT

123

Toward scalable stochastic unit commitment 327

splines that is conditional on the error category. We also t the error distributions using epi-splines.

Our experiments indicate that our models:

Can be generated using data readily available to system operators. Can be produced in minutes of run time for multiple years of input data. Produce errors that in expectation are competitive in the aggregate with load forecasting procedures used by industry. Obtain MAPE values for forecast load that are similar to those found in hind-casting studies, which eliminate weather forecast uncertainty by focusing strictly on actual or observed weather quantities. Explicitly model likely variability of realized load.

While we do not examine multi-stage stochastic unit commitment models in the companion paper [4], we note that our scenarios are tree-structured and consequently can be effectively used in that context.

We continue to rene our methods, specically by focusing on reducing prediction errors associated with peak load periods and adapting the approach to times of the year where temperature is not as strong a predictor of the load. Our approach could also be extended to model price-responsive demand when such programs are widely implemented and data become available to estimate demand-price relationships.

The sets of scenarios generated by the methods we describe serve as input to rigorously test the scalability of stochastic unit commitment solvers, as described in the companion paper. Scenarios generated in a rigorous way from realistic input data are important for ongoing research in stochastic and robust unit commitment, and for sim-

Fig. 7 Scenarios for high variance load forecast day in 2011 for zone CT

123

328 Y. Feng et al.

ulation and related production cost models. The scenarios generated by our method using data from ISO-NE for 2011 are available as an online supplement to the paper.

Acknowledgments Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energys National Nuclear Security Administration under Contract DE-AC04-94-AL85000. This work was funded by the Department of Energys Advanced Research Projects Agency - Energy, under the Green Energy Network Integration (GENI) project portfolio. The authors would like to thank Dr. Eugene Litvinov and his research group at ISO-NE (in particular, Bill Callan) for their assistance with ISO-NE system and data sources.

1. Black, J.: Load hindcasting: a retrospective regional load prediction method using reanalysis weather data. Ph.D. thesis, Department of Mechanical and Industrial Engineering, University of Massachusetts Amherst (2011)

2. Bouffard, F., Galiana, F.D., Conejo, A.J.: Market-clearing with stochastic securityPart 1: formulation. IEEE Trans. Power Syst. 20(4), 18181826 (2005)

3. Carpentier, P., Gohen, G., Culioli, J.C., Renaud, A.: stochastic optimization of unit commitment: a new decomposition framework. IEEE Trans. Power Syst. 11, 10671073 (1996)

4. Cheung, K., Gade, D., Monroy, C.S., Ryan, S.M., Watson, J.P., Wets, R.J.B., Woodruff, D.L.: Scalable stochastic unit commitment, part 2: solver performance. Energy Syst (to appear)

5. Dupaov, J., Grwe-Kuska, N., Rmisch, W.: Scenario reduction in stochastic programming. Math. Program. 95(3), 493511 (2003). doi:http://dx.doi.org/10.1007/s10107-002-0331-0

Web End =10.1007/s10107-002-0331-0

6. Feinberg, E.A., Genethliou, D.: Load forecasting. In: Chow, J.H., Wu, F.F., Momoh, J. (eds.) Applied Mathematics for Restructured Electric Power Systems, pp. 269285. Springer (2005)

7. Feng, Y., Gade, D., Ryan, S.M., Watson, J.P., Wets, R.J.B., Woodruff, D.L.: A new approximation method for generating day-ahead load scenarios. In: Proceedings of the 2013 IEEE Power and Energy Society General Meeting (2013)

8. Gould, S.: Create demand forecast. Tech. Rep. SOP-OUTSCH.0040.0010, ISO New England, Inc (2011)

9. Hong, T.: Electric load forecasting with holiday effect. https://sites.google.com/site/hongtao01/courses/holiday

Web End =https://sites.google.com/site/hongtao01/ https://sites.google.com/site/hongtao01/courses/holiday

Web End =courses/holiday (2012). Presented at AEIC Load Research Workshop 2012, Orlando, FL

10. Hong, T., Gui, M., Baran, M., Willis, H.: Modeling and forecasting hourly electric load by multiple linear regression with interactions. In: Proceedings of the 2010 IEEE Power and Energy Society Meeting (2010)

11. Hong, T., Pinson, P., Fan, S.: Global energy forecasting competition 2012. Int. J. Forecast. (2013). doi:http://dx.doi.org/10.1016/j.ijforecast.2013.07.001

Web End =10.1016/j.ijforecast.2013.07.001

12. Huang, Y., Zheng, Q.P., Wang, J.: Two-stage stochastic unit commitment model including non-generation resources with conditional value-at-risk constraints. Electr. Power Syst. Res. 116, 427438 (2014)

13. Klboe, G., Eriksrud, A.L., Fleten, S.E.: Benchmarking time series based forecasting models for electricity balancing market prices. Energy Syst. (2013). doi:http://dx.doi.org/10.1007/s12667-013-0103-3

Web End =10.1007/s12667-013-0103-3

14. Liu, J., Chen, R., Liu, L., Harris, J.: A semi-parametric time series approach in modeling hourly electricity loads. J. Forecast. 25, 537559 (2006)

15. Morales, J.M., Conejo, A.J., Liu, K., Zhong, J.: Pricing electricity in pools with wind producers. IEEE Trans. Power Syst. 27(3), 13661376 (2012)

16. Papavasiliou, A., Oren, S.: A stochastic unit commitment model for integrating renewable supply and demand response. In: Proceedings of the 2012 IEEE Power and Energy Society Meeting (2012)17. Rios, I., Wets, R.J.B., Woodruff, D.L.: Multi-period forecasting and scenario generation with limited data. Comput. Manag. Sci. (to appear)

18. Royset, J., Wets, R.: Epi-splines and exponential epi-splines: pliable approximation tools. Tech. rep., University of California Davis (2012)

19. Royset, J., Wets, R.: Fusion of hard and soft information in non-parametric density estimation. Tech. rep., University of California Davis (2012)

References

123

Toward scalable stochastic unit commitment 329

20. Ruiz, P., Philbrick, C., Sauer, P.: Modeling approaches for computational cost reduction in stochastic unit commitment formulations. IEEE Trans. Power Syst. 25(1), 588589 (2010)

21. Ruiz, P., Philbrick, R., Zack, E., Cheung, K., Sauer, P.: Uncertainty management in the unit commitment problem. IEEE Trans. Power Syst. 24(2), 642651 (2009)

22. Salazar, J.M., Diwekar, U., Constantinescu, E., Zavala, V.M.: Stochastic optimization approach to water management in cooling-constrained power plants. Appl. Energy 112(0), 1222 (2013). doi:http://dx.doi.org/10.1016/j.apenergy.2013.05.077

Web End =10. http://dx.doi.org/10.1016/j.apenergy.2013.05.077

Web End =1016/j.apenergy.2013.05.077

23. Scharff, R., Egerer, J., Sder, L.: A description of the operative decision-making process of a power generating company on the nordic electricity market. Energy Syst. 5(2), 349369 (2014). doi:http://dx.doi.org/10.1007/s12667-013-0104-2

Web End =10.1007/ http://dx.doi.org/10.1007/s12667-013-0104-2

Web End =s12667-013-0104-2

24. Siface, D., Vespucci, M.T., Gelmini, A.: Solution of the mixed integer large scale unit commitment problem by means of a continuous stochastic linear programming model. Energy Syst. 5, 269284 (2014). doi:http://dx.doi.org/10.1007/s12667-013-0107-z

Web End =10.1007/s12667-013-0107-z

25. Takriti, S., Birge, J., Long, E.: A stochastic model for the unit commitment problem. IEEE Trans. Power Syst. 11(3), 14971508 (1996)

26. Takriti, S., Krasenbrink, B., Wu, L.S.Y.: Incorporating fuel constraints and electricity spot prices into the stochastic unit commitment problem. Oper. Res. 48(2), 268280 (2000)

27. Wets, R., Bianchi, S.: Term and volatility structures. In: Zenios, S., Ziemba, W. (eds.) Handbook of Asset and Liability Management, pp. 2668. Elsevier (2006)

28. Wu, L., Shahidehpour, M., Li, T.: Stochastic security-constrained unit commitment. IEEE Trans. Power Syst. 22(2), 800811 (2007)

29. Zavala, V.M., Constantinescu, E.M., Krause, T., Anitescu, M.: On-line economic optimization of energy systems using weather forecast information. J. Process Control 19(10), 17251736 (2009). doi:http://dx.doi.org/10.1016/j.jprocont.2009.07.004

Web End =10.1016/j.jprocont.2009.07.004

123

Word count: 7942

Show less

Springer-Verlag Berlin Heidelberg 2015

Toward scalable stochastic unit commitment. Part 1: load scenario generation

Content area

Abstract

Full text