Optimal Hour-Ahead Bidding in the Real-Time Electricity Market with Battery Storage using Approximate Dynamic Programming

Introduction

Bidding into the electricity market can be a complicated process, mainly due to the requirement of balancing supply and demand at each point in the grid. To solve this issue, the Independent System Operators (ISOs) and the Regional Transmission Organizations (RTOs) generally use multi–settlement markets: several tiers of markets covering planning horizons that range from day–ahead to real–time. The idea is that the markets further away from the operating time settle the majority of the generation needed to handle the predicted load, while the markets closer to the operating time correct for the small, yet unpredictable deviations that may be caused by issues like weather, transmission problems, and generation outages (see, for example, Shahidehpour et al. (2002), Eydeland and Wolyniec (2003), Harris (2011), for more details). Settlements in these real–time markets are based on a set of intra–hour prices, typically computed at 5, 10, or 15 minute intervals, depending on the specific market in question. A settlement refers to the financial transaction after a generator clears the market, which refers to being selected to either buy or sell energy from the market. If a generator does not clear the market, it remains idle and no settlement occurs. We refer to this situation as being out of the market.

Many ISO’s and RTO’s, such as the Pennsylvania–New Jersey–Maryland Interconnection (PJM), deal with the balancing market primarily through the day–ahead market. PJM’s balancing market clears every 5 minutes (considered to be near “real–time”), but the bids are all placed the previous day. See Eydeland and Wolyniec (2003) and the PJM Energy and Ancillary Services Market Operations Manual for more information. In certain markets, however, it is not only possible to settle in real–time, but market participants can also submit bids each hour, for an hour in the future. Thus, a bid (consisting of buy and sell prices) can be made at 1pm that will govern the battery between 2pm and 3pm. The process of both bidding and settling in real–time is a characteristic of the New York Independent System Operator (NYISO) real–time market and is the motivating example for this paper. Other prominent examples of markets that include a real–time bidding aspect include California ISO (CAISO) and Midcontinent ISO (MISO). In particular, our goal is to pair battery storage with hour–ahead bidding in the real–time market for revenue maximization, a strategy sometimes referred to as energy arbitrage.

It is unlikely that profits from battery/energy arbitrage alone can be sustainable for a company; however, if performed optimally, it can be an important part of a range of profit generating activities (one such example is the frequency regulation market). See Walawalkar et al. (2007) for an economic analysis of using a storage device for both energy arbitrage (using a simple “charge–off–peak and discharge–on–peak” policy) and frequency regulation in the New York area. The analysis shows that in New York City (but not the surrounding areas), there is a “high probability of positive NPV [net present value] for both energy arbitrage and regulation,” but even so, there is still significant risk in not being able to recover the initial capital cost. However, the potential for more efficient and cost–effective technology combined with better control policies can make energy arbitrage feasible in the near future. Other studies on the topic of the value of storage include Sioshansi et al. (2009), Sioshansi et al. (2011), and Byrne and Silva-Monroy (2012).

In our problem, we assume that the goal is to optimally control a 1 MW battery; in practice, a company may operate a fleet of such batteries. Market rules state that we must bid in integer increments, meaning the possible actions at each settlement are to charge, discharge (both at a rate of 1 MW), or do nothing. Hence, our precise problem is to optimize the placement of two hour–ahead bids, a “positive” bid (for a quantity of $+1$ MW) and a “negative” bid (for a quantity of $-1$ MW) that correspond to selling (generation) and buying (negative generation), respectively, over a period of time such that purchased energy can be stored in the finite capacity battery. The goal is to maximize expected revenue. Further, given that our model is tailored to battery storage (inherently small capacity), it is reasonable to assume no price impact (i.e., our bids do not affect the spot prices of electricity). In the real–time market, bidding for the operating hour closes an hour in advance and the hour–ahead bid is fixed for the entire operating hour.

This paper makes the following contributions. We describe, in detail, a mathematical model of the bidding process in the real–time electricity market and formulate the sequential decision problem as a Markov Decision Process (MDP). Along the way, we show the structural properties of the problem (monotonicity of the contribution and value functions) that we utilize in our solution technique. Next, we describe and benchmark a convergent approximate dynamic programming algorithm called Monotone–ADP (M–ADP) (Jiang and Powell (2015)) that can be used to obtain an approximate, but near–optimal bidding policy. We also present a new version of Monotone–ADP utilizing post–decision states that allows us to train bidding policies without any model or knowledge of the distribution of real–time prices (which we call a distribution–free method), allowing our solution technique to be easily adopted in practice. Finally, we present a case study detailing the results of an ADP policy trained using only historical real–time price data from the NYISO. In the case study, we also compare the ADP policy to other rule–based policies, two of which are from the energy arbitrage literature and one from our industry contacts. All proofs can be found in the Appendix.

Literature Review

With renewable energy sources like wind and solar becoming more established, the problem of energy storage is also becoming increasingly important. In this section, we first review studies dedicated solely to storage and then move on to those that consider the bidding aspect. Lastly, we discuss algorithmic techniques similar to our proposed method (Monotone–ADP).

Coupling wind energy with storage has been well–studied in a variety of ways. The paper by Kim and Powell (2011) poses a wind energy commitment problem given storage and then analytically determines the optimal policy for the infinite horizon case. Sioshansi (2011) uses ideas from economics and game theory (i.e., the Stackelberg Model) to make several conclusions, including the finding that the value of storage increases with market–competitiveness. In addition, Greenblatt et al. (2007) finds that for high green house gas (GHG) emissions prices, compressed air energy storage is a better choice as a supplemental generator to wind energy when compared to natural gas turbines. The well–known smoothing effects of energy storage on intermittent renewable sources is studied in the context of wind power output by Paatero and Lund (2005).

Another problem within this realm is the storage of natural gas, which involves optimally controlling injection and withdrawal of gas from a storage facility that is typically underground. Carmona and Ludkovski (2010) uses a technique known as optimal switching to solve a natural gas storage problem; computationally, the value function is approximated using basis functions. In a similar vein, Thompson et al. (2009) formulates a stochastic control problem and numerically solve the resulting integro–differential equation to arrive at the optimal policy. Lai et al. (2010) proposes using an ADP algorithm along with an approximation technique to reduce the number of state space dimensions for natural gas storage valuation.

Other energy storage problems include reservoir management (see Nandalal and Bogardi (2007)) and pairing solar with battery storage (see Barnhart et al. (2013)). It quickly becomes clear that all of these problems are similar; in fact, Secomandi (2010) gives the structure of the optimal policy for trading generic commodities given storage. At its core, energy storage has similarities to an array of classical problems related to operations research, such as resource allocation and inventory control.

There are also many studies that consider the bidding aspect of the electricity markets. One significant difference between many of these studies and our paper is that, rather than placing many bids at once, we consider a sequential, hourly bidding problem. Löhndorf and Minner (2010) considers a day–ahead bidding problem different from ours using an infinite horizon MDP; Conejo et al. (2002) solves a price–taker bidding problem using a deterministic look–ahead policy; Gross and Finlay (2000) formulate a constrained optimization problem for optimal bidding in a competitive power pool; and David (1993) develops both deterministic and stochastic models for bidding under the consideration of other market players. Lastly, Löhndorf et al. (2013) uses approximate dual dynamic programming (ADDP) to solve a day–ahead bidding problem involving hydro storage. Besides the major algorithmic differences from our paper, Löhndorf et al. (2013) also works in a day–ahead setting with individual bids for each hourly subinterval, while we work in an hourly setting with bids that must be simultaneously active for every 5 minute subinterval. Furthermore, in order to have complete information to make the optimal decision and to implement the transition dynamics, the previous bid (placed in the last time interval) is a part of our state variable, which is not the case for Löhndorf et al. (2013). For more details, the literature survey by Wen and David (2000) provides an excellent overview to strategic bidding.

