{'fraction_non_alphanumeric': 0.031668564611463584, 'fraction_numerical': 0.012237225023950256, 'mean_word_length': 4.566914872207688, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Online Decision Making for Trading Wind Energy

This paper proposes and develops a new algorithm for trading wind energy in electricity markets, within an online learning and optimization framework. In particular, we combine a component-wise adaptive variant of the gradient descent algorithm with recent advances in the feature-driven newsvendor model. This results in an online offering approach capable of leveraging data-rich environments, while adapting to non-stationary characteristics of energy generation and electricity markets, and with a minimal computational burden. The performance of our approach is analyzed based on several numerical experiments, showing both better adaptability to non-stationary uncertain parameters and significant economic gains.


Problem statement
Traditionally, the way in which trading wind energy has been considered relied on a two-step approach, with first the predictive modeling of future energy generation (within either deterministic or probabilistic frameworks), and then the use of such forecasts as input to expected utility maximization strategies (possibly also accommodating risk aversion). Although fruitful, these methodologies may be computationally expensive (especially the forecasting part), while the value of final decisions is highly affected by the quality of the forecasts employed. It may then be beneficial to integrate the forecasting and decisionmaking stages, also within a data-rich and non-stationary environment thanks to online learning and optimization. On a different front, electricity markets are amid rapid transformations towards reducing granularity and lead times, facilitating the integration of non-dispatchable energy sources but increasing the computational and adaptability requirements of the offering algorithms.
Online learning can help on both fronts. On the one hand, online learning algorithms free the decision-maker from most assumptions about the wind or market dynamics. On the other hand, online learning algorithms are typically efficient methods capable of adapting to the increasing computational needs. Furthermore, the online learning analysis is based on regret as opposed to the classical maximization of the expected utility, possibly allowing to derive additional insights into the properties of trading strategies.

From renewable energy trading to online learning
Most wind energy is traded in wholesale electricity markets (referred to as forward markets in this paper), where an offer is submitted prior to the actual delivery of energy. However, the stochastic nature of wind energy entails incurring deviations from the original offer. There are countless ways of approaching this problem depending on the market structure and how uncertainty is addressed, and therefore it is infeasible to address such a vast literature. As a starting point, we refer the reader to Morales et al (2014) where the authors study different market variants and strategies assuming a classical stochastic programming framework. We deal, in particular, with markets with a dualprice settlement for imbalances, under which there is no possibility of benefiting from a deviation and where imbalance penalties are asymmetric.
Next, we review relevant publications related to our research. Early works in this area proposed an optimal quantile strategy based on analytical results on the underlying newsvendor model and probabilistic forecasts of the wind energy production (Bremnes, 2004;Pinson et al, 2007). Several authors revisited this strategy. Zugno et al (2013) proposed constraining the offer to avoid sharp losses in single instances. Mazzi and Pinson (2016) employed a reinforcement learning algorithm to track the optimal quotient. Unfortunately, these methods involve several complex and independent components used to produce the probabilistic forecast and the quantile estimation.
Inspired by new advances in decision making under uncertainty in datarich environments, this problem regained interest in recent years. Stratigakos et al (2022) used an ensemble of decision trees that considers the objective function to estimate the energy production. From the modeling perspective, the work of Muñoz et al (2020) is one of the closest to ours, also aligned with the new stream of research that utilizes features to produce problem-specific decisions directly. In a similar problem, they leveraged recent advances in the newsvendor model, recently proposed by Ban and Rudin (2019), and proposed an approach that iteratively solves a linear optimization problem to update the model that produces the decision. Although relatively inexpensive, the time involved in its resolution may pose an issue in markets like the Australian NEM 1 , where 5-minute energy dispatches are carried out. Moreover, this approach seems redundant in the sense that the complete optimization problem is solved for each new step even though consecutive training sets may only differ by a few samples. These particular pitfalls pointed us to the thriving field of online learning in order to find suitable alternative approaches.
All the aforementioned references rely to some extent on stochastic programming, assuming that the uncertain wind production can be modeled by a fix probabilistic distribution or that the available historical samples (data) are independent and identically distributed (i.i.d.) variables. Online learning is a different decision-making paradigm that avoids these assumptions. Online learning studies problems where a recurrent decision is to be made and where the cost of each choice is available to feedback the next one. The target in this paradigm is to minimize the regret, which compares the performance of the algorithm against a benchmark, usually the best single action in hindsight. For an introduction to the subject, we refer the reader to the surveys of Shalev-Shwartz (2011) and Hazan et al (2016). Online Convex Optimization (OCO) focuses on the subset of online problems that face convex losses. Within OCO, we find a group of algorithms that continuously update the decision variables based on gradients or subgradients of a convex objective function, typically performing small but cumulative changes. These algorithms stand out for their versatility to deal with a wide range of problems and for their computational speed, provided that there is a closed-form expression to evaluate the (sub-)gradient (Duchi et al, 2011;Zheng, 2011). These methods offer long-term regret guarantees (Orabona, 2022) and have proven to be helpful in numerous applications on power systems (Gan and Low, 2016;Hauswirth et al, 2017;Colombino et al, 2019;Guo et al, 2021;Yuan et al, 2022). The wellknown online gradient descent, proposed by Zinkevich (2003), is the first in this category, achieving excellent results and giving rise to a wealth of derivative methods.
The strategy followed by the family of online gradient descent algorithms is in sharp contrast with most stochastic programming approaches, which solve an independent optimization problem with a different training set every time a decision has to be updated, e.g., the parameter of the decision rule in Muñoz et al (2020). Under convexity assumptions, an optimal solution can be found to each optimization problem, meaning that no single decision can ever achieve better performance on average in that training set. However, there is no certainty that the out-of-sample performance of such a decision enjoys the same privilege in finite sample sets. Instead, only probability guarantees can be offered even if the samples are i.i.d. (Van Parys et al, 2021).
Indeed, when samples are generated by a non-stationary environment following a process with seasonality or trend, the i.i.d. assumption does not hold and the out-of-sample performance can be very poor. This issue can be partly compensated using a rolling window setting (Bashir and Lehtonen, 2018) that updates the decision vector frequently. However, the dynamic process can evolve significantly within the training set. In that case, conventional stochastic programming approaches might still be hampered by unavoidably lengthy training set containing misleading old samples. On the contrary, online gradient descent algorithms update the decision through a point-wise update that involves the most recent information at the time, which enables capturing seasonality and trends even in short periods. Therefore, online gradient methods not only offer computational advantages but may also outperform their stochastic counterparts in terms of profits. In this line, the illustrative example provided in Section 4.2 and the case study presented in Section 5 corroborate these observations, showing faster tracking of non-stationary environments and an increase in the economic profit against rolling window approaches.

Contributions and structure
Far from being an isolated case, the Australian NEM is just an example of the existing trend pointing to a shortening in lead times and an increase in granularity, which is spreading across electricity markets. Without a doubt, these developments reduce the operational uncertainty, facilitating the integration of stochastic renewable energy sources 2 , although, at the same time, they raise the computational needs and require methodologies that adapt to changes fast. To face these new challenges, we propose an algorithm that combines a featuredriven newsvendor model inspired by Ban and Rudin (2019) with a variant of the online gradient descent algorithm presented in Zeiler (2012). Intending to explore the full potential of the proposed algorithm and in line with this trend, we conceive a case study in which we analyze an hourly forward market that closes just before the start of the next period. This demonstration allows us to use actual data of the Danish Transmission System Operator (TSO), Energinet 3 , and provide a relevant test bench to demonstrate the salient features that online learning algorithms can bring to power system problems. To the best of our knowledge, this is the first paper that analyzes the problem of trading wind energy in the online learning setting. Next, we formulate the contribution of our work which are threefold: • We develop an online offering algorithm for the wind power producer problem merging powerful tools from online learning, decision-making with contextual information and the technical literature of the wind power producer. Results show that this algorithm is computationally inexpensive and achieves substantial economic profits. • We propose a new non-stationary regret benchmark against which we empirically compare our algorithm. • We showcase the ability of the proposed algorithm to adapt to nonstationary scenarios through a concise illustrative example, and demonstrate the superior economic performance and computational efficiency of our approach on a case study that leverages more than five years of realistic data published by the Danish TSO, Energinet. The remaining of the manuscript is structured as follows: Section 2 introduces the problem of a wind power producer offering in the forward market with a dual-price mechanism. Section 3 develops a new offering algorithm based on an adaptive gradient descent algorithm and explores several performance metrics. Section 4 is built upon two illustrative examples that investigate the behavior of an alternative online implementation and the dynamic response of this algorithm in comparison with previous rolling window approaches. Section 5 empirically analyzes the performance of our proposed algorithm in a case study that considers an hourly forward market based on real data retrieved from the Danish TSO, Energinet. Finally, conclusions are duly drawn in Section 6.

Rolling window newsvendor
This section is devoted to introducing the problem of the wind power producer offering in a forward market, which is cleared sometime before their actual production is realized. Therefore, the producer is likely to suffer deviations from her offer that are cleared ex-post in a real-time market under a dualprice financial settlement for imbalances. Furthermore, the offer is assumed to be always accepted, as the marginal operational cost of wind farms is close to zero and therefore this technology is usually prioritized for being scheduled. The market revenue ρ ∈ R of a renewable energy producer in such a context is given by the summation of the amounts obtained in the forward (ρ F ) and in the real-time (also referred to as balancing) markets (ρ B ), i.e., where (a) + = max(a, 0). In addition, the unknown parameters λ F , λ UP , and λ DW ∈ R are the forward, up-regulation and down-regulation prices, respectively. The key decision variable for the renewable energy producer is her offer E F ∈ R + in the forward stage. Note that E ∈ R + denotes the actual realization of her stochastic energy production, which is not known in the forward stage.
In accordance with (1), the revenue in the forward stage (λ F E F ) is modified when the producer deviates from her offer E F . When the production is greater than expected E ≥ E F , the producer needs to buy downward regulation at price λ DW in the real-time stage to compensate the imbalance E − E F > 0. On the contrary, if she produces less than her forward offer E ≤ E F , the renewable energy producer has to buy upward regulation at price λ UP to compensate the imbalance E F − E > 0. Under a dual-price financial settlement, we have λ UP ≥ λ F and λ DW ≤ λ F , with at most one of them different from λ F (Morales et al, 2014, Ch. 7). In accordance with the aforementioned description, let ψ + , ψ − ∈ R + denote penalties for over-or under-production as Using (2) and (3) and the equivalence Note that the first term of (4) is out of the control of the underlying price-taker renewable producer (both λ F and E are uncertain parameters), the profit-maximizing offer can be computed through the minimization of the expected deviation cost, i.e., where E[.] is the expected operator. The optimization program (5) solves an instance of the very well-studied newsvendor model (Qin et al, 2011). Under a price-taker scenario, i.e., when the capacity of the producer is relatively small, an analytical solution to (5) can be computed through (Bremnes, 2004): where F −1 E is the cumulative distribution function (cdf) of the renewable energy production and the overline denotes the expected value of the random variable. The reader is referred to Maggioni et al (2019) for a discussion about the value of right distribution in newsvendor applications.
On top of the fact that the true distribution of the wind production and the optimal quotient are generally unknown, (6) suffers from another major drawback, which is its inability to directly leverage additional information that may be available, e.g., the wind energy forecasts of neighboring areas. In fact, it is usually the case that the renewable energy producer is equipped with a vector of auxiliary information or features x ⊆ X ∈ R p that may help explain the behavior of the uncertain parameters in (5). As proposed by Ban and Rudin (2019), this information can be exploited in newsvendor instances assuming that the optimal offer follows a linear decision rule of the form E F : X → R, E F = x ⊤ q with q ∈ R p being a new decision vector that determines the linear model. Then, considering that a set of historical samples of the form x t ), ∀t ∈ T in is available, we compute the best decision q LP for this set through the following linear program: where |·| denotes the cardinality of a set. Note that this model does not implicitly assume a price-taker scenario. In fact, correlations between penalties and wind features may be captured in systems with a high penetration of this technology. Although the linear structure of the mapping may seem restrictive, more complex relationships can be obtained by transforming the feature space, e.g., using a Taylor approximation (Ban and Rudin, 2019). Next, defining the box projection the optimal offer for a new piece of information x t ′ can be computed as Muñoz et al (2020), when new points are incorporated into the dataset T in , the problem (7) can be iteratively solved to update the value of q LP . In the remaining of the manuscript, we refer to this approach as LP (from Linear Programming).

