An Integrated Forecasting and Optimization Framework for Battery Energy Storage Participation in Spanish Electricity Markets

1  Introduction

The rapid growth of variable renewable energy has amplified short-term variability in net load, substantially raising the system value of flexible resources such as battery energy storage systems (BESS). Utility-scale BESS, standalone or collocated, can provide multiple services: energy arbitrage, imbalance mitigation, and participation in reserve markets. In liberalized European markets, these value streams are typically accessed through a sequence of market sessions (day-ahead, intraday, balancing and ancillary services), with progressively shorter lead times and different settlement rules. Capturing stacked value requires coordinated decisions across sessions and consistent treatment of operational constraints (state-of-charge dynamics, charge/discharge limits) and grid constraints (point-of-connection limits).

This paper aims to understand how multi-product electricity price forecasts, market-rule-constrained battery scheduling, robust uncertainty management, and chronological redispatch can be integrated into a single operational framework for PV+BESS participation in the Spanish electricity market. More specifically, the study investigates: (i) how forecast information from day-ahead, intraday, and secondary regulation products can be translated into feasible battery schedules; and (ii) how residual uncertainty, represented through Γ-budgeted uncertainty sets, affects bidding decisions, imbalance exposure, and state-of-charge management.

In the Iberian market (including Portugal and Spain), the day-ahead (DA) market clears quarter-hourly positions for the next day, followed by several intraday (ID) auctions and balancing mechanisms. Secondary regulation (SR) is procured to manage frequency deviations and is settled with availability and activation components. For an asset operator, these coupled mechanisms give rise to a multi-stage optimization problem in which decisions at each settlement layer depend on forecasts of prices, load, renewable generation, and reserve activation signals.

In practical operation, these decisions cannot be decomposed into independent forecasting, bidding, and redispatch tasks. A price forecast that is accurate on average may still lead to infeasible or economically suboptimal battery schedules if state-of-charge limits, point-of-connection constraints, reserve commitments, and previously cleared market positions are not considered simultaneously. Conversely, an optimization model that assumes perfect or static forecasts may overestimate the value of arbitrage and underestimate imbalance exposure. This coupling becomes especially relevant in multi-settlement markets, where each market gate closure fixes part of the schedule while later sessions only allow incremental corrections based on updated information. Therefore, the operational gap addressed in this work is the lack of a unified decision-support pipeline that transforms multi-product forecasts into market-rule-compliant PV+BESS schedules, updates those schedules sequentially as markets clear, and controls residual forecast uncertainty through a tractable robustness mechanism.

This paper contributes an end-to-end operational framework in which forecasting, deterministic scheduling, robust scheduling, and redispatch are not treated as independent modules but as coupled stages of the same bidding process. The proposed framework integrates four layers: (i) a multi-product forecasting stack for DA, intraday and SR prices; (ii) a deterministic mixed-integer linear programming (MILP) scheduler that enforces Spanish market coupling rules and two imbalance settlement representations; (iii) an adjustable robust optimization (ARO) extension to hedge against price and renewable uncertainty; and (iv) a redispatch simulator that sequentially updates the schedule as markets clear. While individual components like multi-product price forecasting, MILP self-scheduling, and robust optimization have been studied separately, the contribution of this work is their tight integration into a single operational pipeline specifically designed for the Spanish MIBEL market, featuring an empirically validated forecasting layer, explicit modelling of Spanish market rules and settlement logic using ESIOS-derived inputs, and a Γ-tunable robustness mechanism calibrated against commercial benchmarks.

Multi-market participation of storage has been addressed using bi-level formulations for price-taker batteries and market-clearing interactions, as well as single-level self-scheduling models for price-takers. Prior studies highlight the importance of including battery degradation and the coupling between energy and reserve commitments. Khalilisenobari and Wu propose a bi-level model for energy and ancillary service participation considering degradation costs [1], and Ghanaee et al. propose a MILP-based battery cycle-aging model for day-ahead energy and reserve scheduling in hybrid power plants, highlighting the relevance of degradation-aware scheduling for BESS participation in reserve markets [2]. In the context of European multi-settlement markets, multistage stochastic approaches have shown that adding intraday opportunities can increase expected profit and reduce imbalances for renewable-based virtual power plants [3].

On the one side, forecasting is a critical input to bidding and scheduling. In the literature, electricity price forecasting has been tackled with statistical time-series models (e.g., ARIMA), kernel methods (SVR) and modern machine learning and deep learning approaches. Gradient-boosted decision trees, particularly XGBoost, have proven competitive for non-linear relationships and heterogeneous features [4], while hybrid trend-seasonality models such as (Neural)Prophet capture multiple seasonality and holiday effects [5], and LSTM networks provide sequence modelling capacity for complex temporal patterns [6].

Moreover, uncertainty in prices and renewable production motivates stochastic and robust optimization. Robust optimization replaces random variables with uncertainty sets, producing solutions feasible for all realizations within the set. Ben-Tal et al. introduced adjustable robust counterparts for problems with recourse decisions, often approximated with affine decision rules for tractability [7]. Bertsimas and Sim proposed the budget-of-uncertainty concept, enabling a continuous trade-off between robustness and nominal performance while preserving linearity [8]. In energy systems, robust formulations have been applied to day-ahead participation problems under price and PV uncertainty [9] and to wind+BESS participation across energy and reserve products [10].