In the case of real–world problems with large state spaces, backward dynamic programming is typically not a viable solution strategy, so we often use approximate dynamic programming (ADP) techniques. In this paper, we consider a variant of the approximate value iteration (AVI) algorithm (see both Bertsekas and Tsitsiklis (1996) and Powell (2011)) that exploits the monotonicity in certain dimensions of the optimal value function (also known as the cost–to–go function) in order to quickly approximate the shape of the value function. The algorithm, called Monotone–ADP, is analyzed in Jiang and Powell (2015) and was used previously as a heuristic in Papadaki and Powell (2003).

Like monotonicity, convexity/concavity also often arise in applications, and similar algorithms to Monotone–ADP that exploit these structural properties have been studied in Godfrey and Powell (2001), Topaloglu and Powell (2003), Powell et al. (2004), and Nascimento and Powell (2009). In general, the above studies on monotonicity and convexity have shown that it is advantageous to use the structural properties of value functions in ADP algorithms.

Mathematical Formulation

Hence, $P_{(t,t+1]}$ does not become fully known until time $t+1$. Next, let our set of bidding decisions be a finite set $\mathcal{B}$ such that

We remark that in the actual bidding process, the buy bid is a negative number and the criteria for clearing the market is that the bid is less than the negative of the spot price. Due to our bids being for only two quantities ($\pm 1$ MW), the above reformulation of the bidding process is cleaner and more intuitive.

where $e_{m}$ is a vector of zeros with a one at the $m$–th row (and thus, picks out the $m$–th component of the price vector $P$). Note that $q$ is not dependent on time, but in the context of our hour–ahead bidding problem, we use it in the form of $q(P_{(t-1,t]},b_{t-2})$, which is deterministic at time $t$. Figure 1 illustrates the intra–hour behavior.

To define the hourly transition function between $R_{t}$ and $R_{t+1}$, we model each of the individual settlements within the hour and then combine them recursively (since from $t$ to $t+1$, we settle $M$ times). Let $q_{s}\in\{-1,0,1\}^{M}$ be a vector of settlement outcomes and suppose $g^{R}_{m}(R_{t},q_{s})$ represents the amount of resource after the $m$–th settlement. Thus, we have

for $1\leq m\leq M$. The intra–hour resource levels are

Finally, let $g^{R}$ be the hourly transition function, which is defined as a composition of the functions $g^{R}_{M}$ and $q$ in the following way:

The need for an hourly transition function from $R_{t}$ directly to $R_{t+1}$ (rather than simply defining the sub–transitions between the intra–hour settlements) is due to the hourly decision epoch of the problem.

We now consider another dimension to our problem by allowing a limit to be imposed on the number of charge–discharge cycles used by the battery, for the sake of increasing the lifetime of the battery. Battery cycle–life (the approximate number of cycles before capacity diminishes to around 80%), a key issue when considering economic feasibility, varies between the different types of battery technologies and the operating conditions, but are typically in the range of 1000 (e.g., lead–acid) to 5000 (e.g., vanadium redox); for an extensive review, see Yang et al. (2011). In our correspondence with industry colleagues, we found that a common (though possibly somewhat conservative) estimate of battery usage is 300 cycles/year, meaning that most devices can last at least 3 to 4 years. However, the model developed in this paper is for hourly decision making and it would be impractical to solve the model for time–horizons of several years. Note that different batteries technologies degrade in different ways, but in general, degradation occurs slowly (nearly linearly with charge–discharge cycles) at first, but after a point, efficiency drops much more rapidly.

Over a short horizon (on the order of days), the effects of battery degradation is negligible, but we propose the following way for one to impose a sort of artificial limit to the number of trades (charge–discharge cycles) performed. Let $L_{t}\in\mathcal{L}=\{0,1,2,\ldots,L_{\text{{max}}}\}$ be decremented on every discharge of the battery (starting with $L_{0}=L_{\text{{max}}}$) and suppose that when selling to the market at a settlement time $t^{\prime}$ in $(t,t+1]$, the revenue is discounted by a factor of $\beta(L_{t^{\prime}})$ where $\beta:\mathcal{L}\rightarrow$ is a nondecreasing function. Depending on the battery technology, preferences of the operator, and the time–horizon of the model, the choice of $\beta$ may vary greatly; the list below offers a few examples:

Constant: $\beta(l)=c\in$ for all $l\in\mathcal{L}$,

Step: $\beta(0)=0$ and $\beta(l)=1$ for $l\in\mathcal{L}\setminus\{0\}$,

Linear: $\beta(l)=l/L_{\text{{max}}}$ for all $l\in\mathcal{L}$,

Power: $\beta(l)=(l/L_{\text{{max}}})^{\frac{1}{n}}$ for some $n>1$ and all $l\in\mathcal{L}$,

where (4) seeks to very roughly mimic the efficiency degradation of a real battery. We assume that the physical characteristics of the battery are summarized through $\beta$ and the dynamics of $L_{t}$, which we now describe.

which indicates the settlements for which a discharge occurred. Like before, we define the transition function from $L_{t}$ to $L_{t+1}$ using a sequence of sub–transitions. Let $d_{s}\in\{0,1\}^{M}$ be a vector of settlement outcomes (in this case, whether a discharge happened or not) and

for $1\leq m\leq M$. The intra–hour values are

and the hourly transition function $g^{L}$ is defined

The undersupply penalty takes values of either 1 (no penalty) or $-K$ (penalty):

This penalization scheme reflects reality: the NYISO penalizes using a price–proportional penalty of $K=1$ (in addition to lost revenue), the reason being to uphold the market balance. When a market participant reneges on a promise to deliver energy to the market, it must pay the penalty of the quantity times the market price to correct the imbalance; this is equivalent to purchasing that energy from another generator at the market price and delivering to the market.

Hence, we can write the following sum (over the settlements) to arrive at the hourly revenue, denoted by the function $C$:

Note that $C$ is not time–dependent. The timeline of events and notation we use is summarized in Figure 2. The top half of Figure 2 shows the contrast between when bids are placed and when bids are active: $b_{t}$ and $b_{t+1}$ are placed at times $t$ and $t+1$ (arrows pointing up), while $b_{t-1}$ is active for the interval $(t,t+1]$ and $b_{t}$ is active for the interval $(t+1,t+2]$. It also shows that the revenue function $C(R_{t},L_{t},P_{(t,t+1]},b_{t-1})$ refers to the interval $(t,t+1]$. The bottom half of Figure 2 shows an example of the bidding outcomes, i.e., the output of $q(P_{(t,t+1]},b_{t-1})$. Finally, we emphasize that $M$ settlements (and thus, transitions) occur between consecutive values of $R_{t}$ and $L_{t}$ due to the discrepancy between the bidding timeline (hourly) and the settlement timeline (every five minutes).

The problem of optimizing revenue over a time horizon is a sequential decision problem that we can formulate as a Markov Decision Process (MDP). First, suppose the set of state variables associated with the price process $P_{t}$ is denoted $P^{S}_{t}\in\mathcal{P}$, where $\mathcal{P}$ is the space of price model state variables. The MDP can be characterized by the following components:

The state variable for the overall problem is $S_{t}=(R_{t},L_{t},b_{t-1}^{-},b_{t-1}^{+},P^{S}_{t})\in\mathcal{S}$ where $\mathcal{S}$ is the state space. The previous bid $b_{t-1}$ is included because it is the bid that becomes valid at time $t$ for the interval $(t,t+1]$ and is necessary for computing the resource transition function.

The decision is the hour–ahead bid $b_{t}=(b_{t}^{-},b_{t}^{+})\in\mathcal{B}$ that is active for the interval $(t+1,t+2]$.

The exogenous information in this problem is the price process $P_{t}$.

The state transition function or system model $S^{M}$ is given by

The contribution function in this model represents the expected value of the revenue in the interval from $t+1$ to $t+2$ using bid $b_{t}$ given the current state $S_{t}$. Define:

Let $T-1$ be the last time for which a bid needs to be placed (hence, the trading horizon lasts until $T+1$ and the last value function we need to define is at $T$) and let $B_{t}^{\pi}:\mathcal{S}\rightarrow\mathcal{B}$ be the decision function for a policy $\pi$ from the class $\Pi$ of all admissible policies. The following is the objective function for maximizing expected revenue:

where $C_{\textnormal{term}}(S_{T})$ represents a terminal contribution that is nondecreasing in $R_{T}$, $L_{T}$, and $b_{T-1}$.

We can express the optimal policy in the form of a stochastic dynamic program using Bellman’s optimality equation (Bellman, 1957). The optimal value function $V^{*}$ is defined for each $t$ and each state $S_{t}$:

Figure 3 illustrates the above notation. Notice that at any decision epoch $t$, both the contribution and value functions are looking one step ahead, i.e., from $t+1$ onwards, in the form of an expectation. Because of this, the revenue from $t$ to $t+1$ become, in a sense, irrelevant. However, the link between the time periods comes from the dependence of $R_{t+1}$ and $L_{t+1}$ on $R_{t}$, $L_{t}$, and $b_{t-1}$ (and of course, the random prices). In other words, at time $t$, our bid has to be placed for $(t+1,t+2]$ with an uncertain amount of resource, $R_{t+1}$ in the battery. It is important to note that it is precisely because $C_{t,t+2}(S_{t},b_{t})$ does not include the revenue made in $(t,t+1]$ that allows us to show the important structural property of monotonicity for $C_{t,t+2}$ in $b_{t-1}$ (see Proposition 3 in the next section).

We now provide some results regarding the structure of the contribution and value functions. The algorithm (Monotone–ADP–Bidding) that we implement to solve for the optimal value function is inspired by the following monotonicity properties.

The contribution function ${C}_{t,t+2}(S_{t},b_{t})$, with $S_{t}=(R_{t},L_{t},b_{t-1},P_{t}^{S})$ is nondecreasing in $R_{t}$, $L_{t}$, $b_{t-1}^{-}$, and $b_{t-1}^{+}$.

The optimal value function $V^{*}_{t}(S_{t})$, with $S_{t}=(R_{t},L_{t},b_{t-1},P_{t}^{S})$ is nondecreasing in $R_{t}$, $L_{t}$, $b_{t-1}^{-}$, and $b_{t-1}^{+}$.

Algorithmic Technique

The traditional way to solve for the optimal value function in (13) is by backward dynamic programming. Due to the fact that this technique requires us to visit every state (which is computationally difficult), we propose the use of approximate dynamic programming. We first note that both methods require a finite state space. Since $\mathcal{R}$, $\mathcal{L}$ and $\mathcal{B}$ were assumed to be finite, we need to assume, in particular, that $\mathcal{P}$ is also finite or that it is properly discretized.

The idea behind our ADP algorithm, which we call Monotone–ADP–Bidding (see Jiang and Powell (2015)) is to iteratively learn the approximations $\mkern 1.5mu\overline{\mkern-1.5muV\mkern-1.5mu}\mkern 1.5mu^{n}_{t}(S_{t})$ (after $n$ iterations) of $V_{t}^{*}(S_{t})$ that obey the structural property of monotonicity. The algorithm is a form of asynchronous (or approximate) value iteration (AVI), so for each time $t$ in iteration $n$, only one state $S_{t}^{n}$ is visited. In addition, at each step, we perform a monotonicity preservation step to ensure the approximation is in a sense, structurally similar to $V_{t}^{*}$. We show experimentally that failure to maintain monotonicity, despite the availability of convergence proofs, produces an algorithm that simply does not work in practice.

Let $\hat{v}_{t}^{n}(S_{t}^{n})$ be an observation of the value of being in state $S_{t}^{n}$ at iteration $n$ and time $t$. Define the noise term

to be the difference between the observation and the optimal value using the iteration $n-1$ approximation. We remark, for the sake of clarity, that this is not the noise representing the deviation from the true value, $V_{t}^{*}(S_{t}^{n})$. Rather, $w_{t}^{n}(S_{t}^{n})$ is the noise from an inability to exactly observe the optimal value of the maximization: \max_{b_{t}\in\mathcal{B}}\bigl{[}C_{t,t+2}(S_{t}^{n},b_{t})+\mathbf{E}\bigl{[}\mkern 1.5mu\overline{\mkern-1.5muV\mkern-1.5mu}\mkern 1.5mu^{n-1}_{t+1}(S_{t+1})\,|\,S_{t}^{n}\bigr{]}\bigr{]}.

Before we continue, let us define a partial order $\preccurlyeq$ on the state space $\mathcal{S}$ so that for $s=(r,l,b,p)$ and $s^{\prime}=(r^{\prime},l^{\prime},b^{\prime},p^{\prime})$ where $r,r^{\prime}\in\mathcal{R}$, $l,l^{\prime}\in\mathcal{L}$, $b,b^{\prime}\in\mathcal{B}$, and $p,p^{\prime}\in\mathcal{P}$, we have that $s\preccurlyeq s^{\prime}$ if and only if the following are satisfied:

The values of any two states that can be related by $\preccurlyeq$ can be compared using Proposition 4. The main idea of the algorithm is that every observation $\hat{v}_{t}^{n}(S_{t}^{n})$ is smoothed with the previous estimate of the value of $S_{t}^{n}$ and the resulting smoothed estimate $z_{t}^{n}(S_{t}^{n})$ can be used to generalize to the rest of the state space by means of a monotonicity preserving operator, $\Pi_{M}$. Let $s\in\mathcal{S}$ be an arbitrary state that has a current estimated value of $v$. After $z_{t}^{n}(S_{t}^{n})$ is known, $\Pi_{M}$ adjusts the value of $s$ in the following way:

First, we note that if monotonicity is already satisfied, then nothing changes because in the second and third cases of (14), we get that $z_{t}^{n}\vee v=v$ and $z_{t}^{n}\wedge v=v$, respectively. If, however, monotonicity is violated, then the newly observed value $z_{t}^{n}$ prevails and replaces the previous value of $v$. Figure 4 shows an example of this operation for the two bids $b_{t-1}^{-}$ and $b_{t-1}^{+}$. In the illustration, assume that the observations are made for fixed values of $R_{t}$ and $L_{t}$, but note that when we run the algorithm, this adjustment is made over all four dimensions. The figure should be interpreted as a three–dimensional plot of the value function, where all state variables besides $b_{t-1}$ are fixed. Each bid pair, $b_{t-1}=(b_{t-1}^{-},b_{t-1}^{+})$, is associated with a $z$–coordinate value represented by colors in gray scale (darker colors correspond to larger values). In the first and third plots, new observations arrive, and in the second and fourth plots, we see how the $\Pi_{M}$ operator uses monotonicity to generalize the observed values to the rest of the state space.

The stepsize sequence used for smoothing in new observations with the previous approximation is denoted $\alpha_{t}^{n}$, which can be thought of as a (possibly stochastic) sequence in $n$, for each $t$. Furthermore, states that are not visited do not get updated unless the update is made through the operator $\Pi_{M}$, so we also define:

For notational purposes, let us also define the history of the algorithm until iteration $n$ by the filtration

2. Algorithm Description and Convergence

The full description of the algorithm is given in Figure 5.

Monotone–ADP–Bidding can be shown to converge; we reproduce the set of assumptions and the resulting theorem here.

i.e, every state is visited infinitely often (see the Extended Borel–Cantelli Lemma in Breiman (1992)).

The optimal value function $V_{t}^{*}(s)$ and the observations $\hat{v}_{t}^{n}(S_{t}^{n})$ are bounded above and below by $\pm V_{\text{{max}}}$, where $V_{\text{{max}}}$ is a positive constant.

The noise sequence $w_{t}^{n}$ satisfies \mathbf{E}\bigl{[}w_{t}^{n+1}(s)\,|\,\mathcal{F}^{n}\bigr{]}=0.

For each $t\leq T$ and state $s$, suppose $\alpha_{t}^{n}\in$ is $\mathcal{F}^{n}$–measurable and

$\displaystyle\sum_{n=0}^{\infty}\alpha_{t}^{n}(s)=\infty\quad a.s.$,

$\displaystyle\sum_{n=0}^{\infty}\alpha_{t}^{n}(s)^{2}<\infty\quad a.s.$