Online trading
In the Online Convex Optimization (OCO) framework, a decision-maker faces an online problem where iterative decisions are to be made. The cost of each decision is determined by a convex loss function f t : R dz → R unknown beforehand. After a decision z t ∈ Z ⊆ R dz is made, the decision-maker learns f t and pays f t (z t ). Within OCO the Online Gradient Descent (OGD) algorithm, introduced by Zinkevich (2003), has proven to be very effective and versatile (Gan and Low, 2016;Narayanaswamy et al, 2012;Hauswirth et al, 2016;Nonhoff and Müller, 2020;Wood et al, 2021). Starting from an initial value, the OGD performs iterative updates z t based on (sub-)gradients of f t , denoted as g t from hereon. The magnitude of the step is controlled through a variable learning rate η t . On each round, the updated vector is forced to lie within the feasible region Z through the Euclidean projection. In OGD we rely on just the last point learned to obtain a gradient, thus resulting in a computationally inexpensive method, especially if the gradient and projection can be computed through a closed-form expression.
The selection of the learning rate is of paramount importance giving rise to countless variants of the gradient descent algorithm. The original proposal by Zinkevich (2003) presents two main alternatives, namely, a variable and a fixed learning rate. In a dynamic environment, the classical choice η t ∈ R + , η t ∝ t −1/2 is not suitable due to the fact that lim t→∞ = 0, reducing the ability to track changes as t increases. Alternatively, one could select a fixed value η t = η that keeps this capacity unaltered but may lose the fast convergence that the initial high values of η t provide. Regardless of the selection, both choices are scale-dependent and treat each component of the gradient vector equally. To tackle this, McMahan and Streeter (2010) and Duchi et al (2011) propose to use a component-wise adaptive rate η t ∈ R p and η t,p = η( t j=1 g 2 j,p ) −1/2 where g t = [g t,1 , ..., g t,p ] ⊤ . As in the case of η t ∝ t −1/2 , the previous expression is monotonically decreasing (component-wise), again limiting the long-term ability to learn. Aware of this limitation, Zeiler (2012) suggests an exponentially decaying average of the squared gradients to modulate the learning rate based on the most recent information. We leverage this gradient descent variant to implement our algorithm in Section 3.1.
In the online learning community, the de facto metric to evaluate the performance of a series of decision vectors z 1 , ..., z T is the regret R T ∈ R. The regret provides a versatile and, in a sense, normalized metric to compare an algorithm through different problems with the advantage that little assumption is made about the oracle that generates the decisions. Traditionally, the benchmark used to compute regret is the best single action in hindsight that can be obtained as the solution to an offline optimization problem under perfect information. However, in a dynamic environment, this benchmark can be beaten easily. In Section 3.3 we propose an alternative benchmark more suitable for the non-stationary context of the wind power producer problem.