The remainder of the paper is organized as follows. Section 2 describes the data sources and feature engineering process. Section 3 presents the proposed methodology, including the forecasting layer, MILP scheduling formulation, and adjustable robust optimization extension. Section 4 reports out-of-sample forecasting results and illustrative scheduling outcomes. Section 5 discusses implications, limitations, and future directions, and Section 6 concludes.

2  Data Sources and Feature Engineering

The framework developed in this work is designed for operational deployment and therefore relies on publicly available or routinely accessible data. Market variables (day-ahead, intraday and ancillary service prices) are retrieved from the Spanish system operator’s (ESIOS) transparency platform via API [11]. Meteorological variables are obtained from the AEMET OpenData service [12] and fuel price fundamentals are represented by the daily day-ahead natural gas price index from MIBGAS [13].

All data series are aligned to the market time reference and aggregated to the corresponding product granularity. Missing values are handled by a combination of forward filling for slowly varying series (e.g., daily gas price) and interpolation for meteorological data. For operational consistency, the production pipeline re-trains models daily with fixed hyperparameters to include the features of the most recent available data.

For weather, AEMET provides historical data of daily minimum and maximum temperature and, for some locations, daily mean temperature, as well as the corresponding time of each value. To obtain an historical temperature profile usable for model training, we reconstruct a normalized intra-day temperature curve and then re-scale it to match each day’s observed minimum and maximum.

Fig. 1 illustrates the standard intra-day temperature profile used to reconstruct quarter-hourly temperature from daily statistics.

Figure 1: Average day temperature profile.

Calendar features include hour of day, day of week, month and day-type (working day vs. weekend/holiday). To capture structural breaks and extraordinary events, an event calendar is encoded as binary flags (e.g., COVID-related periods, geopolitical events, and other system-wide anomalies). These event features are appended to the feature vector for price models to improve generalization during atypical regimes.

More than a thousand energy-market indicators are provided by ESIOS. Given this large number, feature selection is performed using correlation screening and domain knowledge. Both linear correlation (Pearson) and non-linear dependence (maximal information coefficient, MIC) [14] are used to identify redundant signals. Features with Pearson and MIC correlation above a heuristic threshold (e.g., 0.7) are grouped, and representative indicators are chosen to minimize collinearity while preserving explanatory power.

PV generation forecasts for the hybrid plant can be derived from meteorological inputs using the pvlib Python library [15] and AEMET day-ahead sky-condition forecast.

3  Methods

The proposed methods are designed for a case study that represents a single PV+BESS market participant operating as a price-taker in the Spanish electricity market. The model does not simulate the full wholesale market-clearing process and does not represent individual competing generators. Instead, the rest of the power system is represented through exogenous market and system variables, including cleared market prices, renewable generation indicators, system balance variables, and reserve-related signals obtained from public data sources. The resulting tool is therefore best interpreted as a decision-support and self-scheduling framework for a price-taking PV+BESS operator, rather than as a market-clearing or strategic price-maker bidding model.

The modelled asset is a hybrid PV+BESS plant connected to the grid through a single point of connection. The point-of-connection constraint limits the net import and export of the aggregated PV and battery system. Network constraints beyond this point of connection are not modelled explicitly, which is consistent with the price-taker formulation adopted in this work. The scheduling model therefore focuses on the operational decisions available to the asset owner: charging, discharging, market energy sales and purchases, secondary regulation commitments, and imbalance management.

3.1 Forecasting Methodology

The forecasting layer produces short-term predictions for (i) PV generation (for operational conditions), (ii) day-ahead spot price, (iii) intraday auction prices (multiple sessions) and (iv) secondary regulation price signals. A portfolio of candidate models is evaluated offline, including support vector regression (SVR), autoregressive integrated moving average (ARIMA), gradient-boosted trees (XGBoost), NeuralProphet, and long short-term memory networks (LSTM). Model selection is performed offline via time-series cross-validation, using expanding windows to respect the temporal ordering of observations and prevent look-ahead bias. Hyperparameter tuning is performed once during the offline phase; the selected configuration is then frozen and the model is re-trained daily on a rolling window to incorporate the most recent market behavior. This design prioritizes operational stability over marginal in-sample gains, avoiding the risk of overfitting to recent regime shifts.

Performance is measured using absolute error (AE), root mean square error (RMSE) and relative error (RE). Where commercial benchmarks are available, we report head-to-head comparisons for identical evaluation windows.

The day-ahead price model uses 19 explanatory variables that combine temporal features with fundamental and market-based signals. In addition to calendar and weather variables, the model includes renewable supply proxies (solar and wind), a daily gas price index, and lagged values of multiple ESIOS price indicators from the previous day. Exogenous event flags are also included.

To interpret the trained models, we compute feature importance scores for the gradient-boosted trees. Fig. 2 reports a representative feature importance profile for the day-ahead spot model.

Figure 2: Most representative features for day ahead predictions with XGB.

“1_days_value”, “2_days_value”, and “3_days value” represent the price on the previous 1, 2 and 3 days at the same time, “2_week_avg” the average value of the previous 2 weeks on the same day of the week at the same time, “working”, indicates if it is a working day, and “renew” the amount of renewable energy generation expected for that time.

Intraday prices are forecast per session. The first session (S1) is scheduled shortly after the day-ahead market, and later sessions S2 and S3 occur closer to delivery. Each session model includes the latest available spot and intraday price information, updated renewable and demand signals, and the same calendar/event feature set.