Under Assumptions 1–4, for each $t\leq T$ and $s\in\mathcal{S}$, the estimates $\mkern 1.5mu\overline{\mkern-1.5muV\mkern-1.5mu}\mkern 1.5mu_{t}^{n}(s)$ produced by the Monotone–ADP–Bidding Algorithm of Figure 5, converge to the optimal value function $V_{t}^{*}(s)$ almost surely.

The proof is based on the result for a generalized MDP with a monotone value function in Jiang and Powell (2015). ∎

3. Approximating the Expectation

Our algorithm can be applied to any model of spot prices $P_{t}$, with the caveat that more complex models generally require a higher dimensional state space. These include diffusion models (i.e., Schwartz (1997), Cartea and Figueroa (2005), Coulon et al. (2013)), which often incorporate features such as Ornstein–Uhlenbeck processes, jump terms, and regime switching. Recently, there has also been interest in structural models of electricity prices, where the underlying supply, demand, and bid–stack behavior is taken into account; see Carmona and Coulon (2014) for a survey.

The fact that the state space becomes larger or higher dimensional is indeed a computational difficulty that requires the availability of more powerful computational resources, but the convergence of the algorithm is unaffected (as long as $P_{t}^{S}$ is properly discretized). On the other hand, any model without finite support (or finite, but with large cardinality) necessitates the approximation of an expectation using a sample mean in Step 2a of the Monotone–ADP–Bidding algorithm (see Figure 5). In other words, because the expectation \mathbf{E}\bigl{(}\mkern 1.5mu\overline{\mkern-1.5muV\mkern-1.5mu}\mkern 1.5mu^{n-1}_{t+1}(S_{t+1})\,|\,S_{t}^{n}\bigr{)} of Step 2a is, in general, impossible to compute, we must resort to letting $v_{t}^{n}(S_{t}^{n})$ be the solution to the sample average approximation (see Kleywegt et al. (2002)) problem:

where $S_{t+1}^{j}$ are samples drawn independently from the distribution $S_{t+1}\,|\,S_{t}^{n}$. Suppose we take the observation $\hat{v}_{t}^{n}(S_{t}^{n})$ to be the value of (15). By an interchange of the conditional expectation and the max operator, we see that:

and thus, after conditioning on $\mathcal{F}^{n-1}$ on both sides, we see that $\mathbf{E}[w_{t}^{n}(S_{t}^{n})|\mathcal{F}^{n-1}]$ is biased upward from zero, a contradiction of Assumption 3. When $J$ is large, we can certainly solve the sample average approximation problem in Step 2a and apply the algorithm as is and expect an effective heuristic solution. Practically speaking, our informal tests (using $J=1000$ on a diffusion price model) showed no significant convergence issues. Even so, we cannot claim that such an approximation produces a theoretically sound and convergent algorithm due to the biased noise term. This calls for us to propose another version of Monotone–ADP, for which Assumption 3 can be easily satisfied, without restricting to price process models that facilitate an easily computable expectation of the downstream value.

4. Post–Decision, Distribution–Free Approach

Using the idea of a post–decision state (see Powell (2011)), we can make a small adjustment to the algorithm, so that Assumption 3 is satisfied. In the case of the hourly bidding problem, the post–decision state $S_{t}^{b}$ is the state–action pair $(S_{t},b_{t})$. Oftentimes, post–decision states help simplify the computational aspect of an MDP, but unfortunately, for this problem instance, the post–decision state space is higher dimensional than the pre–decision state space. Let $S_{t}^{b}=(S_{t},b_{t})\in\mathcal{S}^{b}$ and define the post–decision value function

Notice that we can rewrite Bellman’s optimality equation as:

Instead of attempting to learn $V_{t}^{*}$, the idea now is to algorithmically learn the post–decision value function $V_{t}^{b}$ using the relation (23) and to implement the policy by solving

Not surprisingly, the post–decision value function $V_{t}^{b}$ also satisfies a monotonicity property, over six dimensions.

The post–decision value function $V^{b}_{t}(S_{t}^{b})$, with $S_{t}^{b}=(R_{t},L_{t},b_{t-1},b_{t},P_{t}^{S})$ is nondecreasing in $R_{t}$, $L_{t}$, $b_{t-1}^{-}$, $b_{t-1}^{+}$, $b_{t}^{-}$, and $b_{t}^{+}$.

Let $\mkern 1.5mu\overline{\mkern-1.5muV\mkern-1.5mu}\mkern 1.5mu_{t}^{b,\,n}$ be the iteration $n$ approximation of the post–decision value function, $S_{t}^{b,\,n}$ be the state visited by the algorithm in iteration $n$, $\hat{v}_{t}^{b,\,n}(S_{t}^{b,\,n})$ be an observation of $V_{t}^{b}(S_{t}^{b,\,n})$ using the iteration $n-1$ approximation, $w_{t}^{b,\,n}(S_{t}^{b,\,n})$ be the observation noise, $\mathcal{F}^{b,\,n}$ be a filtration defined analogously to $\mathcal{F}^{n}$, and $\Pi_{M}^{b}$ be the monotonicity preservation operator on $\mathcal{S}^{b}$ defined analogously to $\Pi_{M}$. More precisely,

where $s=(r,l,b_{1},b_{2},p)\preccurlyeq^{b}s^{\prime}=(r^{\prime},l^{\prime},b_{1}^{\prime},b_{2}^{\prime},p^{\prime})$ with $r,r^{\prime}\in\mathcal{R}$, $l,l^{\prime}\in\mathcal{L}$, $b_{1},b_{1}^{\prime},b_{2},b_{2}^{\prime}\in\mathcal{B}$, and $p,p^{\prime}\in\mathcal{P}$ if and only if

The new algorithm for post–decision states is shown in Figure 6, and a set of analogous assumptions are provided below. We remark that by definition, V_{T-1}^{b}(S_{T-1}^{b})=\mathbf{E}\bigl{[}C_{\textnormal{term}}(S_{T})\,|\,S_{T-1}^{b}\bigr{]}; thus, we only need to loop until $T-2$ in Step $2$.

For all $s\in\mathcal{S}^{b}$ and $t\leq T$,

The optimal post–decision value function $V_{t}^{b}(s)$ and the observations $\hat{v}_{t}^{b,n}(S_{t}^{b,n})$ are bounded above and below, by $\pm V_{\text{{max}}}$.

The noise sequence $w_{t}^{b,n}$ satisfies \mathbf{E}\bigl{[}w_{t}^{b,n+1}(s)\,|\,\mathcal{F}^{b,n}\bigr{]}=0.

The advantage to applying this revised algorithm is that even when we cannot compute the expectation in Step 2a and must rely on sample paths, we can still easily satisfy Assumption 7 (the unbiased noise assumption), unlike in the pre–decision case. To do so, we simply use:

where we transition from $S_{t}^{b,n}$ to $S_{t+1}^{n}$ using a single sample outcome of prices $P_{(t,t+1]}$. Hence, the noise term $w_{t}^{b,n}(S_{t}^{b,n})$ is trivially unbiased.

Step 2a of the algorithm (Figure 6) can thus be rewritten as:

Under Assumptions 4–7, for each $t\leq T$ and $s\in\mathcal{S}^{b}$, the estimates $\mkern 1.5mu\overline{\mkern-1.5muV\mkern-1.5mu}\mkern 1.5mu_{t}^{b,n}(s)$ produced by the post–decision version of Monotone–ADP–Bidding Algorithm of Figure 6, converge to the optimal post–decision value function $V_{t}^{b}(s)$ almost surely.

Before discussing the proof, we state two necessary lemmas (proofs in Appendix A). The idea of the first lemma is attributed to Tsitsiklis (1994).

Define deterministic bounding sequences $L_{t}^{k}$ and $U_{t}^{k}$ in the following way. Let $U^{0}=V^{*}+V_{\text{{max}}}\cdot e$ and $L^{0}=V^{*}-V_{\text{{max}}}\cdot e$, where $e$ is a vector of ones. In addition, $U^{k+1}=(U^{k}+HU^{k})/2$ and $L^{k+1}=(L^{k}+HL^{k})/2$. Then, for each $s\in\mathcal{S}^{b}$ and $t\leq T-1$,