Online newsvendor
In this section, we particularize the gradient descent introduced in the previous paragraphs to the context of the wind power producer offering in a forward market, incorporating elements of the rolling window problem presented in Section 2. We name the resulting algorithm OLNV (from OnLine NewsVendor). Contrary to the rolling window approach, the OLNV algorithm updates q based on the information provided by the last realization. The objective function (7) when the set T in reduces to one sample yields The OLNV method requires computing a gradient of the objective function, for which we analyze two alternative procedures in the following paragraphs.
The first approach is inspired by the work of Zheng (2011) on the pinball loss, a particular case of the objective function found in newsvendor models. Since the pinball loss is not strictly differentiable, the authors propose an alternative smooth approximation to ensure that computing gradients is always possible. Note that in our case the objective function (9) is not differentiable at E t = x ⊤ t q. Therefore, we first propose to circumvent this issue extending the approach in Zheng (2011) to the more general expression (9) that considers arbitrary (positive) penalties as where α > 0 is a parameter that controls the approximation and where higher values of this parameter result in smoother functions. The function N V t,α is convex in q and upper bounds N V t for any value of q as proved in Proposition 1 and 2 in Appendix A, respectively. Then, we derive a closed-form solution to obtain gradients of (10), yielding The second approach deals directly with the objective function as formulated in (9). Even though the original objective is not strictly differentiable, a variant of the OLNV algorithm is readily applicable to subdifferentiable functions, provided that a subgradient can be computed instead (Orabona, 2022). In this case, the mapping that returns a subdifferential of (9) is given by Note that, is a legitimate subgradient belonging to ∂N V t (q). For the sake of simplicity and reproducibility, the implementation of our algorithm returns zero whenever this condition is fulfilled.
Once a gradient as in (11) or a subgradient as in (12) has been computed, the key step of OLNV is to update q t using a multidimensional learning rate where • denotes the element-wise product, depending on the implementation of OLNV, and Π is a projection operator defined as Π : R p × X → R p . Precisely, Π maps its arguments into the solution of the following optimization problem: The feasible set in (14) is defined by the set-valued mapping Q : Note that, for any input x the output of Q is a convex region bounded by two parallel hyperplanes. As the Euclidean norm is used, a unique solution is guaranteed to exist for any instance of (14). Generally, the Euclidean projection of a point into a convex set requires solving a convex optimization problem, however the definition of Q allows us to find a closed-form expression, yielding This reduces the resolution of the optimization problem (14) to evaluating the above expression. Even though the operator Π guarantees the feasibility of q t under the realization x t−1 , we need to resort to (8) setting E F t = π(x t , q t ) to guarantee E F t remains feasible for any new arbitrary x t . The last remaining aspect is to compute the vector η t following the ideas in Zeiler (2012). Consider the running average of the squared gradient component is denoted as g 2 t,p where g t = [g t,1 , ..., g t,p ] ⊤ . Then, where ρ ∈ [0, 1) is a decay constant. Then, the learning rate applied to each component becomes where ǫ ∈ R + helps better condition the denominator and η > 0 is a constant. We use the update given by (16) and (17) in the proposed OLNV algorithm with the values ǫ = 10 −6 and ρ = 0.95, as originally suggested in Zeiler (2012). The benefits of this update is twofold. On the one hand, OLNV adapts the learning rate vector to the scale of every feature. On the other, OLNV is capable to track the most recent dynamic between the uncertain vector [E t , ψ + t , ψ − t ] and the feature vector x t . The complete OLNV algorithm particularized for the feature-driven wind power producer problem is compiled in Algorithm 1.