Secondary regulation forecasts focus on the price of the reserve band (availability price), the related activation signals (up and down) and the final prices (up and down). Inputs include recent reserve prices, system balance indicators and calendar features. These forecasts are used as inputs to the scheduler when evaluating participation in the SR market and determining reserve capacity commitments.

The market session schedule and bid horizons used in the scheduling and redispatch workflow are summarized in Table 1.

3.2 Scheduling and Bidding Optimization

The optimization layer solves a self-scheduling problem for a PV+BESS asset participating as a price-taker in multiple markets. The model produces hourly/quarter-hourly bids for the day-ahead market (DA), adjustments for intraday auctions (ID1, ID2, ID3), and reserve commitments for secondary regulation (SR). The formulation is a mixed-integer linear program (MILP) designed for daily operation. The model determines optimal up- and down-bidding strategies across markets while ensuring physical feasibility, inter-market consistency, and compliance with operational constraints.

The BESS is modelled with charge and discharge power variables (ct,dt), binary indicators to prevent simultaneous charging and discharging, and a state-of-charge recursion. State of Charge (SOCt) is bounded between minimum (SOCmin = 0.1) and maximum (SOCmax = 0.9) limits. The point of connection (POC) enforces export (10 MW)/import (10 MW) limits and accounts for PV generation (PV peak = 20 MW), battery operation (max charge-discharge = 10 MW) and market transactions. Other relevant parameters are found in Table 2.

Let T denote the discretized planning horizon, and let M = {DA, ID1, ID2, ID3, SR} denote the set of markets. For each technology u ∈ U = {PV, BESS}, the decision variables pm,u,tdown≥0,pm,u,tup≥0 represent energy sold and purchased in market m at time t, respectively. Market quantities are expressed as energy over each scheduling interval. Battery charge and discharge variables are expressed in power units and converted into energy using the interval duration Δt. This convention ensures consistency between market commitments and physical battery operation. For the SR market, additional availability and activation variables are defined.

The objective function maximizes the total expected profit over all markets:

max(ΠE+ΠSR−ΠIB)(1)

being:

ΠE=∑t∈T∑u∈U∑m∈Mno−SRλm,t(pm,u,tup−pm,u,tdown),

ΠSR=∑t∈TλtSR(ptSR,up−ptSR,down)+∑t∈T(λtSR,act,upptSR,upπtSR,up−λtSR,act,downptSR,downπtSR,down),

ΠIB=∑t∈T(ϑIB+pIB,t++ϑIB−pIB,t−),

where:

•   λm,t: energy price at time t in market m

•   λtSR,act,up, λtSR,act,down: activation price for energy up/down in SR market at time t

•   ϑIB+, ϑIB−: positive and negative imbalance penalizations at t

•   pm,u,tup, pm,u,tdown: sold and purchased energy committed at time t, in market m for technology u, respectively

•   ptSR,up, ptSR,down: up/down energy committed in SR market

•   πtSR,up, πtSR,down: up/down activation ratio in SR market at time t

•   pIB,t+, pIB,t−: positive and negative imbalances at t

•   Mno−SR: set of markets without accounting for secondary reserve

Inter-market consistency is ensured through cumulative program variables:

progtDA=∑u∈U(−pDA,u,tup+pDA,u,tdown),(2)

progtID1=∑u∈U(−pID1,u,tup+pID1,u,tdown)+progtDA,(3)

progtID2=∑u∈U(−pID2,u,tup+pID2,u,tdown)+progtID1,(4)

and analogously for ID3 program. Upward bids correspond to energy injection into the system and may involve battery energy storage system (BESS) discharge. Conversely, downward bids indicate energy withdrawal from the system and may be associated with BESS charging.

The physical feasibility of aggregated trading actions is enforced through point-of-connection (PoC) limits and battery dynamics are modelled as:

0≤dt≤vtd,maxit,(5)

0≤ct≤Λtc,max(1−it),(6)

SOCt=SOCt−1+(γRTEct−dtγRTE)emax,(7)

SOCmin≤SOCt≤SOCmax,(8)

where:

•   dt, ct: battery charge and discharge power at time t

•   it: binary variable for discharge or charge state at time t

•   vtd,max, Λtc,max: maximum discharge and charge power at time t

•   SOCt, SOCmin, SOCmax: state of charge at time t, and min/max boundaries

•   γRTE: round trip efficiency of the BESS, symmetrically shared in charge/discharge

•   emax: maximum capacity of the BESS

The system-wide energy balance, including SR activation and imbalance variables, and PV production (θtS) is formulated as:

pIB,t+−pIB,t−=θtS+(∑m,upm,u,tdown+ptSR,downπtSR,down+dt)−(∑m,upm,u,tup+ptSR,upπtSR,up+ct)(9)

The optimization process is executed as follows: (i) The DA schedule defines a baseline energy position. Then, (ii) intraday auctions are modelled as incremental adjustments to the DA position, respecting the rule that intraday participation is conditional on having a prior DA position. (iii) Cleared offers from earlier sessions are fixed in later re-optimizations. (iv) A redispatch simulator executes this sequential logic, re-solving the MILP before each market gate closure using updated forecasts.

Two imbalance settlement representations are supported in the present approach. In the priced imbalance model, deviations between scheduled net injection and realized net injection are settled at an imbalance price (with separate prices for positive and negative deviations). In the penalized imbalance model, positive deviations are penalized at a fixed cost rate, representing a conservative approximation when imbalance prices are unfavorable or unknown. The choice affects risk appetite and is useful for sensitivity analysis. Both representations are consistent with the Spanish imbalance settlement rules under Royal Decree 1747/2003 and its subsequent amendments, which distinguish between deviations that help the system (settled at a more favorable price) and deviations that worsen system balance (settled at the imbalance price or penalized at a fixed rate).