$U^{k}$ and $L^{k}$ both satisfy the monotonicity property: for each $t$, $k$, and $s_{1},s_{2}\in\mathcal{S}^{b}$ such that $s_{1}\preccurlyeq^{b}s_{2}$,

With Lemmas 1 and 21, we can proceed to show convergence of the post–decision state version of Monotone–ADP using the general steps to prove convergence of Monotone–ADP for pre–decision states taken in Jiang and Powell (2015). The steps are as follows:

Given a fixed $k$ and a state $s\in\mathcal{S}^{b}$ such that $s$ is increased finitely often by the monotonicity preservation operator $\Pi_{M}^{b}$, then we can show that for any sufficiently large $n$,

There exists at least one such state, i.e., the minimal state $(0,0,(b_{\text{{min}}},b_{\text{{min}}}),P_{t}^{S})$. Repeat the argument for states that are decreased finitely often by $\Pi_{M}^{b}$.

Next, we must show that states $s$ that are affected by $\Pi_{M}^{b}$ infinitely often also satisfy (22). This leverages the fact that the result has already been proven for states that are affected finitely often. The idea is that if all states immediately less than $s$ (i.e., $x$ is immediately less than $y$ if $x\preccurlyeq^{b}y$ and there does not exist $z$ such that $x\preccurlyeq^{b}z\preccurlyeq^{b}y$) satisfy (22), then $s$ satisfies (22) as well. Lemma 21 and an induction argument are used in this part of the proof.

Finally, combining Lemma 1 along with the fact that all post–decision states $s\in\mathcal{S}^{b}$ satisfy (22), it is easy to see that from a type of squeeze argument,

Note that both Steps (1) and (2) require Assumption 7, hence the focus that we have placed on it in this paper. ∎

5. Stepsize Selection

The selection of the stepsize $\alpha_{t}^{n}$, also known as a learning rate, can have a profound effect on the speed of convergence of an ADP algorithm. A common example of stepsize rule that satisfies Assumption 4 is simply:

where $N(S_{t}^{n},n)=\sum_{m=1}^{n}\mathbf{1}_{\{S_{t}^{m}=S_{t}^{n}\}}$ is the number of visits by the algorithm to the state $S_{t}^{n}$. The issue is that this method weighs all observations equally, even though we know that the error can be extremely large in early iterations of any ADP algorithm. See Chapter 11 of Powell (2011) for an overview of the numerous available stepsize rules.

After some experimentation, we found that the bias–adjusted Kalman Filter (BAKF) developed in George and Powell (2006), performed better than simpler alternatives. The main idea behind BAKF is to choose $\alpha_{t}^{n}$ such that the mean squared error to the true value function is minimized; we omit the details and refer interested readers to the original paper.

Benchmarking on Stylized Problems using Pre–Decision Monotone–ADP

In this section, we present results of running Monotone–ADP–Bidding and traditional approximate value iteration on a tractable problem (i.e., the optimal solution is computable) in order to show the advantages of using $\Pi_{M}$. In this section, we consider both four and five dimensional versions of the sequential bidding problem. We first describe some simplifications to make benchmarking possible.

In order to benchmark the algorithm against a truly optimal solution, we make some simplifying assumptions (to be relaxed in the following section) so that backward dynamic programming can be used to compute an optimal solution. First, we suppose that $P_{t}$ has finite support and that $M=1$, so that the exact value of \mathbf{E}\bigl{[}V_{t+1}(S_{t+1})\,|\,S_{t}\bigr{]} can be computed easily. When $M$ is larger, we can only compute an approximation to the expectation, due to the fact that an exponential in $M$ number of outcomes of the price process need to be considered for an exact result.

In addition, in the numerical work of this paper, we take the traditional approach and choose $C_{\textnormal{term}}(s)=0$; however, we remark that this may not always be the best choice in practice. See Section 6.3 for further discussion on the issue of selecting a terminal contribution function.

To test the approximate policies, we compute a value of the policy in the following way. For a particular set of value function approximations $\mkern 1.5mu\overline{\mkern-1.5muV\mkern-1.5mu}\mkern 1.5mu$, the set of decision functions can be written as

be a sample outcome of the revenue. We report the empirical value of the policy, which is the sample mean of F\bigl{(}\mkern 1.5mu\overline{\mkern-1.5muV\mkern-1.5mu}\mkern 1.5mu,\omega\bigr{)} over $1000$ sample paths $\omega$.

First, we consider a four dimensional variation of the bidding problem, where $S_{t}=(R_{t},L_{t},b_{t-1}^{-},b_{t-1}^{+})$. In particular, we assume that the price process has no state variables. Several versions of this problem are explored by altering the parameter values: a typical size for the batteries under consideration for the energy arbitrage application is $R_{\text{{max}}}=6$ MWh, but we also allow values of $R_{\text{{max}}}=12$ MWh and $R_{\text{{max}}}=18$ MWh for variety. The decision space is fixed in the following way: we set $b_{\text{{min}}}=15$ and $b_{\text{{max}}}=85$, and discretized linearly between $b_{\text{{min}}}$ and $b_{\text{{max}}}$ for a total of 30 possible values in each dimension of the bid. The price process $P_{t}$ has the form

where the sinusoidal (representing the hour–of–day effects on price) deterministic component is

and $\epsilon_{t}\in\{0,\pm 1,\pm 2,\ldots,\pm 20\}$, a sequence of mean zero i.i.d. random variables distributed according to the discrete pseudonormal distribution with $\sigma_{X}^{2}=49$ (a discrete distribution where the probability masses are defined by the evaluating at the density function of $\mathcal{N}(0,\sigma_{X}^{2})$ and then normalizing). We consider both the cases where the battery age does and does not matter (by setting $\beta(l)=1$), in effect introducing an irrelevant state variable. When aging does matter, the aging function we use is $\beta(l)=(l/L_{\text{{max}}})^{\frac{1}{6}}$, which provides a roughly linear decline in efficiency from 100% to around 70%, followed by a much steeper decline. Lastly, in Problem 4, we considered a uniform distribution for the noise, while the remaining problems used pseudonormal noise. In line with the operation procedures of the NYISO, the undersupply penalty parameter $K$ is set to $1$ in our simulations — this means that if one is unable to deliver energy to the market, then the penalty is precisely the current spot price (essentially, we are paying another generator to produce the energy instead). The different problems instances, labeled $A_{1}$–$F_{1}$, along with their state space cardinalities are summarized in Table 1.

2. Numerical Results for Variation 1

We first evaluate the effectiveness of Monotone–ADP–Bidding versus approximate value iteration, a traditional ADP algorithm (exactly the same as Monotone–ADP–Bidding with $\Pi_{M}$ removed); the results for Variation 1 are given in Table 2.

Figure 7 gives a quick visual comparison between the two types of approximate value functions, generated by approximate value iteration and Monotone–ADP–Bidding. We remark that after $N=1000$ iterations, the value function approximation in Figure 7(b) obtained by exploiting monotonicity has developed a discernible shape and structure, with a relatively wide range of values. The result in Figure 7(a), on the other hand, is relatively unusable as a policy.

We notice that as the cardinality of the state space increases, the value of monotonicity preservation becomes more pronounced. This is especially evident in Problem $F$, where after $N=1000$ iterations, Monotone–ADP–Bidding achieves 45.9% optimality while traditional approximate value iteration does not even reach 10%. Although this finite state, lookup table version of approximate value iteration (for lookup table) is also a convergent algorithm (see Proposition 4.6 of Bertsekas and Tsitsiklis (1996)), its performance is markedly worse, especially when the state space is large. Because it exploits the monotone structure, Monotone–ADP–Bidding has the ability to quickly attain the general shape of the value function. Figure 8 illustrates this by showing the approximations at early iterations of the algorithm. These numerical results suggest that the convergence rate of the ADP algorithm is substantially increased through the use of the monotonicity preserving operation.