Regularization through average penalty anchoring
In a mature electricity market under a dual-price mechanism, it is common that ψ + t = ψ − t = 0 throughout a significant number of hours, meaning that load and generation are close to equilibrium in the system. From (9), it is evident that, in this situation, the producer experiences no cost no matter the deviation from the actual production. Moreover, the gradients computed through (9) are zero and therefore the variable vector q t is not updated, wasting information about the relationship between E t and x t . On a different front, when penalties are different from zero, they typically exhibit random behavior with sharp spikes representing highly imbalanced scenarios which, in turn, yields destabilizing updates of the vector q t . To tackle both issues, we propose performing the following convex transformation of the original penalties: where 0 ≤ µ ≤ 1 and ψ + , ψ − ∈ R + are the historical average penalties. This convex transformation is inspired by the concept of constraining the optimal offer around the point forecast developed in Zugno et al (2013), but unlike them, we do not impose hard constraints on the decision vector q t . Instead, we smooth the objective function using as anchor the sample average optimal market quantile determined by the average market penalties ψ + and ψ − . To do so, we consider a convex combination of the original objective function (7a) with an additional term that minimizes such a quantile, Then, by means of (18) and (19), the original objective structure is recovered, i.e., Therefore, by replacing ψ + t , ψ − t with ψ ′+ t , ψ ′− t in the original objective function, we regularize the learning procedure adding no extra computational cost. When the available samples are not sufficient to provide reliable estimators of the true ψ + and ψ − , the producer can resort to assume a balanced market with penalties ψ + = ψ − = 1. Thus, by the selection of µ < 1, provided that ψ + , ψ − > 0, the algorithm leverages the information contained in samples with both penalties equal to zero, potentially accelerating the convergence and obtaining smoother updates through the gradient. The same reasoning applies to the smooth objective function.

Performance evaluation
In order to assess the economic performance of an algorithm over a set of samples (E t , ψ − t , ψ + t , x t ), ∀t ∈ T oos , we can use the average deviation cost. To ease the notation, we set T = |T oos |. Consider that we are equipped with a vector of offers E F 1 , ..., E F T obtained from (7) and (8) or from Algorithm 1 after leaning one by one the samples belonging to T oos . We compute the average deviation cost as The value of this metric per se gives limited information about how a particular method is performing. A natural benchmark is the score obtained when a forecast of the wind energy production (in the sense of minimizing the root mean square error) is directly used as an offer in the market. We refer to this method as FO (from FOrecast). Let N V oos FO be the deviation cost incurred by FO. Leveraging this quantity, we redefine the original metric in relative terms, i.e., In this manner, the metric is expressed as a percentage of the possible improvement where a value of 100% means perfect performance with zero deviation cost.
On a different front, in online learning problems the customary performance measure is the regret. Traditionally, the regret compares a sequence of decision q 1 , ..., q T against the best single vector in hindsight q H . The latter is computed ex-post solving a problem analogous to (7) once the whole collection of samples belonging to T oos is known. Let Q H be the intersection of all feasible sets Given the assumption of a non-stationary environment, outperforming a constant q H can be a relatively easy task even though it is determined under perfect information. On the other side of the spectrum, one may consider the worst case regret (Besbes et al, 2015) interchanging the sum and minimum, i.e., where the second term of (25) computes the best individual decision q H t ∈ arg min q∈Q(xt) N V t (q). The regret computed in this way can be very pessimistic and unrealistic. Note that in the context of the wind power producer, it is always possible to find a value for q such that E t − x ⊤ t q = 0, and therefore (25) readily reduces to the summation of the original objective function R w T = t∈T oos N V t (q t ). Alternatively, Zinkevich (2003) proposed to compare the performance of online algorithms against a sequence of arbitrary decisions u 1 , ..., u T , u t ∈ Q(x t ), We refer to this approach as dynamic regret. This formulation allows to define a metric with an adjustable difficulty between the previous benchmarks. Note that (24) and (25) are special cases of (26) with u t = q H ∀t and u t = q H t ∀t, respectively. Then, the question is how to choose a reasonable series of comparators u t to use against OLNV. To this end, we propose dividing T oos in k adjacent partitions of equal length l, except possibly the last one. Without loss of generality, by assuming T − kl = 0, we have T oos Finally, the sequence of comparators that we propose to use in this paper is u t = q H i , ∀t ∈ T oos i . We will empirically investigate the regret performance of OLNV in the case study presented in Section 5.