Eq. (9) states the energy balance of the system, and links power variables ptDA and θtS with the charge and discharge of the BESS to define power imbalances at every hour of D+1. In other words, the imbalance of the system is given by the difference between the system generated energy and the system committed energy.

Secondary regulation is modelled as a reserve band commitment (up/down capacity) that must be sustainable given SOC limits and POC constraints. The objective includes an availability payment term (capacity price times committed reserve) and, optionally, an expected activation component based on historical activation ratios or scenario assumptions. Reserve commitments reduce the energy arbitrage flexibility and require maintaining SOC margins, creating an endogenous trade-off between arbitrage and ancillary service revenue.

3.3 Adjustable Robust Optimization Extension

To address uncertainty in market prices and renewable production, the deterministic framework is extended using Adjustable Robust Optimization (ARO) based on Γ-robustness principles. Uncertain parameters are modelled as belonging to bounded intervals around nominal forecasts.

For each uncertain price parameter λm,t, the uncertainty set is defined as:

λm,t∈[λ¯m,t−λ^m,t,λ¯m,t+λ^m,t],(10)

where λ¯m,t is the nominal forecast and λ^m,t denotes the maximum deviation (heuristically selected as an absolute or proportional value). A budget parameter Γm limits the number of simultaneous worst-case deviations.

Based on [8], the term of objective function related to the market m is:

∑tλ¯m,tpm,t−∑t(Γmzm,t+qm,t),(11)

subject to new constraints:

zm,t+qm,t≥λ^m,tpm,t,zm,t≥0,qm,t≥0.(12)

The ARO formulation augments the deterministic model with dual variables (zm,t, qm,t) for each uncertain parameter family, embedding risk control directly into the optimization problem. The variables zm,t and qm,t are auxiliary dual variables introduced by the Bertsimas–Sim robust counterpart to model market-price uncertainty. The variable zm,t represents the protection level associated with the uncertainty budget Γm for market m at time t, while qm,t captures the residual protection required when the adverse price exposure λ^m,tpm,t exceeds this protection level. Therefore, the term Γmzm,t+qm,t represents robust protection cost, or worst-case revenue loss, associated with price uncertainty at market m and time t.

The ARO counterpart is implemented as an optional robust scheduling mode and can be enabled for each uncertainty source: market prices, SR activations and solar production. Hence, the deterministic initial formulation detailed in Section 3.2 includes new terms and constraints aligned with enabled ARO versions. The resulting model remains a mixed-integer linear program with increased dimensionality but preserves tractability while offering structured protection against price volatility and renewable uncertainty.

The mathematical programming model is formulated in Pyomo framework for Python [16,17] which provides a flexible and abstract implementation. The solver used is HiGHS solver [18], given the speed of computation references.

4  Results

The framework is evaluated using Spanish market data from 2015 onwards for training, with 2023 selected for detailed out-of-sample testing. For selected periods, forecasts are compared against commercial data for benchmark. To reflect operational deployment, models are trained on data available up to each forecast date and evaluated on subsequent horizons. XGBoost achieved the best accuracy (lower RMSE) among the evaluated models.

Price benchmarks are reported for: (i) day-ahead spot prices (Table 3); (ii) the first intraday session S1 (Table 4); and (iii) the secondary regulation band price (Table 5). For intraday sessions affected by market redesigns, accuracy is reported for each session definition where sufficient historical data exists.

The proposed models achieve RMSE reductions of up to 31% relative to the commercial benchmark, with larger improvements during high-volatility periods. Improvements are more pronounced in high-volatility periods, suggesting that the feature set (including weather, renewables and gas) and non-linear model capacity provide robustness to changing fundamentals.

Fig. 3 provides an illustrative example of the day-ahead spot forecast against the real price series.

Figure 3: Market price prediction for day ahead market in 2023.

For the first intraday session (Table 4), the proposed model reduces RMSE by 47%–60% vs. the benchmark. This is operationally significant because S1 occurs shortly after the day-ahead market and provides the earliest opportunity to correct the day-ahead position as new information becomes available.

Secondary regulation band price results (Table 5) show modest RMSE improvements (5%–8%). The relatively low relative errors reflect the smaller dynamic range of this price signal compared to energy prices. Even small improvements may matter economically because reserve revenues depend on both capacity and the probability of activation.

Forecast accuracy directly affects the optimization layer by shaping expected arbitrage spreads and reserve revenues. Lower price forecast error reduces the risk of over-committing the battery in DA and then facing unfavorable corrections or imbalances. The ARO extension provides an additional control knob: increasing the uncertainty budget Γ shifts schedules toward more conservative positions, typically reducing cycling intensity and preserving SOC margins to accommodate adverse price or PV realizations. In operational practice, Γ can be adapted to forecast confidence (e.g., higher Γ during extreme weather or exceptional events) and to the operator’s risk tolerance.