With the effectiveness of Monotone–ADP–Bidding on Variation 1 established, we now examine its computational benefits over backward dynamic programming. A comparison of CPU times between Monotone–ADP–Bidding and backward dynamic programming is shown in Figure 9, where the horizontal axis is in log–scale. Once again, we notice the order of magnitude difference in computation time for the exact solution and for the near–optimal ADP solution. Indeed from Table 3, we see that we can achieve very good solutions using an ADP approach while cutting computational resources by over 93%. In the most drastic case, Problem $F$ (over 150,000 states), a 95% optimal solution is achieved using only 4% the amount of computational power.

3. Variation 2

Briefly, we also consider a problem with a more complex price process: a Markov Regime–Switching model with two regimes, denoted by the process $X_{t}$. We represent the normal regime as $X_{t}=0$ and the spike regime as $X_{t}=1$. Let $S(t)$ be a deterministic seasonal component, $\epsilon_{t}$ be discrete, independent and identically distributed (i.i.d.) random variables representing noise in the normal regime, and $\epsilon^{s}_{t}$ be discrete, i.i.d. random variables representing noise in the spike regime. The price process can be written as:

Also, we define the transition probabilities of the (time–inhomogenous) Markov chain $X_{t}$:

Because $X_{t}$ only takes two states, let $p(t)=p_{0,1}(t)$ (the probability, at time $t$, of moving from the normal regime into the spike regime) and $q(t)=p_{1,0}(t)$ (the probability, at time $t$, of returning to the normal regime). The state variable for this problem is five–dimensional: $S_{t}=(R_{t},L_{t},b_{t-1}^{-},b_{t-1}^{+},X_{t})$. In order to generate a small library of test problems, we considered two versions of the seasonal component:

for $i\in\{1,2\}$ and $f_{1}(x)=\sin(x)$ and $f_{2}(x)=\cos(x)$. We roughly model the fact that price spikes tend to occur more frequently when demand is high. Since demand is often modeled using sinusoidal functions, we use the following for $p(t)$ (the probability of moving from the normal regime to the spike regime) when the seasonal component is $S_{i}(t)$:

for some parameter $\alpha_{p}\leq 1$, representing the maximum probability of moving to the spike regime: $p_{t}(t)\in[0,\alpha_{p}]$. In these numerical results, $q(t)$, the probability of returning to the normal regime, is always modeled as a constant $\alpha_{q}$. Moreover, both $\epsilon_{t}$ and $\epsilon_{t}^{s}$ have support $\{-10,-9,-8,\ldots,+39,+40\}$ and are distributed according to the discrete pseudonormal distribution (with parameters $(\mu_{X},\sigma_{X})=(0,7)$ and $(\mu_{X},\sigma_{X})=(15,20)$, respectively). The skewed support allows us to model the preponderance of upward spikes in electricity spot prices. The remainder of the parameters vary across the test problems and are summarized in Table 4 below.

We ran both Monotone–ADP–Bidding and traditional approximate value iteration for 10,000 iterations on each of the test problems. The results of the benchmarking are summarized in Table 5 below (for brevity, we omit plots of the approximate value function and computation times and instead state that they are very comparable to those of Variation 1). It is clear that, once again, Monotone–ADP–Bidding provides significantly better solutions than approximate value iteration, particularly in the early iterations.

Case Study: Training and Testing an ADP Policy Using Real NYISO Data

In this section, we use the distribution–free, post–decision state version of Monotone–ADP to produce bidding policies for the New York City zone of the NYISO, with the goal of demonstrating the idea of training using only historical data as “sample paths.” The case study uses two full years of 5–minute real–time price data obtained from the NYISO, for the recent years of 2011 and 2012. See Figure 10 below for a visual comparison.

The concern with the historical prices is that we must satisfy Assumption 6; i.e., we must assume that the unknown stochastic process $P_{t}$ is bounded. This is not an unreasonable assumption if we allow the bound to be high, say $3,000, which is consistent with the prices in our data set. We remark again that, in order to satisfy Assumption 7, we use:

in Step 2a of Figure 6, where the transition from $S_{t}^{b,n}$ to $S_{t+1}^{n}$ is accomplished using a single sample from historical data. The remaining assumptions are satisfied for the same reasons as before.

There are many sensible ways to choose training data for a specific operating time horizon. In this paper, we consider two commonsense methods: 1) using historical samples from the same month of the previous year to train a policy for the current month (“ADP Policy 1”), and 2) using samples from the previous month to train a policy for the current month (“ADP Policy 2”). The rationale for the first method is that the price process may behave similarly in the same month across years (though factors like weather, natural gas prices, etc, should be considered before assuming that such a statement is true), and the rationale for the second method is to simply use the most recent data available. We train an hourly bidding policy that has a horizon of one day ($T+1=24$) and the post–decision state variable for this case study is five–dimensional:

where the bids are linearly discretized between $b_{\text{{min}}}=0$ and $b_{\text{{max}}}=150$ into 15 grid points in each dimension. Although it is difficult to discern from Figure 10, 98.2% of the prices in our data set are below $150. To have a lower dimensional state variable for more reasonable runtimes, we elect to assume$P_{t}^{S}=\{\}$and$\beta(l)=1$(it is also not typical for a battery manufacturer to provide an expression for$\beta(l)$; an accurate model for$\beta(l)$would require estimation from empirical data, which is outside the scope of this paper). Conversations with industry colleagues suggested that, for this application, it is reasonable to model a 1 MW, 6 MWh battery. Since$M=12$, we choose$R_{\text{{max}}}=72$, giving us a state space of size$|\mathcal{S}^{b}|=3.6$ million states, much larger than that of the benchmark problems in the previous section. The remaining details are summarized in the list below.

Due to the fact that the characteristics of the spot prices can be very different on weekends (see e.g., Coulon et al. (2013)), we considered weekdays only. In a true application, it would be important to train a separate policy for weekends.

In order to have a larger data set for our simulations, our main assumption is that spot prices of a particular hour are identically distributed across weekdays of the same month, allowing us to train and test on a large set of sample paths.

We train a daily value function for each month of the year. In essence, we combine the data for the weekdays of each month to produce a policy that is valid for any given weekday of the same month.

As before, we set the undersupply penalty parameter $K$ to 1.

The real–time prices from 2011 are used as training data and the prices from 2012 are used simultaneously as training and test data: for each month of 2012, we generate two policies, one trained using data from the same month in 2011 and the other trained using data from the previous month. The revenues generated by these policies are given in Table 6, where the evaluation method from Section 5 is used. The results correspond to running the algorithm for $N=100{,}000$ iterations. Note that because the post–decision version does not compute an expectation, each iteration is significantly faster than that of the pre–decision version, but in general, requires more iterations. The results show that ADP Policy 1 (training on data from the same month of the previous year) narrowly outperforms ADP Policy 2 (training on data from the previous month) in most cases. Although our MDP only optimizes for revenue in expectation, we nevertheless report that the (0.05–quantile, 0.95–quantile) of daily revenue for ADP Policy 1 is ($60.37,$474.24) with a median of $174.55. For ADP Policy 2, we have a (0.05–quantile, 0.95–quantile) of daily revenue of ($44.13, $453.30) with a median of$154.12. These results confirm that the policies consistently generate revenue.

2. Comparison to Standard Trading Policies

Taking advantage of the trend that lower prices occur at night, we split the operating day hours into two intervals $1$ to $h^{*}$ and $h^{*}+1$ to 24, with $h^{*}>6$. The intervals are then sorted using average historical prices. If hour $h$ of the first interval has one of the six lowest prices, then it is designated a buy interval. Similarly, if hour $h$ of the second interval has one of the six highest prices, then it is a sell interval. All other hours are idle intervals. When placing a bid $b_{t}$, we consider the hour $h$ corresponding to $(t+1,t+2]$: if hour $h$ is a buy interval, we choose $b_{t}=(b_{\text{{max}}},b_{\text{{max}}})$; if hour $h$ is a sell interval, we choose $b_{t}=(b_{\text{{min}}},b_{\text{{min}}})$; and if hour $h$ is an idle interval, we choose $b_{t}=(b_{\text{{min}}},b_{\text{{max}}})$. This policy essentially guarantees (with the possible exception of spike situations where prices exceed $b_{\text{{max}}}$) that we fill up the battery in the interval from 1 to $h^{*}$ and then empty it in the interval from $h^{*}+1$ to 24. With some tuning, we found that $h^{*}=12$ provided the highest valued policies.