Illustrative examples
This section analyzes several illustrative examples to gain insight into the behavior of OLNV. The first case compares the two alternative implementations introduced in Section 3.1 and discusses their main properties. As a result of this analysis, we select the subgradient objective function as the default procedure to perform the update of q t in OLNV. On a different front, one of the key features of online learning algorithms is their tracking ability, given the chronological order in which the updates are performed. In the second illustrative example, we deal with alternating penalty scenarios, demonstrating the salient properties of OLNV to adapt to a changing environment in a concise case.

Comparing the smooth and subgradient implementations
This illustrative example aims to elucidate whether the smooth approximation presented in (10) provides any advantage over the direct subgradient implementation of OLNV and determine the implementation to be used in the rest of the numerical experiments.
To this end, we consider a simplified setting with a single feature, a forecast of the actual wind production that we also use as the baseline of the method FO, and a single regressor variable q t ∈ R. No intercept is considered to ease the representation and analysis of q t . We sample the feature from a uniform distribution x t ∼ U (10, 90) (MW) and the true wind generation series is built adding a normal noise E t = x t + ǫ t with ǫ t ∼ N (0, 6) (MW). We produce a dataset of 1 year (8760 samples) through this process. Given that the penalties ψ + t and ψ − t are difficult to simulate, we compute them based on real dayahead and regulation prices of the Danish DK1 bidding zone. We retrieve data corresponding to the year 2017 from the data portal of the Danish TSO, Energinet 4 . Four implementations of Algorithm 1 are executed, three of them computing gradients of the smooth objective function through (11) with α = 0.05, 5 and 20 and the last one leveraging subgradients of the original cost mapping as in (12), to which we refer to as ∂. All instances are initialized with q 1 = 1, which means that the first offer produced by FO and OLNV are the same. In this section we do not use any convex transformation of the prices, i.e., µ = 1, and we set η = 0.005. We run the OLNV algorithm throughout the dataset, performing updates of q t every hour.
In this section, we accompany the numerical results with some theoretical analysis. The function N V t,α (q) approximate well the original function N V t (q) when |E t − x ⊤ t q| → ∞ as shown in Proposition 3 in Appendix A. Then, the question that remains is the behavior of both functions in a neighborhood of δ} with δ > 0 and δ close to zero. Let q 1 and q 2 be two vectors with E t − x ⊤ t q 1 ≤ 0, E t − x ⊤ t q 2 ≥ 0 and q 1 , q 2 ∈ ϕ. The subgradient that OLNV computes for each vector changes dramatically with ∂N V t (q 1 ) = ψ − t x t and N V t (q 2 ) = −ψ + t x t , which may result in diametrically opposed updates of the vector q for similar values of x t or q t . Conversely, N V t,α is everywhere differentiable, which ensures a smooth change of ∇N V t,α (q) for similar values of q − T and x t .   Table 1: Average absolute value |g| and standard deviation σ of the (sub-) gradients and the metric N V oos (%) computed for three smooth (α) and one subgradient (∂) implementations of the OLNV. Figure 1 shows a sample of ∂N V t and ∇N V t,20 that correspond to the subgradient and gradient of the smooth objective function with α = 20, respectively. Note that only N V t and N V t,20 are represented for the sake of clarity. It can be observed that most of the spikes in the case of N V t,20 are comparatively lower due to the aforementioned smoothing effect in the neighborhood of E t − x ⊤ t q = 0. This observation is aligned with the decreasing value of the standard deviation of the (sub-)gradients σ collated in Table 1 as α increases.
On the contrary, the mean absolute value of the (sub-)gradients, denoted as |g|, follows the opposite evolution. To understand the rationale behind this evolution, we provide Figure 2 showing three instances of the original and smooth losses. In all cases, we corroborate that N V t,α upper bounds N V t by a finite amount as proved in Proposition 2 in Appendix A. However, plot 2a evidences that the minimum of N V t,α is biased. This is true whenever ψ + t = ψ − t , a common situation in markets with a dual-price settlement for imbalances. Furthermore, when one penalty is equal to zero, the minimum is never attained.
In this line, the gradient computed through eq. (11) almost always introduces an error strictly greater than zero compared to the true value returned by eq. (12). The value of this error is given by the following expression: The imperfect approximation of N V t,α distort the magnitude and even the sign of the gradients, causing a long-term drift of q t that increases with the smoothing parameter α as shown in Figure 3. Finally, the last row of Table 1 presents the N V oos (%) obtained by each implementation with respect to FO. In this table, it is clear that N V oos (%) deteriorates when α increases. In summary, the smooth approach increasingly dampens the evolution of the decision vector for higher values of α but at the expense of a biased q t and important economic losses. Therefore, we conclude that the smooth approximation does not provide any substantial advantage over the subgradient implementation in this application, given that the producer is neutral to risk and volatility (only concerned with the economic return), and there is no technical constraint that encourage a smooth evolution of q. As a result of this analysis, we use subgradients to implement the OLNV method throughout the rest of the paper.