To illustrate the end-to-end impact of forecast quality on scheduling decisions, consider the day-ahead case. A 30% reduction in RMSE translates directly into narrower uncertainty bands around expected arbitrage spreads. When spreads are better estimated, the MILP scheduler can commit the BESS more confidently to peak export periods, reducing the need for conservative SOC margins and unlocking additional revenue. The ARO layer, considered for the solar prediction, DA and SR, then acts as a second line of defense: by setting Γ > 0, the operator explicitly trades a fraction of expected profit for protection against the residual forecast error, particularly during high-volatility sessions where even the improved model may underperform.

Examples of the resulting market and battery schedule are shown in Figs. 4 and 5, including the price trajectory, PV production, BESS charge/discharge and the derived net energy plan. Fig. 4 shows a clear intraday structure in power positions, with net buying concentrated around midday (roughly hours 10–16), likely driven by higher solar generation, and more balanced or selling positions in early and late hours. Ancillary and imbalance components contribute consistently to downside exposure, particularly during peak activity, while day-ahead and intraday trades dominate the positive positions. Overall, the strategy appears more exposed during high-activity periods, with a mix of markets used to manage net positions across the day.

Figure 4: Hourly energy flows by market, prices, and PV production for the representative 24 h schedule. Positive values indicate buy bids, while negative values indicate sell bids.

Figure 5: Battery charge/discharge operation.

The sequential redispatch simulator is essential in multi-settlement markets because it respects the chronological information structure: only decisions for future sessions should be adjustable when new information arrives, while cleared transactions must remain fixed. This prevents unrealistic look-ahead bias and enables fair backtesting of integrated strategies. A more detailed image of the battery charge-discharge operation, together with the representation of the SOC can be observed in Fig. 5.

The optimized 24-h schedule shows that the BESS is operated as a flexible market asset rather than as a simple PV-shifting device. Instead of following a single daily pattern in which the battery charges only when PV generation is available and discharges later, the schedule combines several charge and discharge actions throughout the horizon. This behavior reflects the multi-market nature of the optimization problem: upward market positions can be supported by BESS discharge, whereas downward positions can be supported by BESS charging. As a result, the battery dispatch is shaped jointly by arbitrage opportunities, intraday corrections, secondary regulation value, point-of-connection limits, and imbalance exposure, rather than by PV availability alone.

The SOC trajectory in Fig. 5 confirms that these market decisions remain physically feasible, with the battery kept within its operational limits throughout the horizon. When secondary regulation participation is attractive, the scheduler must also preserve sufficient upward or downward flexibility so that reserve commitments remain deliverable. Therefore, part of the battery capacity is implicitly allocated to operational headroom or footroom, rather than being used only for energy arbitrage. This illustrates the value of the integrated framework for battery operators: it coordinates arbitrage, redispatch, reserve provision, and imbalance-risk management within a single feasible schedule.

5  Discussion

The results presented in Section 4 confirm that integrating multi-product price forecasting with a structured optimization layer yields measurable improvements over commercial benchmarks and produces schedules that are both physically feasible and market-rule compliant. This section elaborates on the practical implications of these findings, analyses the interactions between the framework components, and outlines directions for future research.

Although the proposed framework provides an integrated representation of forecasting, scheduling, robust optimization, and redispatch, its practical implementation also raises several operational challenges. The computational burden is particularly relevant in multi-market participation settings, especially when all market products are represented at 15-min resolution. In this configuration, both the deterministic MILP and the adaptive robust optimization formulation may become computationally demanding; in some cases, the HiGHS solver may require long solution times or fail to return a feasible solution within the operational time window. Therefore, practical deployment requires solver time limits, fallback deterministic schedules, warm-start strategies, and careful tuning of the temporal resolution and robustness budgets.

Potential model-coupling mismatches must also be considered. In the proposed framework, these mismatches are mitigated through backtesting validation, consistency checks, and parameter recalibration, which help align the behavior of the forecasting, optimization, and redispatch components. However, the robustness of this coupling has not yet been fully validated in a live production environment, where data delays, forecast drift, and market-design changes may affect performance. By contrast, data sharing is not expected to represent a major limitation, since the framework assumes a specified and validated data layer that provides each model with timely access to the required inputs for inference and decision recommendation.

5.1 Forecasting Performance

The day-ahead price model achieves RMSE reductions of up to 31% relative to the commercial benchmark, with the largest improvements concentrated in high-volatility periods. This pattern is consistent with the hypothesis that fundamental drivers, such as renewable supply proxies, gas prices, and weather variables carry predictive information precisely when markets are most stressed, i.e., when price spikes or negative prices occur due to excess renewable generation or unexpected demand drops. Standard benchmarks that rely primarily on lagged price signals tend to underperform in these regimes because they lack the causal structure needed to anticipate sharp deviations from recent history.

The first intraday session (S1) shows even larger gains (47%–60% RMSE reduction). This is operationally significant for two reasons. First, S1 occurs shortly after the DA market closes, providing the earliest opportunity to revise positions as renewable generation forecasts are updated. Second, S1 liquidity in MIBEL has historically been lower than DA liquidity, meaning that price discovery is noisier and a better model provides a more durable informational advantage. The ability to enter S1 with a more accurate price forecast allows the scheduler to capture correction trades that partially offset DA forecast errors before the system enters the penalized imbalance regime.

For the secondary regulation (SR) band price, improvements are modest in absolute terms (5%–8% RMSE reduction) but may be commercially meaningful. Reserve revenues depend not only on the band price but also on activation frequency: a small forecasting edge on the band price, combined with a better estimate of activation probability, can shift the expected SR revenue by several percent points, enough to affect the binary decision of whether to participate in the SR market at all. Future work could focus on jointly forecasting the band price and the activation probability as correlated outputs, potentially using a multi-task learning architecture that captures the structural dependence between these two signals.

