Energy markets facilitate the balancing of electricity generation (supply) and demand while ensuring non-discriminatory access. Understanding energy market dynamics is essential to improving grid e- ciency and resilience and optimizing the development of new energy conversion and storage technologies. Accurate energy price forecasts are essential for many energy storage technologies to be pro table from price arbitrage. In this paper, we apply the novel spatial-temporal dimensionality reduction method of Dynamic Mode Decomposition (DMD) to forecast 6587 locational marginal prices in the California Inde- pendent System Operator (CAISO) on the Day-Ahead Market (DAM). DMD is a promising equation-free modeling technique in systems with inherent periodic tendencies in time such as nancial markets and uid dynamics. Yet we show, for the rst time, that DMD cannot reliably forecast day-ahead energy prices due to the so-called standing wave problem. Instead, we nd Augmented DMD (ADMD) overcomes these limitations is a fast and accurate price forecaster. We benchmark DMD, ADMD, and backcasting forecasting methods for optimal price arbitrage with an energy storage system. We nd, using ADMD, a market-connected energy storage system can capture up to 92% of allowable revenues in rolling horizon simulations. Lastly, we show ADMD is up to 1000-times faster than time-series forecasting methods (i.e., ARIMA) which requires orders of magnitude less data than deep/machine learning techniques. 1 1. Introduction Wholesale energy markets are ubiquitous across the globe and provide a competitive mechanism to bal- ances electrical energy supply and demand and ensure reliable operation of large-scale power systems.[41] Energy prices encapsulate complex trade-o s between costs and exibility of diverse generation technolo- gies (e.g., coal, nuclear, natural gas, solar, wind, hydro), demand elasticity (e.g., demand response), and operational limits of high-voltage transmission systems. These constraints interact with uncertainties from demand, renewable availability, unplanned equipment outages to create strong spatial-temporal variations in energy prices.[12, 40] Figure 1 illustrates stochastic and periodic variations in wholesale prices for six nodes in the day-ahead market run by the California Independent System Operator (CAISO). In this paper, we analyze hourly time-series data obtained from http://oasis.caiso.com/ for all 6587 nodes (locations) during the calendar year 2015. Descriptive statistics for these data are reported in [12, 13]. The economics of modern electricity generation and consumption is driven by these price variations. For example, the economic case for energy storage systems is based, in part, on energy arbitrage. Through \buy low, sell high" operating policies, storage systems can monetize time uctuations in wholesale energy prices. In their seminal work, Walawalkar et al. [48] use mathematical optimization to estimate the maximum possi- ble revenue from historical market data and compare energy storage technology options. This now-standard approach has been applied by many technologies including batteries[12, 43, 16, 30, 44], and concentrated solar thermal with energy storage[13, 25, 22, 27, 14]. Likewise, time-varying energy prices incentivize large energy consumers to either directly participate in wholesale markets or demand response programs.[1, 34]. Early work focused on optimizing the schedule of large energy consumers such as our mills, air separation, and combined heat and power plants subject to time-varying electricity prices.[ mitra2013optimal , 4, 21, 5, 29] Recent work has focused on providing ancillary services with aluminum smelters, air separation plants, heat ventilation and air conditioning systems, and chemical plants.[ dowling2018economic , 50, 49, 32] Yet many of these approaches assume market forecasts are available. Figure 1: Day-ahead market (DAM) electricity prices for six nodes (locations) in the California from January 25th28th, 2015. DAM prices exhibit strong daily trends typi ed with valleys (minimums) at midday and at night and peaks (maximums) in the morning and evening. Both time-series and machine learning techniques have been used to extract trends from massive market datasets and make accurate forecasts. For example, Conejo et al. demonstrate autoregressive integrated moving average (ARIMA) models can predict 2002 Spanish DAM market prices with 10% mean squared error (MSE).[9] Likewise, Garcia et al. used generalized autoregressive conditional heteroskedasticity (GARCH) methods to predict Spanish and California DAM market prices from 2000 with 10% MSE.[18] These classical time-series models require tedious hypertuning (e.g., feature selection, residual analysis), especially to capture periodic trends. As such, recent work focuses on applying machine learning algorithms for DAM forecasting [2, 17, 45, 6, 33, 35, 39, 3, 35, 26]. For example, Lora et al. demonstrated 8% MSE for forecasting 2001 Spanish DAM prices using a weighted nearest neighbor approach. Similarly, Catal~ ao et al. used a feed forward neural network to forecast 2002 Spanish and 2000 California day-ahead market prices, achieving 8% and 3% MSE, respectively. They nd that the feedforward neural network is less computationally intensive than an ARIMA forecast. However, neural network methods are known to struggle with noise and recapitulating multiscale dynamics[39, 46]. To address these issues, Sha e-khah et al. combined wavelet transformations, ARIMA, 2 and neural networks to forecast DAM prices in Spain from 2002 [39] and demonstrated 5% MSE. Similarly, Tan et al. combined ARIMA, GARCH, and wavelet transformations to achieve near 1% MSE in forecasts for the Spanish and PJM DAM [46]. Wavelet transformations add robustness to noise and provide a systematic framework to capture diverse timescale[46, 3, 8]. All of these frameworks build individual forecasting models for each node in the electricity market, which can cause high computational times, necessitate model tuning at each node, and o ers no fundamental insights about the aggregated dynamics of energy markets. In this paper, we propose the rst application of Dynamic Mode Decomposition (DMD) to forecasting day-ahead energy market prices. DMD is a so-called equation-free dimensionality reduction technique to extract complex, nonlinear, periodic behavior from experimental and computational simulations.[37, 38, 31, 47, 24, 38, 36] DMD is best-known for its ability to discover low-dimension spatial and temporal insights from computational uid dynamics, although emerging applications include video processing and nancial (stock) markets[24]. For example, Mann and Kutz[28] proposed a stock trading algorithm build on spatiotemporal information encoded in DMD modes and demonstrate 5 times higher growth than the S&P 500 index via retrospective analysis. Likewise, Cui and Long[10] showed how DMD-based trading beat a benchmark moving-average strategy in a bull market for Chinese A-stocks. Hua et al.[20] identi ed four new reproducible modes (cycles) by analyzing stock prices between 1992 and 2014 for major companies. All of these works underscore that DMD can successfully capture the underling nonlinear dynamics in stock markets from data alone. In this paper, we show that Augmented DMD is an accurate short-term forecaster energy prices in wholesale electricity markets. ADMD requires only a modest amount of training data, much less than machine learning methods, and is orders of magnitude faster than ARIMA methods. The paper is organized as follows: Section 2 provides a mathematical primer on DMD and Augmented DMD (ADMD). Section 3 explores singular value decomposition (SVD) rank reduction, interprets dominate DMD modes, shows how DMD su ers from the so-called standing wave problem. Section 4 demonstrates the superior performance of ADMD. Section 5 demonstrates how ADMD forecasting can capture up to 92% available revenue in a simple energy arbitrage optimize case study. Finally, Section 7 summarizes conclusions and future work. 2. Mathematical Primer on DMD In this section, we review the mathematical foundations of DMD in the context of energy market forecasts. See Kutz et al.[24] and Schmid et al.[38] for additional details. 2.1. DMD Algorithm DMD starts with an input matrix Xt2Rmnwhere the rows ( m) ofXtare the spatial data features and the columns ( n) are discrete time measurements. For day-ahead energy market analysis, each row is a node (location in the transmission network), each column is an hour, and each element is a price. The price matrix Xtis then partitioned into two matrices X2Rm(n1)andX02Rm(n1), as shown in Eqs. (2) and (3). Here each xi2Rm1is a single time step of spatial data in price matrix Xt. Xt= [x1xn] (1) X= [x1xn1] (2) X0= [x2xn] (3) The overall goal of DMD is to compute a best t linear operator, A2Rmmin Eq. (4), which advances the energy price dynamics by one timestep (one hour). AXX0(4) Ais a nite dimensional approximation to the in nite dimensional Koopman operator, K, which exactly maps any measurement function ffrom one state vector, xi, to another state vector, xi+1: 3 Kf(xi) =f(xi+1) (5) For practical analysis, we desire a low rank approximation to the linear operator A. This can be accom- plished via singular value decomposition (SVD) of the data matrix X2Rm(n1): X=UV(6) HereU2RmmandV2R(n1)(n1)are unitary matrices and is the conjugate transpose. Matrix 2Rm(n1)has singular values on the diagonal ordered 12:::n1. Whenm>n1, which is often the case in DMD analysis, the bottom of are rows of zeros: U= [u1;u2;:::;um];=2 6666666666410::: 0 02::: 0 ............ 0 0:::  n1 0 0::: 0 ............ 0 0::: 03 77777777775;V= [v1;v2;:::;vn1] (7) Retaining only the rlargest singular values and vectors gives a rank rapproximation to X. The choice ofris an important tuning parameter in DMD work ows (see Section 3). The maximum number of singular values isrmax= min(m;n1). XXr=UrrV r=rX i=1uiiv i Ur= [u1;u2;:::;ur];r=2 641::: 0 ......... 0:::  r3 75;Vr= [v1;v2;:::;vr](8) Next, we calculate Ausing Eq. (4) and the pseudoinverse X1. Substituting Eqs. (6) and (8) gives: AX0X1=X0(UV)1=X0V1UX0Vr1 rUr (9) Often it is more convenient to use the low-rank projection ~A2Rrr: ~A=U rAUr=U rX0Vr1 r (10) which advances the low-dimensional projected states, ~x2Rr1, from timestep ktok+ 1: ~xk+1=~A~xk (11) High-dimensional states are reconstructed via projection with Ur2Rmr: xk=Ur~xk (12) Because ~Ais a linear operator, the price trajectories can be calculated analytically using an eigendecom- position: ~AW =W (13) Here the columns of W2Crrare eigenvectors and the diagonal of are the corresponding eigenvalues. The spatial modes, = [ 1 r]2Cmr, are then calculated via Eq. (14), and represent how prices depend on geographic location in the electric transmission network. Next, the discrete time eigenvalues, , are converted to continuous time eigenvalues, 2Crr, via Eq. (15). Here iandjare row and column indices, respectively, such that i=jare the diagonals of and . Thus encodes price time dynamics that are independent of spatial behavior. 4 =X0Vr1 rW (14) ij=(ln( ij) tifi=j 0 otherwise(15) The low-dimensional time dynamics, D(t)2Cr1, are de ned in Eq. (16). The vector b2Cr1is the initial condition of the low-dimensional system, which is obtained by solving Eq. (17), often via least-squares. D(t) =e tb (16) b=x1 (17) This concludes the decomposition of the original price data. A key advantage of DMD is that forecasting for any time and every node (location) can be quickly computed via Eq. (18). x(t) =D(t) (18) In DMD literature, reconstruction of the original data refers to evaluating Eq. (18) for all of the times in the original original Xtinput matrix. Recall that DMD assumes a low-dimensional linear system governs the data. When this assumption is violated, the reconstructed prices can di er dramatically from the actual (training) prices. We further discuss the relationship between reconstruction error and forecasting error in Section 3. 2.2. Augmented DMD Algorithm A well-documented challenge for standard DMD is modeling simple dynamics of standing waves, e.g., x(t) = sin(t) [24]. The fundamental problem for this example is that the input matrix Xis rank one, and as such DMD identi es a single mode with a positive real eigenvalue. But a sine wave is the solution to the di erential equation  x+x= 0; a conjugate pair of complex eigenvalues are needed to correctly recapitulate a sine wave. For typical days, the energy price closely resembles a simple oscillating wave (Fig. 1), which o ers one explanation for the mediocre reconstruction and forecasting errors with standard DMD. Augmented DMD (ADMD) overcomes the so-called standing wave problem by stacking time-shifted copies of input data to form the augmented input matrix which has higher rank than the original input matrix. This is analogous to applying a nite di erence rules to estimate unobserved derivates needed to model higher order systems (e.g., sine wave) and di erencing data to achieve stationary for ARIMA time-series modeling [7]. To implement ADMD, we stack l <n time-shifted copies of the data to form the augmented input matrix. The augmented matrix is then partitioned into two matrices, similar to Eqs. (2) and (3): Xaug=2 6664x1x2 xnl x2x3xnl+1 ............ xlxl+1 xn13 7775; X0 aug=2 6664x2x3xnl+1 x3x4xnl+2 ............ xl+1xl+2 xn3 7775(19) The DMD algorithm prescribed in Eqs. (2) - (18) is applied to augmented matrices Xaug2R(ml)(nl) andX0 aug2R(ml)(nl)in place of XandX0, giving eigenvalues augand modes aug. The rstmrows of augcorrespond to the current (not shifted) time and are used to forecast x(t). For ADMD, the maximum number of singular values is rmax= min(ml;nl) 3. DMD Results and Discussion We now apply DMD to CAISO energy price data. Our analysis focuses on truncation level, window size, consistency, forecasting, interpretation of DMD modes. Throughout this paper, we measure percent error using the`2norm: % error =jjXtXDMDjj2 jjXtjj2100 (20) where XDMD are either reconstructions or forecasts calculated with Eq. (18). 5 3.1. Optimal Singular Value Trunctionation to Minimize Reconstruction Error An essential tuning parameter in DMD analysis is the level of singular value truncation, especially for noisy data[24]. We begin by systematically exploring the impact r, the number of singular values retained in Eqs. (8) - (10) on reconstruction and forecasting error for DAM energy price data. It is well-known that retaining too many singular values, especially for noisy datasets, can corrupt DMD results.[24] Fig. 2 shows this trade-o holds for DAM price data for all nodes for three days ( Xt2Rm=6587n=96). Too much truncation leads to decaying reconstruction with damped dynamics, whereas not enough truncation leads to an exponentially growing reconstruction (diverges to 1) with exaggerated dynamics. The optimal truncation level was determined via sensitivity analysis and corresponds with the minimum reconstruction error of 12%, computed via Eq. (20). We reiterate these ndings are consistent with DMD literature for noisy datasets. Figure 2: Three DMD reconstructions for DAM energy prices from all nodes for January 1st to 4th, 2015 such that Xt2Rm=6587n=96. The mean price for all nodes is shown in black (dashed) and the mean DMD reconstruction is shown in red (solid). Too much truncation, r= 70 (70=95 = 73:7%), causes the DMD reconstruction to decay to constant value as t!1 . Too little truncation, r= 91 (91=95 = 96:8%) causes the DMD reconstruction to diverge as t!1 . But with the best truncation, r= 82 (82=95 = 82:3%), the DMD reconstructions recapitulates the dominant dynamics of the energy price signals. We now explore how the optimal singular value truncation level, r, changes for di erent datasets. Fig. 3 shows the sensitivity analysis of reconstruction error versus rfor three four-day windows considering all nodes (Xt2Rm=6587n=96). For all three datasets, we clearly see the trade-o between not enough and too much truncation. Moreover, the optimal truncation level (lowest reconstruction error averaged over all nodes) varies between r= 81 andr= 91, which is 81/96 = 84% and 91/96 = 95% of the full rank number of singular values. This means that standard DMD fails to extract coherent low-rank structures from the DAM energy price data. Moreover, each individual training window requires a di erent rto recover the optimal reconstruction . We will revisit these two observations in Section 4 when we discuss ADMD. Next, we quantify how DMD performance varies for individual nodes. Fig. 4 shows the distribution of reconstruction errors across all nodes for three 4-day training windows. These errors are approximately normally distributed for the January 1st- 4thdataset, but skewed-left for the January 11th- 14thdataset and bimodal for January 6th- 9thdataset. All three datasets show statistical outliers, which are likely due to random phenomena in the market clearing (price setting) process. Next, we compute average reconstruction errors for individual nodes using the optimal r(determined via sensitivity analysis) for all eight-day training windows (Xt2Rm=6587n=192) spanning January 1stthrough 30th. Fig. 5 visualizes means reconstruction errors in a histogram (all nodes) and overlaid on a map ( 2000 nodes with known latitude and longitude coordinates). We nd 75% of nodes have an average reconstruction error between 12% and 13% for January 2015. Reconstruction errors follow a right-skewed distribution, with 1.9% of nodes with 20% to 30% error. From these results, we six geographically separated nodes in California such that nodes 1 and 2 have great (less than 11%), nodes 3 and 4 have average (12%), and nodes 5 and 6 have poor (greater than 16%) reconstruction errors. Table 1 highlights ocial CAISO node names, location, and mean reconstruction 6 Figure 3: Comparison of reconstruction error aggregated over all 6587 nodes (%) versus the number of retained singular values ( r) for three 4-day training windows. The ?symbol marks the minimum reconstruction error: (left) r= 81 corresponds with 11.95% error; (middle) r= 85 gives 18.46% error; (right)r= 91 gives 16.72% error. error. For these six nodes, we reproduce the sensitivity analysis from Fig. 3. The reconstruction error versus rplots are only 1% di erent from those shown in Fig. 3. This suggests that the optimal rfor the entire dataset (all nodes) is near optimal for individual nodes. Moreover, there are no obvious geographic trends in the distribution of reconstruction errors in Fig. 5 (middle). Figure 4: Histogram of reconstruction errors tabulated per node for three training windows using r= 81 (left), 91 (middle), and 92 (right), m= 6587, and n= 96. While each window has a distinct distribution of errors, a majority of errors are between 5% and 15%. Each window shows a few nodes, which are outliers, with reconstruction errors above 15%. 7 Figure 5: (Left) Distribution of the mean reconstruction error over calendar year 2015 tabulate for each node using m=6587,n=96, and the optimal rfor each window. (Center) Spatial distribution of mean reconstruction errors for 2220 nodes in CAISO GPS coordinates publicly available on oasis.caiso.com . (Right) Location for six example nodes detailed in Table 1. Node Name Mean % Error Map Index Location General Performance CORNING 6N102 10.76 1 Corning, CA Great PRAXAIR 1B1 10.93 2 Antioch, CA Great PARADISX 6N004 12.23 3 Paradise, CA Average CERTTEED 1N002 12.23 4 Madera Acres, CA Average URBAN 6N005 16.87 5 San Diego, CA Poor PALA 6N001 17.01 6 Pala, CA Poor Table 1: We select six geographically distinct nodes, shown in Fig. 5 (right), ordered from lowest to highest mean reconstruction error. The node names are the unique identi ers used by CAISO. 3.2. Optimal Training Window Size to Minimize Reconstruction Error Besidesr, selecting the window size, i.e., the dimensions of the input data Xt2Rmn, is another important choice in DMD analysis including nancial market forecasting[28]. We repeat sensitivity analysis of reconstruction error versus rfor January 1stto 4th(n= 244 = 96) with m= 20, 100, 1000, 2000, 3000, and 6587 randomly selected nodes. Figure 6 shows that, in general, as mincreases, the minimum reconstruction decreases. This suggests that standard DMD is extracting useful multimode dynamics and correlations such that more data (larger m) improves recapitulation of trends (lower error). With only 20 nodes, reconstruction error is between 25% and 40%, which suggests there is insucient data to extract meaningful trends. We will revisit this nding in Section 4. Next, we explore how the time horizon nimpacts reconstruction error. Training windows for 2 through 15 days (n2f48;72;96;;360g) were evaluated for a rolling horizon for a whole year (2015) using all 6587 nodes in Figure 7. Each box and whisker plot in Figure 7 reports the reconstruction error with the optimalr(lowest average reconstruction error for all nodes). For the year 2015 data set, we see the median reconstruction error gradually increases with n. The lowest median error is 16% for n= 48 (2 days), which is the shortest time horizon considered. These results are consistent with intuition. DAM prices are strongly periodic with well-documented seasonal, day of the week, and daily trends.[15] A two-day window captures two full cycles. But too long of a training window attempts to learn one monolithic linear operator Ain Eq. (4) to explain dynamics across too many days. We conclude by highlight our ndings to maximize m (consider all nodes) and minimize n(short time horizon) are consistent with literature that suggests DMD inputs should be skinny , i.e., withm>>n [24]. 8 Figure 6: Sensitivity analysis of reconstruction error with respect to truncation level ( r=rmax) form= 20 tom= 6587. The nodes included in the analysis were randomly selected for m< 6587. The results show there is a clear bene t for selecting large m. Figure 7: Sensitivity analysis reveals the impact of training window size for n= 48 (2 days) to n= 360 (15 days) on reconstruction error. The center of each box marks the median reconstruction error calculated across all nodes. The edges of the box mark the quartiles. The whiskers mark the smallest and largest reconstruction errors for an individual node, after excluding outliers. n= 48 gives the smallest median reconstruction error of 16%. 3.3. From Optimal Reconstruction to Optimal Forecasts We now explore how insights for specifying r,m, andnto optimize reconstruction translate to forecast- ing. Figure 8 shows sensitivity analysis of forecasting and reconstruction errors versus truncation level r. For all three data sets considered, we observe the optimal rfor reconstruction and forecasting are similar. 9 To generalize this trend, we conduct sensitivity analysis over the entire year with 4-day training windows in a rolling horizon framework. On average, there was only a 4% di erence between rvalues that minimized reconstruction error and day 1 forecast errors. Di erences were 0.4% and 0.8% for reconstruction error versus day 2 and day 3 forecast errors, respectively. From these results, we propose a simple heuristic: truncate DMD the with best rfor reconstruction to produce for the (near) optimal forecast. Figure 8: Reconstruction and forecasting errors are compared for three 4-day training windows ( n= 96, m= 6587) for many di erent number of retained singular values r. For the rst window, training is done with data from January 1stthrough 4th. January 5th, 6th, and 7thcorrespond to the day 1, 2, and 3 forecasts, respectively. Next, we develop guidelines for the minimum training horizon length nto create stable multi-day DMD price forecasts. Computing optimal bids for participating in day-ahead energy markets, including CAISO, often require 48-hour or longer forecasts [12]. We perform a sensitivity analysis of all possible 2 to 15- day training windows in the year 2015 and record forecast errors for 1 to 4 days after the end of each training window. Fig. 9 shows forecasting error averaged over all of these windows. When the training window is shorter than the forecasting horizon, the forecasting error is large. For example, forecasting 4 days into the future using only a 1-day training horizon gives 50% average errors. In contrast, forecasting 4 days into the future with a 5-day or longer training horizon gives average errors less than 20%. This trend is consistent for all forecasting horizons considered and leads to another heuristic: choose the training horizon length nto be longer than the desired forecasting horizon. Beyond this threshold, the average forecasting error only gradually changes as nincreases. For shorter forecasts, e.g., 1 day, longer training windows slightly increase forecasting error. But for longer forecasts, e.g., 4 days, longer training windows slightly decrease forecasting errors. We emphasize that while reconstruction error is minimized with shorter training horizons (see Fig. 7), forecasting requires longer horizons. We speculate that unplanned outages, transmission constraints, uncertain electricity demands, and other factors add noise to the energy prices. This noise causes reconstruction error to increase with training window length. But, for accurate forecasts, more data is required to separate meaningful dynamics from noise, and hence longer training horizons are required. 10 Figure 9: Mean forecasting errors for 6587 nodes ( m= 6587) calculated with a rolling horizon for 1, 2, 3 and 4 days into the future. The forecasting horizon starts on January 17 thand continues through the end of the year. The SV truncation rwas performed for each windom 3.4. DMD Mode Analysis We now interpret the DMD modes, similar to previous nancial market analyses [24, 28, 20]. Speci cally, temporal modes (eigenvalues) capture price periodicity and spatial modes weight di erent locations (nodes) in the electricity market. Figure 10 shows spatial-temporal DMD modes for a 4-day training window n= 96 without truncation. In the top three plots, we see eigenvalues ordered by the prominence of the singular values. The rst six singular values explain 99% of the variability in the dataset. Interestingly, as the most important singular values (lowest number) correspond to the lowest frequency oscillations. Beyond singular value 60, this linear trend between singular value number and oscillation frequency ceases. We hypothesize the least prominent singular values are localized uctuations in energy prices caused by external factors. We nd the real eigenvalue components are both positive and negative. We do not see a clear trend with real eigenvalue components and the number of retained singular values r. The middle three plots show how the three most prominent DMD modes 1, 2, and 3are spatially distributed. We see all three modes place the least weight on prices from the San Fransisco area. The bottom three plots explore the maximum, minimum, and median contributions to 1. We see the minimum contribution has price spikes on the last day. This makes sense, as we expect the most prominent singular value to capture the average trends of the entire CA system and deemphasize localized events. 11 Figure 10: Visualization and interpretation of the DMD modes. (a)The fraction of variability that can be attributed to each singular value. The horizontal line marks 1% variability. We see that retaining the r= 6 largest singular values in the truncation explains over 99% of the variability in the dataset. (b)and(c) The imaginary and real components of the eigenvalues associated with each singular value ordered from largest to smallest. (d),(e),(f)Three most prominent spatial modes, 1, 2, and 3, respectively. Larger absolute values (red or blue) indicate a larger weight for the speci c location in the DMD mode. (i) Time-series price data from January 4thto 7thfor corresponding to the nodes with the maximum (most positive), minimum (most negative), and median, respectively, weights in 1. 4. ADMD Results and Discussion As discussed in Section 3, standard DMD, even with optimal truncation, is only able to achieve a 16% reconstruction error for the CAISO dataset from 2015. We now show that ADMD overcomes these limitations and consistently achieve less than 5% reconstruction error. 4.1. ADMD is Less Sensitivity to Truncation First, we explore the sensitivity of ADMD to truncation level rfor threen= 96 training windows containing 100 randomly selected nodes and augmentation level l= 48. Figure 11 clearly shows ADMD is more robust with 1.5 to 3 times smaller reconstruction error than DMD. ADMD reconstruction error reaches a plateau around r= 15%, suggesting ADMD successfully identi es a low-rank structure. To explain this, we veri ed the condition number of Xis 10,000 times larger than Xaug. These results support our hypothesis that standard DMD struggles to extract low-rank trends because of the standing wave problem. 12 Figure 11: Sensitivity analysis of reconstruction error versus truncation level rshows that ADMD consistently outperforms DMD. For this comparison, m= 100,l= 48, andn= 96 such that X2R96100 andXaug2R480048. DMD performs much better in Figure 6 with m= 6587 whereas this Figure considersm= 100 to benchmark against ADMD. 4.2. ADMD Performs Best Augmenting Over At Least 1 Day Next, we explore the impact of augmentation level lin Eq. 19 on ADMD performance. In Fig. 12, we see reconstruction error is not robust with l= 12 andl= 18. In contrast, reconstruction error reaches a plateau around truncation level r=rmax= 15% with l= 24,l= 36, and l= 48. We hypothesize this is because the dominant price dynamics in energy markets have 24-hour periodicity. Based on these results, we recommend l= 36,l= 48, or a similar value. Figure 12: Sensitivity analysis of ADMD reconstruction error with respect to the truncation level r=rmax wherermax= min(ml;nl). Results show l24 is needed for a robust low rank approximation with a small reconstruction error. 4.3. ADMD Does Not Bene t from Training Many Nodes Simultaneously Next, we determine if ADMD bene ts from analyzing multiple nodes simultaneously. Recall, DMD performs best when the data matrix Xis \tall and skinny", i.e., there are many more nodes than timesteps. As consequence, we show in Section 3.2 DMD works best with m= 6587 (all nodes considered). But as 13 augmentation level lincreases, the number of rows in Xaugincrease and the number of columns decrease. This means ADMD can be applied to a single node. Figure 13 compares reconstruction error for a single node and 100 randomly selected nodes over with r= 8 andl= 48 for 96 hours of training data. Surprisingly, we nd reconstruction error only improves 0.65% when forecasting many nodes simultaneously. This suggests there is less valuable information available in the spatial dynamics than previously postulated. Figure 13: Comparison of ADMD with m= 1 (red) and m= 100. Each panel is a di erent 48-hour training window. There is little bene t to reconstructing and forecasting multiple nodes simultaneously. 4.4. ADMD Outperforms DMD for Forecasting Next, we compare forecasting error from DMD and ADMD for a single node. Fig. 14 shows forecasts for January 16thand 17th, 2015 (48 hours) using the previous 96 hours for training. DMD forecasts were constructed with all 6587 nodes, whereas ADMD was only trained using the node of interest (same as Fig. 14) withr= 8 andl= 48. The results show the ADMD forecast closely tracks the realized price with around 10% forecasting error. In contrast, the DMD forecast fails to predict major temporal features and is out of phase. 14 Figure 14: Comparison of DMD and ADMD forecasts using n= 96 training data (January 12ththrough 15th, 2015) and two-days of testing data (January 16thand 17th, 2015). The ADMD forecast is superior and most importantly captures the phase of the testing data. 5. Optimal Market Arbitrage for Energy Storage Systems Accurate price forecasts are important for many aspects of grid reliability, including optimal control of individual energy systems such as conventional generators, grid-connected storage, microgrids, and virtual power plants. In this section, we benchmark the performance of deterministic model predictive control of a price-taker battery energy storage system self-scheduling into a wholesale energy market using di erent price forecasting techniques including DMD and ADMD. 5.1. Optimization Formulation We consider a storage system with time-varying energy hold-up Eiat timei. The system can purchase (charge) or sell (discharge) energy to and from the market at rates cianddi, respectively, during time interval [i;i+ 1). The energy storage operator seeks to maximize net revenue as follows: max c;d;EnX i=0xi(dici)) (21a) s:t: E i+1=Ei+cipdip8i2f0;:::;n1g (21b) wherexiis the market price forecast during time interval [ i;i+ 1). The hold-up is constrained by an energy balance and subject to round trip eciency . Variable bounds, parameter values, and units are speci ed in Table 2. Figure ??illustrates the interaction between the energy storage system and the energy market. Fig. 15 compares the solution to Eq. (21) for the market node near Chino, CA. Price 48-hour forecasts were generated with ADMD using the previous 96 hours of prices for the node of interest as training data withr= 8 andl= 48. As expected, the optimal solution is to sell energy when the forecasted price is lowest and buy energy when the forecasted price is the highest. Although price forecasts are consistently biased (e.g., higher than the actual price) the periodicity (e.g., when prices increase and decrease) is correct, leading to an e ective control strategy. 15 Variables Description Units Dimensions Range c Charging Rate MW Rn[0,1] d Discharging Rate MW Rn[0,1] E Battery Charge MWh Rn[0,4] Parameters Description Units Dimensions Rangep Round Trip Eciency - R10.88 n Length of Forecast h Z E0 Initial Battery Charge MWh R1[0,4] x Price Forecast $/MWh Rn Indicies Description Units Dimensions Range i Time index h Z [0,n] Table 2: Optimization variables and parameters used in the battery optimal control problem. 16 Figure 15: Example battery optimization solution for Chino, CA for January 7thand 8th, 2015. (Top) The historical market price and ADMD forecast both show 4 distinct price maximums which correspond with early morning (6 am) and late afternoon (3 or 4 pm). (Middle) As expected, the optimal arbitrage policy is to charge the battery (buy from the market) during low price periods and discharge (sell) at higher prices. (Bottom) The state-of-charge (SOC) in the battery is the integral of the charge and discharge signals. The black ?marks the SOC at the end of the rst day ( E24). In the rolling horizon framework, this value is the xed initial condition E0for the subsequent day. 17 5.2. Backcasting Benchmark We use simple backcasting as a benchmark forecasting strategy[43]. In backcasting, prices for the n most recent observed days as used as forecasts for the next ndays. We compare 48-hour rolling horizon forecasting errors for backcasting, DMD (48 hours of training data, optimal truncation), and ADMD (96 hours of data using, r= 8,l= 48) averaged over all 6587 nodes in CAISO for January 5th, 2015 to December 7th, 2015. Surprisingly, we found median forecasting errors of 17%, 26%, and 31% for backcasting, ADMD, and DMD respectively. Below we show the lowest median forecasting error does not always correspond with the maximum realized revenue. Sioshansi et al observe backcasting is extremely e ective because it recapitulates the periodicity of day-ahead energy market prices. [42] 5.3. Rolling Horizon Computational Results We now compare the performance of DMD, ADMD, and backcasting for energy arbitrage by repeatedly solving Eq. (21) in a rolling horizon framework for Januaryth, 2015 to December 9th, 2015. Each instance of Eq. (21) is rst solved with a 48 hour prediction horizon. The computed policy is then simulated over the rst 24 hours and pro ts are recorded with the real (historical) price. The simulation window is then advanced 24 hours. The value for E24for the previous day is copied to E0for the next day. We start by considering all 6587 nodes in the CAISO market. For DMD, we use a 48-hour training window (Xt2R486587) with optimal truncation. For ADMD, we generate forecasts individually for each node using r= 8 andl= 48 with 96 hours of training data (thus 6587 unique instances of Xaug2R4848). Backcasting was performed using the previous 48 hours of data. Figure 16 compares the market revenues for optimal arbitrage with perfect information, backcast, DMD, and ADMD price forecasts. The perfect information scenario, which uses historical prices in Eq. (21), represents the maximum available revenue. We see hot spots in the perfect information map, which show the greatest potential for arbitrage is in central California as well as the coasts. We nd that DMD performs poorly, yielding only 24% of the perfect information revenue. This is expected as DMD su ers from the so-called standing wave problem. In contrast, ADMD and backcasting capture 81% and 84% of the perfect information revenue, respectively. Backcasting and ADMD capture the hot spot in central California but miss the opportunities along the coast. These results are summarized in Table 3. Next we identi ed 2209 nodes (2209 / 6587 = 33.5% of the original data set) for which the prices were always between 0 $ / MWh and 200 $ / MWh during calendar year 2015. We refer to these as Select Nodes . We nd that ADMD and backcasting recover 92% and 88% of potential revenues, respectively, for only select nodes. This is a modest increase compared to all nodes. DMD, however, recovers less relative pro t when only considering select nodes. We hypothesize this is because DMD bene ts from forecasting as many nodes as possible. Backcasting DMD ADMD All Nodes 83.7% 24.3% 80.8% Select Nodes 88.3% 21.7% 92.2% Table 3: Comparison of net revenues, relative to the perfect information case, for each forecasting method for January 5thto December 9thin 2015. Select nodes are locations where the market price was always between 0 $ / MWh and 200 $ / MWh. 18 Figure 16: Simulated net pro ts for energy arbitrage for all nodes in the CAISO day-ahead market in 2015. DMD forecasting greatly underperforms relative to ADMD and backcasting. ADMD and backcasting perform well through most of CAISO, but fail to capture the largest revenue opportunities across the coast. The perfect information reference shows the maximum available revenue. 6. Computational Timing Results A key advantage of DMD and ADMD is computational speed. Training is quick, often limited by the SVD calculation, and forecasting is almost instantaneous using the explicit formula in Eq. (18). To demonstrate these bene ts, we performed simple computational benchmarks using an Auto-Regressive Integrated Moving Average (ARIMA) model implemented in the StatsModel package. We considered 48 hours of training data and ARIMA model form ( p;d;q ) = (7;1;1). Training and forecasting with an ARIMA model took 33 seconds per node on average. In contrast, it took 11.2 seconds in total for ADMD ( r= 8,l= 48,n= 96) to train and forecast for all 6587 nodes, which is 0.002 seconds per node. On a per-node basis, ADMD is 20,000 times faster than ARIMA. Our custom implementation of DMD using SciPy is available on GitHub ( https://github. com/dowlinglab/dmd-energy-markets ). These benchmark tests were performed on a 2014 MacBook Pro (3.0 GHz Intel Core i7 processor, 16 GB of RAM) using an Anaconda Python environment. 7. Conclusions and Future Work In this paper, we extend the success of DMD in forecasting nancial markets[28, 10] to wholesale elec- tricity markets. We nd that standard DMD is unable to extract coherent low-rank structures in California Independent System Operator (CAISO) day-ahead market prices for the calendar year 2015. Instead, we 19 nd Augmented DMD (ADMD), which is equivalent to di erencing in a time-series model, overcomes the so-called standing wave problem and can consistently reconstruct and forecast with errors less than 5% with a low-rank approximation. We benchmark DMD, ADMD, and backcasting in a rolling horizon model predictive control case study for energy arbitrage with a market connected battery. We show that ADMD outperforms backcasting by capturing 92% versus 88% of available net market revenue at nodes without extreme prices. We also nd that ADMD is up to 20,000 times faster than ARIMA for training and forecasting. There are several opportunities to further enhance DMD for energy market forecasting. Denoising pro- cedures can greatly improve the performance of DMD with experimental data and should be exploring for energy price forecasting.[11, 19] Multiresolution DMD may better identify and recapitulate dynamics at disparate time-scales in electricity markets.[23] Both of these methods may help explain why ADMD as currently implemented does not bene t from forecasting multiple nodes simultaneously. We also plan to more thoroughly benchmark DMD against other forecasting techniques, including emerging machine learning methods, to understand trade-o s between data requirements, computational time, need for model tuning, and forecasting error. Real-time markets (RTMs) set prices in 5 to 15-minute increments and are much more volatile than day-ahead markets.[13] Studying the performance of DMD for RTMs and other geographic regions is another future opportunity.