Dynamic behavior
In this illustrative example, we compare the tracking ability of OLNV and LP approaches in a non-stationary environment. Similar to the previous case, we assume that the producer has access to a unique feature and considers a model with a single regressor. Again, we sample the forecast from a uniform distribution x t ∼ U (10, 90) (MW) and the true wind generation series is built adding a normal noise E t = x t +ǫ t with ǫ t ∼ N (0, 6) (MW). Instead of the real DK1 data, we consider two possible scenarios with penalties ψ + t = 1, ψ − t = 3 and ψ + t = 3, ψ − t = 1, alternating every two months. This process generates 8 months of data (5760 hours) using the last 4 months (2880 hours test set. The start of the test set is aligned with the beginning of a twomonth scenario with ψ + t = 1 and ψ − t = 3. The rolling window approach is implemented solving the optimization problem (7) with a set of historical samples T in . Then, we leverage (8) to cast an offer based on the context E F t = π(x t , q LP t ). The coefficient q LP t is refreshed every 24 hours solving problem (7), following a rolling window. The reason for a 24-hour update is twofold: it is aligned with the original proposal in Muñoz et al (2020) and we empirically checked that there was little economic gain to be had with a lower update. The resolution time in the case of an hourly update, for example, took 24 times longer. As we later discuss, DR follows a rolling window approach that produces small changes in the training set, resulting in similar w DR t . We train four versions of the LP model with |T in | = 720, 1440, 2160 and 2880 (1, 2, 3, or 4 months), denoted as LP-1M to LP-4M, respectively. We use the first four months of the dataset to construct the initial training sets. Although the concept of training is substantially vaguer in OLNV, the last month of the training set is used to update the value of q t , originally initialized as q 1 = 1, to resemble a model that has been operating for some time. Figure 4 depicts the evolution of the single regressor q t in the test set, together with the optimal q * in each penalty scenario. In the first two months, the higher value of ψ − t penalizes offers above the true production E F t > E t and, consequently, the optimal strategy is to underestimate E F t with q * < 1. In the final months, the opposite is the case.
As one may expect, the evolution of the decision vector of LP models is smoother than in the case of OLNV, given that the first one considers many historical samples to perform the update. However, Figure 4 shows that the trajectory of q t produced by the rolling window models LP-1M to LP-4M is substantially lagged with respect to the change in the penalty scenario (emphasized by different background colors). This delay increases with the length of the training set to the point that LP-4M completely overlooks it.
Note that the length of the training set in LP-4M and the period of the penalty scenarios are identical. Therefore, the number of samples that penalize underor overproduction is equal and remains constant. As a result, LP-4M offers no incentive to overestimate or underestimate the forecast, yielding the same value as FO (neglecting slight deviations due to the finite sample and noise). Figure 4 also shows that OLNV is substantially faster tracking the optimal q * given that its updates are free from misleading past samples. In effect, the LP problem (7) determines the decision q t with the best performance on average in the training set, assuming that all the samples in the set are equally probable representations of future outcomes. Conversely, OLNV only uses the most recent information to perform a point-wise update that swiftly captures changes in the environment. The tracking capability of each model has an impact on its economic performance. Table 2 summarizes the out-of-sample N V oos (%) obtained by each approach in the test set. In line with the previous analysis, LP-4M obtains the same performance as FO. The other three LP methods experience decreasing N V oos (%) as the length of the training set and the lag of q t increase. Finally, the adaptability of OLNV clearly outperforms the LP approaches.
In this simplified example, we could have analyzed LP models with a shorter training set, probably resulting in reduced lag and better performances. However, in a realistic situation with a huge feature space and random penalties, months of data are typically required to capture the underlying relationships and generalize well in the out-of-sample set (Muñoz et al, 2020). Therefore, the length of the training set of the LP models has to be selected as a trade-off; enough data is required to learn a policy that generalizes well, but shorter sets capture dynamics better. On the contrary, the OLNV approach completely avoids this dichotomy, providing a fast and effective method that adapts to uncertain parameters generated by non-stationary environment.

Case study
Electricity markets are in the midst of a rapid development towards reducing the time between market transactions and the actual exchange of electricity. Examples of this transformation are given, i.e., by the reduction of the electricity lead time (Australian NEM or the Californian CAISO 5 ) or by the development of new intraday markets (OMIE intraday markets or NordPool : Evolution of q produced by five models over the test set. The blue and orange shaded periods correspond to the penalty scenarios ψ + t = 1, ψ − t = 3 and ψ + t = 3, ψ − t = 1, respectively. The entry q * corresponds to the best single vector for each penalty scenario. ELBAS 6 ). Inspired by this trend, we analyze a case study that considers an online forward market that takes place every hour followed by a balancing market with a dual-price settlement for imbalances. The closure of the forward market happens just before the start of the next period. We assume that the wind power producer constantly participates in the market and her offer is always accepted.
This section first describes the data utilized in this case study. Then, several benchmark methods are proposed to compare against OLNV. Finally, in the last part of this section, we analyze the numerical results obtained in the case study, including the regret, economic performance and computational cost of OLNV.

Data and experimental setup
This case study uses historical data compiled by the Danish TSO, Energinet.dk, since it includes market prices and several wind power forecasts that can be leveraged as quality features. We collect the true and day-ahead forecast issued by Energinet for the on-and offshore wind power production of both DK1 and DK2 Danish bidding zones together with the day-ahead and regulation prices of DK1 for the period 01/07/2015 to 06/04/2021 (mm/dd/yyyy). The day-ahead spot and regulation prices are mapped into hourly penalties through equations (2) and (3) and some small negative values, obtained due to rounding errors, are filtered out.
The raw wind power forecast series are also processed to be used in our case study. Given that the installed capacity of the four wind categories shown in Table 3 evolves differently over the dataset, we independently normalize each series to lie between 0 and 100 MW, a figure that can easily represent the capacity of a standard power plant. According to the Danish TSO, the raw wind power forecasts are issued between 12 to 36 hours ahead, although the exact time is difficult to know because no timestamp is provided. To overcome this issue, we use a standard ordinary least square regression model to produce enhanced forecasts with an accuracy comparable to an hour-ahead forecast and, therefore, suitable for our case study. We feed each raw wind power forecast into an independent linear regression model together with the last three lags of the true historical wind realization of the pertaining series. Finally, we use the first 6 months of our dataset to independently train each of the four predictive models, one per column of Table 3.   Table 4 compares the root mean square error (RMSE) of the original and improved out-of-sample forecast against the naive benchmark provided by the first lag of each series (the wind power production of the previous hour), also known in the literature as persistent estimator. Results show that the improved hour-ahead series significantly outperforms the original and persistent estimator. Hence, these enhanced forecasts are suitable for our case study. As a byproduct of this table, it is interesting to note that the wind forecasts issued by the Danish TSO have coherent RMSE being offshore harder to predict than onshore and DK2 harder than DK1 due to a reduced installed capacity and coverage area.
Once we have processed the wind power production series, we explain how we use them in our case study. The stochastic generation of the wind power producer offering in the market is simulated using the normalized onshore wind data series of the Danish DK1 bidding zone, which is consistent with the bidding zone of the imbalance penalties utilized. The four hour-ahead forecasts of the wind power production of DK1 and DK2 are available to the producer as contextual information. Although additional wind forecasts of neighboring bidding zones could have been used as features, we restrict ourselves to the ones produced by the Danish TSO to avoid potential inconsistencies regarding the issuing time that could cast doubt on the results obtained (Muñoz et al, 2020).
Given that our goal is to reduce the imbalance cost incurred by the wind power producer, we also consider several price-related features to be used as contextual information. To this end, we include the first lag of the imbalance penalties ψ + t−1 and ψ − t−1 in the vector of contextual information. As commented in Section 2, it is well known that the ratio between the penalties provides valuable information about the optimal decision of the newsvendor model and, therefore, we add the series r t−1 = ψ + t−1 /(ψ + t−1 + ψ − t−1 + υ) where υ = 10 −5 is a constant that helps better condition the denominator. Finally, we add a column of ones that enable one of the regressors to become an intercept to complete our feature set.
As a summary, let E on1 their findings, we choose a training set of |T in | = 4320 (6 months) and a rolling window step of 24 hours. In addition, we analyze a rolling window model, called LP, that solves exactly (7) and (8) using the full vector of available contextual information. This method is the one utilized in the illustrative example 4.2 but with different inputs. Given the similarities with the other rolling window approach LP2, we also choose a training set length of 6 months and a rolling window step of 24 hours.
Finally, we discuss a benchmark not implementable in practice inspired by the static regret metric defined in (24). We assume perfect information of the whole out-of-sample dataset and leverage problem (7) once more to compute the best linear model in hindsight determined by the vector q H . Once this optimal single vector is computed, the whole sequence of offers is determined through E F t = π(x t , q H ). We name this benchmark FX from FiXed. Next, we discuss the implementation of OLNV in this case study. The OLNV algorithm does not need to solve an optimization model but requires initializing two parameters. To choose µ and η, we perform a grid search on the chunk of data spanning 07/01/2015 to 12/31/2015. After analyzing Table 5, we select the values µ = 0.7 and η = 0.001 that obtain the best N V oos (%). We assume a balanced penalty anchor ψ + = ψ − = 1. Then, we initialize the OLNV regressor associated with the onshore DK1 forecast to 1 and the rest of the values to 0.01. In other words, we start the online offering with a strategy very close to FO, mainly relying on the forecast of the wind energy production. We use the next 6 months (01/01/2016 to 06/30/2016) to update (initialize) q OL with the aim of having a fair comparison against LP and LP2. The performance of all the methods presented in this section is evaluated using the dataset spanning from 07/01/2016 to 06/04/2021 (5 years with 43,200 samples). The optimization models LP, LP2, and FX are implemented with the Python package Pyomo (Bynum et al, 2021) and solved through the optimization solver CPLEX 7 , whereas the implementation OLNV is developed by the authors based on standard Python packages and uploaded to an open repository 8 .

Numerical results
Next, we discuss the results obtained in this case study. We start examining the regret suffered by OLNV over the aforementioned out-of-sample dataset with a length of D = 43, 200 hours (60 months). Let T oos We assess the average dynamic regret R d T /T for each sequence T oos j , j = 1, ..., D/l with partitions length l = 2160, 4320, 8640 hours (3, 6, 12 months). As an additional case, we compute the evolution of the static regret for a sequence T oos j , j = 1, ..., 20 with a step l = 2160 hours (3 months). To this end, in each step we refresh the best single action in hindsight as q H j = arg min q∈Q H j t∈T oos j N V t (q) and u t = q H j ∀t. The four aforementioned regret series are depicted in Figure 5. As expected, the average dynamic regret incurred by OLNV deteriorates quickly as l decreases since lower values of l translate in a more challenging benchmark closer to the the worst-case regret defined in (25). Nevertheless, Figure 5 clearly shows that OLNV achieves a sublinear static regret, i.e., lim T →∞ sup R s T /T ≤ 0. This is also the case for the dynamic regret with partitions of length l ≥ 6 months, proving the ability of OLNV to track dynamic environments.
The economic gains obtained by each method are assessed through the N V oos (%). The average values achieved over the evaluation dataset are collated in Table 6. First, note that all methods outperform the naive FO strategy of offering the DK1 forecast, obtaining positive values and demonstrating that this set of features contributes to reducing the deviation cost.
The LP2 method is developed in a context where recent lags in the penalties are not available. Indeed, the lack of penalty-related features translates into a modest score, showing the evident benefits of disclosing recent information in electricity markets, i.e., reducing the lead time. Even though FX determines the optimal q H in hindsight (i.e., under perfect information), its choice is limited to a single vector for the whole horizon. The fact that several approaches perform better than FX proves the dynamic behavior of the uncertain parameters and the need for updating the decision vector. Therefore, it does not come as a surprise that LP improves the first two approaches as it leverages the full vector of features and periodically updates q LP t . However, the superior adaptability of OLNV allows it to obtain the best score, achieving an additional 7.6% compared to LP and a total 38.6% deviation cost reduction compared to FO. The latter figure translates into an extra 25,930.22 e/year on average for a wind power producer with a capacity of 100 MW.
Finally, the last row of Table 6 summarizes the computational time corresponding to four approaches. The FX method requires little time as it only solves a single optimization problem for the whole horizon. This contrasts with the significant amount of time required by the constant re-optimization of LP and LP2. It is noteworthy that even though OLNV produces 24 times more updates of the vector q t , the time invested is several orders of magnitude lower. In conclusion, OLNV is up to the challenge of the electricity markets transformation achieving significant cost reduction together with exceptional computational performance.

Conclusions
This paper develops a new algorithm, named OLNV, combining a variant of the online gradient descent with recent advances that extends the newsvendor model to consider contextual information directly. The component-wise update of the learning rate enables the use of features with different scales seamlessly. In non-stationary environments, conventional stochastic approaches may consider misleading old samples in their training sets. On the contrary, our algorithm tracks the most recent information of the gradients, adapting the learning rate to follow the dynamics of the uncertain parameters and potentially obtaining higher profits. The closed-form expressions derived to compute the projection into the feasible region and a gradient of the objective function yield a computationally efficient algorithm that can be used in time-demanding problems. We envision the use of OLNV in future electricity markets that evolves toward continuous offering with reduced lead time. In particular, we apply this algorithm to the wind power producer problem offering in an hourly forward market with a dual-price settlement for imbalances. Several numerical experiments are carried out to assess the properties of the proposed OLNV algorithm. In the first illustrative example, we compare the behavior of two alternative implementations, namely, a subgradient approach and a smooth approximation of the original newsvendor function. The numerical and theoretical analysis provided in this example indicates that computing subgradient on the original objective function proves more profitable since it avoids update errors that may be introduced by the smooth approximation. Consequently, we determined that the subgradient implementation was the most suitable to this application and used it throughout the rest of the numerical experiments. Nevertheless, the smooth approximation could be utilized in other applications where other technical concerns advice a smooth update.
The second example shows the adaptability of the OLNV algorithm to nonstationary environments, clearly outperforming other stochastic approaches that optimize (using mathematical programming techniques) over a training set of past information. This superior performance is justified by the pointwise update that only uses the most recent information. Our case study, built upon real data of the Danish TSO Energinet, shows that OLNV achieves a 38.6% cost reduction against using a point forecast as offer and 7.6% compared to a state-of-the-art method. These significant improvements contribute to accelerating the integration of renewable energy technologies. Furthermore, we empirically analyze several dynamic definitions of regret, showing the desired sublinear convergence against most benchmarks.
Although this research focused on wind energy producers, OLNV is readily applicable to managing a portfolio of variable renewable energies with zero marginal cost, including wind, solar and other technologies. Similar algorithms can be developed when the producer's portfolio includes other assets such as loads, thermal power plants, or energy storage facilities, replacing the aggregated source of uncertainty, i.e., the variable net energy production, by a lineal decision rule. In this case, the projection step on the feasible region would likely involve solving a quadratic optimization program that can still be efficiently solved with modern solvers, provided that the feasible region is convex. Another attractive front is extending the OLNV algorithm to address inter-temporal constraints, observing a similar note with regard to the feasible region as in the previous case.
Future work also includes delving into the theoretical guarantees that this algorithm offers in terms of regret. On a different front, a wealth of other algorithms within the field of online learning can be applied to this problem, potentially bringing additional benefits such as faster convergence rates or improved performance. Similarly, variable selection techniques could help determine the subset of the available feature streams that provide the most economic value, whereas nonlinear mapping, i.e., kernels or generalized additive models (GAMs), can extend the regression capabilities of the method. Another exciting line of research concerns the risk analysis of the producer, where other metrics can be used instead of the expected value to create risk-averse strategies.