It is also important to note that the evaluation windows used for benchmark comparison (Tables 2–4) span summer 2023, a period characterized by moderate renewable penetration and relatively stable gas prices following the 2022 energy crisis peak. Performance during more extreme periods, such as the negative-price episodes observed in spring 2024 due to excess solar generation, would provide a more demanding test of the model’s generalization capability. Expanding the evaluation to include such stress periods is a priority for future validation.

5.2 Optimization Layer

The MILP scheduler enforces Spanish market coupling rules through cumulative variables, ensuring that intraday adjustments are always incremental relative to previously committed positions. This design faithfully mirrors the actual settlement logic of MIBEL, where intraday bids are conditional on having a prior DA position and are settled as adjustments rather than independent quantities. Correctly modelling this coupling is non-trivial but essential: an uncoupled formulation would systematically overestimate the degrees of freedom available at each intraday gate, leading to over-optimistic backtesting results.

The two imbalance settlement representations (priced and penalized) serve complementary purposes. The priced model is appropriate for assets that can accurately forecast imbalance prices and whose deviations are sufficiently small to be settled at market rates. The penalized model is a conservative alternative that approximates worst-case imbalance costs, useful when imbalance price forecasts are unreliable or when the asset operator prefers a deterministic cost bound. In practice, operators may wish to switch between representations depending on market conditions: the penalized model provides stronger downside protection during high-volatility periods, while the priced model captures gains during stable regimes. A natural extension would be to endogenize this choice by modelling imbalance prices as an additional uncertain parameter within the ARO framework.

Secondary regulation modelling introduces an explicit trade-off between reserve revenue and energy arbitrage flexibility. Committing SR capacity reduces the effective SOC window available for energy trading, and the SOC margins required to sustain the committed reserve band limit the depth of charge/discharge cycles. The results show that this trade-off is sensitive to the assumed activation probability: if activation is underestimated, the asset commits too much reserve and loses arbitrage revenue; if activation is overestimated, the asset maintains unnecessarily large SOC margins. Improving the activation forecast, for example, by conditioning on real-time system balance indicators available from ESIOS is a key lever for improving SR participation decisions.

The ARO extension with budgeted uncertainty sets (Γ-robustness) provides a tractable mechanism for controlling conservativeness without solving a stochastic program with explicit scenarios. The key practical advantage is that Γ can be set by the operator as a single scalar per market product, allowing intuitive control over the risk-return trade-off: higher Γ produces more conservative positions (lower expected profit, lower downside risk), while Γ = 0 recovers the deterministic baseline.

However, the calibration of uncertainty sets remains an open challenge. In the current formulation, the maximum deviation λ^m,t is treated as a fixed parameter, typically estimated from historical forecast residuals. A more principled approach would derive λ^m,t directly from the predictive distribution of the forecasting model, for instance from quantile regression or conformal prediction intervals, creating a tight coupling between the forecasting and optimization layers. This would allow the uncertainty set to adapt dynamically to the current forecast confidence, tightening when the model is confident (e.g., stable weather, normal demand) and widening when uncertainty is high (e.g., storm events, system stress).

An important limitation of the current ARO formulation is that it handles each uncertainty source (market prices, SR activations, and solar production) independently. In practice, these quantities are correlated: low solar output on a cloudy day tends to coincide with reduced system-wide renewable production and higher energy prices. Ignoring this correlation leads to uncertainty sets that are either too pessimistic (when worst cases are assumed to occur simultaneously across all sources) or too optimistic (when correlations are not accounted for and individual budgets are set conservatively). Future work could explore ellipsoidal or factor-model-based uncertainty sets that capture these dependencies while preserving tractability.

Decision-focused (or task-based) learning, as recently explored for reserve market participation [19], offers a complementary direction: rather than minimizing forecast error, the forecasting model is trained to minimize the downstream scheduling regret, i.e., the gap between the profit achieved with the model’s predictions and the profit achievable with perfect foresight. Preliminary evidence in related settings suggests that task-trained models can reduce scheduling regret even when their point forecast accuracy is comparable to, or slightly worse, than standard models, because they implicitly learn to be accurate where accuracy matters most for the optimization.

5.3 The Redispatch

The sequential redispatch simulator is a critical component for realistic backtesting. By re-optimizing before each market gate closure using only information available at that time, the simulator enforces the chronological information structure that real operators must respect. This prevents look-ahead bias (a common source of overoptimistic results in energy system backtesting) and ensures that committed transactions from earlier sessions are fixed in subsequent re-optimizations, accurately reflecting the sequential and binding nature of electricity market commitments.

In operational deployment, several additional considerations arise. First, computational tractability is essential because the MILP must be solved within the time windows between gate closure and bid submission. Although many deterministic instances can be solved within practical time limits, the ARO augmentation and the use of 15-min resolution across multiple market products can substantially increase computational burden. Therefore, practical implementations should impose solver time limits, monitor infeasibility, warm-start from previous redispatch solutions when possible, and fall back to deterministic or conservative schedules if the robust formulation does not converge within the available time. Second, data pipeline reliability is critical: the framework depends on timely retrieval of ESIOS market data, AEMET meteorological forecasts, and MIBGAS gas price indices. Any latency or gaps in these feeds must be handled gracefully, with fallback strategies (e.g., persistence forecasts or conservative default positions) to ensure that the scheduler always produces a valid bid. Third, model drift monitoring should be embedded in the production pipeline: as market structures evolve (e.g., the ongoing MIBEL intraday market redesign), feature availability and price dynamics change, potentially degrading model performance. Continuous performance tracking with automated re-calibration triggers is recommended.