The second set of policies are again based on the idea of pairing periods of low prices with periods of high prices, but with more flexibility than policies $A_{1}$ and $A_{2}$. Instead, we sort all hours of a given day using average historical prices and designate the $k^{*}$ lowest priced hours as buy intervals, corresponding to $b_{t}=(b_{\text{{max}}},b_{\text{{max}}})$ and the $k^{*}$ highest priced hours as sell intervals, corresponding to $b_{t}=(b_{\text{{min}}},b_{\text{{min}}})$. The remaining hours are idle intervals, meaning we set $b_{t}=(b_{\text{{min}}},b_{\text{{max}}})$. Again using historical prices, at time $t$, we estimate the level of resource $\hat{R}_{t+1}$ at the beginning of the next hour as the average of the outcomes of $R_{t+1}$ over historical sample paths. When encountering a buy interval with $\hat{R}_{t+1}>60$ (nearly full battery) or a sell interval with $\hat{R}_{t+1}<12$ (nearly empty battery), we place the idle bid instead. Finally, if we detect that we have more energy in storage than can be sold in the time left until the end of horizon, we place sell bids thereafter. We report results for the tuned parameter $k^{*}=10$.

Let $\alpha<0.5$ be the parameter to our final set of policies. For each hour $h$, we compute the empirical quantiles of the historical prices at $\alpha$ and $1-\alpha$, denoted $q_{\alpha}$ and $q_{(1-\alpha)}$, respectively (note the suppressed dependence on $h$). When bidding at time $t$, we again estimate $\hat{R}_{t+1}$ using historical data. For times when the battery is estimated to be nearly full, we place the bid $b_{t}=(b_{\text{{min}}},q_{(1-\alpha)})$. Similarly, if the battery is nearly empty, we place the bid $b_{t}=(q_{\alpha},b_{\text{{max}}})$. For anything inbetween, we simply bid $b_{t}=(q_{\alpha},q_{(1-\alpha)})$, with the hope of consistently buying low and selling high. We implement the same logic for when we hold more energy than the maximum that can be sold in the time remaining and initiate a sell–off. In the numerical results below, we use $\alpha=0.1$. Smaller values of $\alpha$ correspond to the notion of reserving the battery for only the highest valued trades.

The results of running policies $A_{i}$, $B_{i}$, and $C_{i}$ are given in Table 7.

Given that they were afforded more nuanced actions than simply buy and sell, perhaps not surprisingly, Policies $C_{i}$ outperformed the rest. However, we also notice that, unlike the other policies, Policy $C_{2}$ generated large negative revenues in July–12 and November–12. Comparing Policy $C_{1}$ against ADP Policy 1 and comparing Policy $C_{2}$ against ADP Policy 2, we see the revenues generated are still a disappointing $68.5\%$ and $55.3\%$, respectively, of the ADP revenues, suggesting that it is difficult, even after tuning, for simple rule–based heuristics to perform at the level of a well–trained ADP policy that considers downstream value. Moreover, the months of July–12 and November–12 (during which Policy $C_{2}$ posted negative revenues) suggest that the ADP strategy is more robust to the differences in training data when compared to Policy $C_{i}$. A possible driving force behind Policy $C_{2}$’s failure to generate revenue during these months is that the training data from June–12 and October–12 has largely differing characteristics (e.g., many spikes) from the testing data in July–12 and November–12 (see Figure 10).

3. Additional Insights

Applying Monotone–ADP–Bidding to real data from the NYISO has given us several insights into the topic of energy arbitrage. First, we note that for both ADP Policy 1 and ADP Policy 2 (see Table 6), the largest revenues were generated in the months of May, June, July, presumably due to changes in weather. The difference between the revenues generated in the months of highest and lowest revenue, June and February, is more drastic than one might expect: \textnormal{Jun Revenue}-\textnormal{Feb Revenue}=\9{,}499.37$for ADP Policy 1 and$\textnormal{Jun Revenue}-\textnormal{Feb Revenue}=\$8{,}872.08$ for ADP Policy 2. These results suggest that perhaps energy arbitrage should not be a year–round investment, but rather one that is active only during months with potential for high revenue. As Sioshansi et al. (2009) concludes, when it comes to the value of energy storage, it is important to consider various revenue sources.

Costs of energy storage can be as low as $160 kWh-1 today, and it is reasonable to expect that they will continue to decrease. As mentioned earlier, with optimal storage control strategies and decreased capital costs, energy arbitrage can soon become profitable on its own, but as it currently stands, storage costs are still relatively high compared to potential revenue. Therefore, it is also imperative that the precise storage needs of our trading/bidding policies are well–understood; it may be the case that in some months, one would choose to dedicate the entire battery to frequency regulation, while in high revenue months, the better strategy may be to use some of the capacity toward arbitrage. It is clear that some policies, such as Policies$A_{i}$, are designed with fully utilizing the available storage in mind, but for more complex policies such as those generated by Monotone–ADP–Bidding, the usage profiles are not obvious. Figure 11 shows the empirical distribution for the storage level of the battery (on an hourly basis) throughout the test data set. Note that for presentation purposes we have scaled the plot so that the bar at$R_{t}=0$ is cut off; due to its designation as the initial state (and final state as well for most sample paths), its probability is skewed to 0.09 and 0.10, for the two plots respectively. The high probability at 1 MWh is likely explained by the fact that it corresponds to full hourly charge, the minimum amount of energy needed to avoid the possibility of an undersupply penalty.

The 0.9– and 0.95–quantiles for ADP Policy 1 occur at 3.00 MWh and 3.75 MWh, and for ADP Policy 2, they are 3.16 MWh and 3.75 MWh. This means for our (relatively short) daily trading horizon, a 6 MWh battery is unnecessary—a 33% smaller device with 4 MWh storage would have sufficed and delivered similar results at a steep discount in capital cost. However, if a longer trading horizon, say, weekly (allowing us to take into account the low prices on the weekends), is desired, it would be necessary to train a policy using a sufficiently large battery and then using simulation to determine the effective amount of storage needed by the policy. In summary, with today’s substantial capital costs, it would be prudent to do an analysis of a policy’s storage needs.

Lastly, we discuss the issue of choosing $C_{\textnormal{term}}(s)$ in a practical implementation of the algorithm. Due to the daily cycles present in the real–time market, there is likely to be little additional value added in expending computational resources towards developing a bidding policy whose horizon lasts much longer than a few days or a week. In fact, from our conversations with industry colleagues, we envision that a bidding policy such as ours have a daily horizon that is used repeatedly day after day, with the policy retrained periodically (perhaps weekly). For such a usage scenario, it is important to correctly choose $C_{\textnormal{term}}(s)$, because leftover energy has value that can be capitalized on even after the true horizon of the policy. We suggest the following practical methods for determining the functional form of $C_{\textnormal{term}}(s)$:

Given the knowledge that the same policy is to be reused, in an effort to prevent the forced “sell–off” type behavior that is expected when $C_{\textnormal{term}}(s)=0$, it is reasonable to choose $C_{\textnormal{term}}(s)$ to structurally resemble $V^{*}_{0}(s)$ (i.e., up to constant shifts). One strategy for accomplishing this is to first compute $V^{*}_{0}(s)$ using a zero terminal contribution, and then re–solving the dynamic program using the previously computed $V^{*}_{0}$ as the terminal contribution. This process can be iterated until the resulting policies (not the value functions themselves) are observed to converge. Our (informal) implementation of this procedure shows that the desired behavior of not forcing the storage to zero at the end of the time horizon is indeed attained.

After training an initial policy, we can determine, by inspecting the resource paths, a point in time where the storage level is empty or very low (e.g., immediately after a period of high prices). The horizon of the problem can then be redefined so that $T$ corresponds to this point in time and a new policy (with zero terminal contribution) can be trained. Essentially, we hope that the forced sell–off is translated to a point in time where a natural sell–off would have likely occurred.

Conclusion