5.4 Considerations and Transferability

The framework is designed specifically for the Spanish MIBEL market but is architecturally transferable to other European multi-settlement markets. The key market-specific elements are: (i) the session schedule and gate closure times encoded in Table 1; (ii) the imbalance settlement rules; (iii) the SR product definition (availability band plus activation); and (iv) the data sources (ESIOS, AEMET, MIBGAS). Adapting the framework to, say, the French EPEX market or the Italian GME would require updating these elements and retraining the forecasting models on the corresponding data, but the overall architecture (forecasting layer, MILP scheduler, ARO extension, redispatch simulator) would remain intact.

An important consideration for transferability is the degree of market coupling between countries. In the Iberian Peninsula, Portugal and Spain share a common day-ahead price through market coupling with the rest of Europe via EUPHEMIA. Cross-border flows and interconnection constraints can induce price differences (Iberian spread) that represent an additional source of uncertainty for asset operators. Explicitly modelling this spread, and its interaction with renewable generation on both sides of the border, could be a valuable extension for operators with cross-border assets or import/export positions.

From a regulatory perspective, the ongoing implementation of the EU Clean Energy Package (Directive 2019/944) is progressively harmonizing imbalance settlement rules, reserve product definitions, and market time unit granularity across European member states. This convergence will likely reduce the market-specific engineering effort required to adapt the framework and increase the relevance of transferable approaches such as the one presented here.

6  Conclusions

This paper presented an end-to-end decision-support framework for PV+BESS participation in the Spanish electricity market, integrating multi-product price forecasting, deterministic MILP scheduling, adjustable robust optimization, and a sequential redispatch simulator. The framework is designed for daily operational deployment using publicly available data from ESIOS, AEMET, and MIBGAS, and is validated against a commercial forecasting benchmark on 2023 out-of-sample data.

The key findings are as follows. On the forecasting side, XGBoost models combining fundamental drivers (renewable proxies, gas prices, weather) with calendar and event features achieved consistent improvements over the benchmark: up to 31% RMSE reduction for day-ahead prices, 47%–60% for the first intraday session, and 5%–8% for the secondary regulation band price. The largest gains occur in high-volatility regimes, where fundamental-driven models outperform lag-based benchmarks by capturing the causal structure of price spikes and negative-price episodes.

On the optimization side, the MILP scheduler correctly models Spanish market coupling through cumulative variables and supports two imbalance settlement representations, enabling sensitivity analysis across different risk profiles. The ARO extension with Γ-budgeted uncertainty sets provides a tractable, operationally intuitive mechanism for trading expected profit against downside protection: higher robustness budgets produce more conservative SOC trajectories and reduce cycling intensity, limiting imbalance exposure during periods of elevated forecast uncertainty. The redispatch simulator enforces chronological feasibility across market sessions, preventing look-ahead bias and enabling fair backtesting of integrated bidding strategies.

From a practical standpoint, the results support several operational recommendations. First, multi-product forecasting covering DA, intraday, and SR prices jointly is more valuable than single-product forecasting because it allows the optimizer to evaluate cross-market trade-offs consistently. Second, the robustness budget Γ should be adapted dynamically to forecast confidence: a higher Γ is warranted during extreme weather events, market redesign transitions, or other atypical regimes where model generalization is less certain. Third, the penalized imbalance representation provides a useful conservative baseline for risk-averse operators or during periods when imbalance price forecasts are unreliable.

Looking ahead, three extensions stand out as high-priority. First, tighter coupling between the forecasting and optimization layers through decision-focused learning could potentially reduce scheduling regret by training price models to be accurate where it matters most for the bidding problem, rather than minimizing average forecast error. Since decision-focused learning is not implemented in the current framework, a systematic comparison with the XGBoost-based forecasting pipeline is left for future work. Second, probabilistic uncertainty sets derived from quantile regression or conformal prediction would allow the ARO layer to adapt automatically to current forecast confidence, replacing the static deviation parameters used in the current formulation. In addition, future versions of the framework should explicitly account for the correlation between price, PV-production, and reserve-activation uncertainty. Treating these uncertainty sources independently may lead to overly conservative or overly optimistic aggregate uncertainty representations, which limits the robustness of the current framework for production-level deployment. Third, incorporating battery degradation costs into the objective function would provide a more complete economic picture, particularly for assets with limited cycle life where aggressive arbitrage strategies may erode long-term value.

The modular design of the framework, with clearly separated forecasting, scheduling, and robustness layers, facilitates both incremental improvements to individual components and full adaptation to other European multi-settlement markets. As renewable penetration continues to grow and electricity price volatility intensifies, integrated forecasting-optimization pipelines of this type will become increasingly essential for the profitable and reliable operation of flexible storage assets.

Acknowledgement: This work has been supported by the Agencia Estatal de Investigación of the Spanish Ministerio de Ciencia, Innovación y Universidades under grant CPP2021-008632. Anaïs Pepey is acknowledged for providing support in BESS technical representations. Jose María Plana, Ben Hill, Joaquim Palau, Regina Enrich and Eduardo Aguilar are acknowledged for their continuous support and understanding of electricity markets, power systems and business insights.

Funding Statement: This work was supported by the public-private collaboration project DIOPTIBA (CPP2021-008632), funded by MCIN/AEI/10.13039/501100011033.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Pol Torres; methodology, Pol Torres, Marcelus Fabri, Lluís Millet and Diego Gallego; predictive models, Pol Torres, Alejandro Clavera, Diego Gallego and Pau Plana; Optimization models, Marcelus Fabri, Giulia Zarpellon, Carlos Molina, Lluís Millet and Pau Plana; software integration, Diego Gallego; writing of original draft, Pol Torres; writing, review and editing, Pol Torres, Carlos Molina, Diego Gallego. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data not available due to commercial restrictions. Due to the nature of this research, participants of this study did not agree for their data to be shared publicly, so supporting data, materials, and code is not available.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Khalilisenobari R, Wu M. Optimal participation of price-maker battery energy storage systems in energy and ancillary services markets considering degradation cost. Int J Electr Power Energy Syst. 2022;138:107924. doi:10.1016/j.ijepes.2021.107924. [Google Scholar] [CrossRef]

2. Ghanaee E, Pérez-Díaz JI, Fernández-Muñoz D, Nájera J, Blanco M. Accurate MILP-based battery cycle aging model for day-ahead energy and reserve scheduling in hybrid power plants. J Energy Storage. 2026;155:121426. doi:10.1016/j.est.2026.121426. [Google Scholar] [CrossRef]

3. Silva AR, Pousinho HMI, Estanqueiro A. Energy market participation of flexible resources with a multi-stage stochastic approach. Energy. 2022;258:124856. [Google Scholar]

4. Chen T, Guestrin C. XGBoost: a scalable tree boosting system. arXiv:1603.02754. 2016. [Google Scholar]

5. Triebe O, Hewamalage H, Pilyugina P, Laptev N, Bergmeir C, Rajagopal R. NeuralProphet: explainable forecasting at scale. arXiv:2111.15397. 2021. [Google Scholar]

6. Hochreiter S, Schmidhuber J. Long short-term memory. Neural Comput. 1997;9(8):1735–80. doi:10.1162/neco.1997.9.8.1735. [Google Scholar] [PubMed] [CrossRef]

7. Ben-Tal A, Goryashko A, Guslitzer E, Nemirovski A. Adjustable robust solutions of uncertain linear programs. Math Program. 2004;99(2):351–76. doi:10.1007/s10107-003-0454-y. [Google Scholar] [CrossRef]

8. Bertsimas D, Sim M. The price of robustness. Oper Res. 2004;52(1):35–53. doi:10.1287/opre.1030.0065. [Google Scholar] [CrossRef]

9. Correa-Florez CA, Michiorri A, Kariniotakis G. Robust optimization for day-ahead market participation of smart-home aggregators. Appl Energy. 2018;229:433–45. doi:10.1016/j.apenergy.2018.07.120. [Google Scholar] [CrossRef]

10. Fernández-Muñoz D, Pérez-Díaz JI. Optimisation models for the day-ahead energy and reserve self-scheduling of a hybrid wind-battery virtual power plant. J Energy Storage. 2023;57:106296. doi:10.1016/j.est.2022.106296. [Google Scholar] [CrossRef]

11. Red Eléctrica de España. ESIOS API documentation. [cited 2026 Jan 1]. Available from: https://www.esios.ree.es/es/pagina/api. [Google Scholar]

12. Agencia Estatal de Meteorología (AEMET). AEMET OpenData. [cited 2026 Jan 1]. Available from: https://opendata.aemet.es/centrodedescargas/inicio. [Google Scholar]

13. MIBGAS. Daily indices (PVB day-ahead). [cited 2026 Jan 1]. Available from: https://www.mibgas.es/en/market-results/daily-indices. [Google Scholar]

14. Zhou P, Zhang Z, Yan Y, Zhao S. Unknown type streaming feature selection via maximal information coefficient. In: Proceedings of the 2022 IEEE International Conference of Data Mining Workshops (ICDMW); 2022 Nov 28–Dec 1; Orlando, FL, USA. p. 650–7. doi:10.1109/ICDMW58026.2022.00089. [Google Scholar] [CrossRef]

15. Holmgren WF, Hansen CW, Mikofski MA. Pvlib Python: a Python package for modeling solar energy systems. J Open Source Softw. 2018;3(29):884. doi:10.21105/joss.00884. [Google Scholar] [CrossRef]

16. Bynum ML, Hackebeil GA, Hart WE, Laird CD, Nicholson BL, Siirola JD, et al. Pyomo—optimization modeling in Python. 3rd ed. Berlin/Heidelberg, Germany: Springer; 2021. [Google Scholar]

17. Hart WE, Watson JP, Woodruff DL. Pyomo: modeling and solving mathematical programs in Python. Math Program Comput. 2011;3(3):219–60. doi:10.1007/s12532-011-0026-8. [Google Scholar] [CrossRef]

18. Huangfu Q, Hall JAJ. Parallelizing the dual revised simplex method. Math Prog Comp. 2018;10(1):119–42. doi:10.1007/s12532-017-0130-5. [Google Scholar] [CrossRef]

19. Terada LZ, Magalhães MM, Cortez JC, Soares J, Vale Z, Rider MJ. Hybrid robust and decision-focused learning for secondary reserve markets. Sustain Energy Grids Netw. 2025;42:101773. [Google Scholar]