In this paper, we describe an hour–ahead bidding and battery arbitrage problem for a real–time electricity market (e.g. NYISO’s real–time market). We then formulate the problem mathematically as an MDP and show that the optimal value function satisfies a monotonicity property, a structural result that can be exploited in order to accelerate the convergence of ADP algorithms. The algorithm that we employ is called Monotone–ADP–Bidding and uses monotonicity to infer the value of states nearby an observed state. When benchmarked against a traditional approximate value iteration algorithm, we found that the improvements in terms of solution quality were drastic. Furthermore, the ADP algorithm can reach near–optimal solutions without the need for significant computational time and power (which an exact solution technique like backward dynamic programming certainly requires); in fact, our empirical results show that near–optimal solutions can be generated using less than 10% of the computational resources necessary for backward dynamic programming. We also describe and sketch the proof of convergence for a distribution–free method where we can train value functions with Monotone–ADP–Bidding using historical spot prices — this removes the need for us to perform the difficult task of specifying and fitting an accurate stochastic model of spot prices. In our case study, the method is tested on two large data sets: the 5–minute real–time prices from the NYISO from the years of 2011 and 2012. The policies from Monotone–ADP–Bidding help us conclude that energy arbitrage may be most valuable if practiced in a select few, high revenue months. Finally, the ADP policies consistently generated more revenue than several rule–based heuristic strategies that we considered, confirming that an ADP approach that approximates future value is worthwhile.

Appendix A Proofs

it is clear that the transition from $g_{m}^{R}$ to $g_{m+1}^{R}$ is nondecreasing in the value of $g_{m}$ and nonincreasing in the value of $e_{m}^{\intercal}q_{s}$. Thus, a simple induction argument shows that for $r_{1},r_{2}\in\mathcal{R}$ and $q_{1},q_{2}\in\{-1,0,1\}^{M}$ where $r_{1}\leq r_{2}$ and $q_{1}\leq q_{2}$,

The result follows from the fact that $q(P,b)$ is nonincreasing in $b$. ∎

is nondecreasing in $g_{m}^{L}$ and nonincreasing in $e_{m}^{\intercal}d_{s}$. Like in Proposition 1, induction shows that for $l_{1},l_{2}\in\mathcal{L}$ and $d_{1},d_{2}\in\{0,1\}^{M}$ where $l_{1}\leq l_{2}$ and $d_{1}\leq d_{2}$,

The result follows from the fact that $d(P,b)$ is nonincreasing in $b$. ∎

First, we argue that the revenue function $C(r,l,P,b)$ is nondecreasing in $r$ and $l$. From their respective definitions, we can see that $\gamma_{m}$ and $U_{m}$ are both nondecreasing in their first arguments. These arguments can be written in terms of $r$ and $l$ through the transition functions $g_{m}^{R}$ and $g_{m}^{L}$. Applying Proposition 1 and Proposition 2, we can confirm that $C(r,l,P,b)$ is nondecreasing in $r$ and $l$. By its definition,

Again, applying Proposition 1 and Proposition 2 (for $m=M$), we see that the term inside the expectation is nondecreasing in $R_{t}$, $b_{t-1}^{-}$, and $b_{t-1}^{+}$ (composition of nondecreasing functions) for any outcome of $P_{(t,t+1]}$ and $P_{(t+1,t+2]}$. Thus, the expectation itself is nondecreasing. ∎

Define the function V_{t}^{b}(S_{t},b_{t})=\mathbf{E}\bigl{[}V^{*}_{t+1}(S_{t+1})\,|\,S_{t},b_{t}\bigr{]}, often called the post–decision value function (see Powell (2011)). Thus, we can rewrite the optimality equation as:

The proof is by backward induction on $t$. The base case is $t=T$ and since $V^{*}_{T}(\cdot)$ satisfies monotonicity for any state $s\in\mathcal{S}$ by definition. Notice that the state transition function satisfies the following property. Suppose we have a fixed action $b_{t}$ and two states $S_{t}=(R_{t},L_{t},b_{t-1},P_{t}^{S})$ and $S_{t}^{\prime}=(R_{t}^{\prime},L_{t}^{\prime},b_{t-1}^{\prime},P_{t}^{S})$ where $(R_{t},L_{t},b_{t-1})\leq(R_{t}^{\prime},L_{t}^{\prime},b_{t-1}^{\prime})$. Then, for any realization of the intra–hour prices $P_{(t,t+1]}$ (by Propositions 1 and 2),

with $R_{t+1}\leq R_{t+1}^{\prime}$ and $L_{t+1}\leq L_{t+1}^{\prime}$, implying that $S_{t+1}\leq S_{t+1}^{\prime}$. This means that the transition function satisfies a specialized nondecreasing property. Using this and supposing that $V_{t+1}^{*}(\cdot)$ satisfies the statement of the proposition (induction hypothesis), it is clear that $V_{t}^{b}(S_{t},b_{t})$ is nondecreasing in $R_{t}$, $L_{t}$, and $b_{t-1}$. Now, by the previous proposition, we see that the term inside the maximum of (23) is nondecreasing in $R_{t}$, $L_{t}$, and $b_{t-1}$ for any action $b_{t}$. Hence, we can take the pointwise maximum and retain monotonicity; the inductive step is complete. ∎

Previously in the proof of Proposition 4, we argued that $V_{t}^{b}(S_{t}^{b})$ is monotone in $R_{t}$, $L_{t}$, $b_{t-1}^{-}$, and $b_{t-1}^{+}$. To see the monotonicity in $b_{t}^{-}$ and $b_{t}^{+}$, first fix an outcome of $P_{(t,t+1]}$ and $b_{t},b_{t}^{\prime}\in\mathcal{B}$, with $b_{t}\leq b_{t}^{\prime}$. Observe that if we let $(R_{t+1},L_{t+1},b_{t},P_{t+1}^{S})=S^{M}(S_{t},b_{t},P_{(t,t+1]})$, then $(R_{t+1},L_{t+1},b^{\prime}_{t},P_{t+1}^{S})=S^{M}(S_{t},b_{t}^{\prime},P_{(t,t+1]})$, with only the bid dimensions changed. Therefore,

Thus, by Proposition 4, for a fixed $S_{t}$, any outcome of the price process $P_{(t,t+1]}$, and $b_{t}\leq b_{t}^{\prime}$,

Hence, after taking expectations, we get the desired result: $V_{t}^{b}(S_{t},b_{t})\leq V_{t}^{b}(S_{t},b_{t}^{\prime})$. ∎

We first show that $H$ satisfies the following properties:

$V\leq V^{\prime}\Longrightarrow HV\leq HV^{\prime}$.

$V^{*}$ is a unique fixed point of $H$, i.e., $HV^{*}=V^{*}$.

$HV-\eta e\leq H(V-\eta e)\leq H(V+\eta e)\leq HV+\eta e$, for $\eta>0$.

Statement (i) is trivial and follows directly from the monotonicity of the max and expectation operators. Statement (ii) follows from the fact that the finite horizon dynamic program exhibits a unique optimal value function (and thus, post–decision value function as well) determined by the backward recursive Bellman equations. Statement (iii) is easy to see directly from the definition of $H$. Now, applying Lemma 4.6 of Bertsekas and Tsitsiklis (1996) gives us the desired limit result. ∎

also obeys the monotonicity property. We point out that there is a small difference between the operator $H$ and $h_{t}$ in that $H$ operates on vectors of dimension $T\cdot|\mathcal{S}^{b}|$. To verify monotonicity, $s_{1},s_{2}\in\mathcal{S}^{b}$ such that $s_{1}\preccurlyeq^{b}s_{2}$. For a fixed sample path of prices $P$, let $S_{t+1}(s_{1},P)$ and $S_{t+1}(s_{2},P)$ be the respective downstream pre–decision states. Applying Propositions 1 and 2, we have that $S_{t+1}(s_{1},P)\preccurlyeq^{b}S_{t+1}(s_{2},P)$. For any fixed $b_{t+1}\in\mathcal{B}$, we apply the monotonicity of the contribution function $C_{t+1,t+3}$ (Proposition 3) and the monotonicity of $Y$ to see that