Optimal Spot Market Participation of PV + BESS: Impact of BESS Sizing in Utility-Scale and Distributed Configurations

1. Introduction

1.1. Motivation

In recent decades, policymakers have intensified efforts to reduce carbon emissions, leading to significant transformations particularly in the energy production sector, which accounts for more than three-quarters of total global greenhouse gas emissions [1]. A key outcome of these efforts has been the widespread deployment of Renewable Energy Sources (RES) worldwide, gradually replacing conventional thermal power plants. This transition has been driven by substantial public and private investments, as well as by international agreements aimed at facilitating the energy transition such as the Paris Agreement [2] and the European Green Deal [3]. Between 2015 and 2024, the effects of these global commitments became evident, as the share of electricity generated from RES surpassed 30% in 2023, a notable increase from 23% in 2015 [4]. This progress has been primarily driven by the expansion of solar photovoltaic (PV) and wind energy, with PV alone accounting for 11% of total electricity production in Europe in 2024 [5].

This ongoing transition is also driving the decentralization of electricity generation, with increasing numbers of small- and medium-scale production units connected to distribution grids. This trend is particularly pronounced in the European solar power sector, where distributed installations accounted for 58% of new PV capacity added in 2024 [6]. Consequently, decentralized energy systems are gaining prominence in European regulatory frameworks. In Italy, the integration of Distributed Energy Resources (DERs) into electricity markets has been formalized through the Testo Integrato del Dispacciamento Elettrico (TIDE) [7], which establishes market-based mechanisms to compensate DERs similarly to utility-scale generation units. This national framework aligns with broader European legislation, notably Directive 2019/944 [8], which sets out specific guidelines for DERs participation in electricity markets. Central to these regulations is the role of aggregators, who coordinate and manage multiple DERs within Virtual Power Plants (VPPs), enabling their efficient and effective market participation.

However, the increasing penetration of non-programmable renewable resources (NP-RES) introduces new challenges for grid stability. The inherent uncertainty associated with these energy sources complicates the real-time balance between electricity generation and consumption. Additionally, the low inertia of NP-RES reduces the system’s ability to counteract frequency fluctuations, thereby impacting overall grid reliability. As a result of the increased demand for real-time balancing energy and of the greater difficulty in finding flexible resources, the costs of ancillary services are expected to increase in future years [9]. These services, purchased by the Transmission System Operator (TSO), are essential for maintaining voltage and frequency stability and preventing grid disturbances or blackouts [10].

One of the most effective solutions to address these challenges is the integration of Battery Energy Storage Systems (BESS) with NP-RES, enhancing system flexibility by decoupling energy generation from consumption throughout the day. BESS can be installed in two main configurations: front-of-the-meter (FTM) and behind-the-meter (BTM) [11]. FTM batteries, typically large-scale systems connected directly to transmission or distribution networks, contribute to grid stability by enabling energy shifting, supporting peak demand reduction, and facilitating renewable energy integration [12]. These systems also participate in wholesale markets and deliver critical ancillary services such as frequency regulation, voltage control, and black-start capabilities [13]. With advanced inverters, they can enhance system inertia and fault response, helping stabilize increasingly decarbonized power systems [14]. In contrast, BTM batteries are deployed at residential, commercial, or industrial sites and primarily serve end-users [12]. They help optimize self-consumption of onsite renewable generation, reduce electricity bills via peak shaving and time-of-use arbitrage, and ensure supply reliability during outages [15]. Additionally, BTM systems can support local grid operations by providing flexibility services to Distribution System Operators (DSOs), easing network constraints, and deferring infrastructure upgrades [12].

Given the central role that PV systems play in energy transition—both as utility-scale assets and as distributed generation aggregated through VPPs—it becomes increasingly important to assess the economic and operational value of integrating BESS with PV installations. Such coupling not only has the potential to improve the financial viability of PV projects, but also helps reduce output uncertainty and enhance the reliability of the electrical system.

1.2. Literature Review and Research Gaps

Several studies focus on the development of optimization models for the participation of PV systems in electricity markets, aiming to enhance their integration into the power system. Due to the inherent variability of solar irradiance throughout the day, scenarios generation to represent possible PV production profiles is a crucial starting point in many of these works. Scientific literature presents a wide range of methodologies for this purpose, with Monte Carlo methods being among the most commonly employed. Zheng et al. [16] explore various approaches for modelling the power output of PV systems under diverse meteorological conditions, emphasizing that Monte Carlo methods are particularly effective in generating representative probabilistic production scenarios for DERs. To mitigate the computational complexity of handling large scenario sets while preserving key statistical attributes, Shi et al. [17] propose an approach that integrates a reduced Monte Carlo method with k-means clustering. Similarly, Falabretti et al. [18] apply this combined methodology to construct input scenarios for a stochastic mixed-integer linear programming (MILP) model, designed to optimize the bidding strategy of a VPP in both the Day-Ahead Market (DAM) and the real-time market. By explicitly accounting for PV production uncertainty, their study demonstrates significant improvements in profitability while effectively mitigating imbalance costs.

Given the uncertainty in PV power output, other studies apply stochastic or robust optimization models to maximize revenues from PV participation in electricity markets, demonstrating superior performances over traditional deterministic approaches. Stochastic optimization aims to find solutions that perform well on average by incorporating the probability distributions of uncertain parameters. In contrast, robust optimization seeks solutions that remain effective under worst-case scenarios within a defined uncertainty set, without requiring probabilistic information [19]. For instance, Nemati et al. [20] propose a regret-based flexible robust optimization model for a VPP composed of wind, solar PV, and flexible demand, simulated in a southern Spain scenario. The VPP participates simultaneously in the day-ahead and secondary reserve markets, while accounting for imbalance penalties. Their model allows the operator to control the level of conservatism using a monetary regret limit rather than traditional uncertainty budgets. Results show that the model is economically effective and that it can adapt to different levels of forecast error and market penalties. In contrast Silva et al. [21] develop a multistage stochastic decision-aid algorithm to optimize the bidding strategy of a VPP composed of wind, solar PV, and storage across day-ahead, intraday, and balancing markets. Unlike previous studies focusing mainly on day-ahead and balancing markets, their approach explicitly includes intraday market participation to better manage forecast uncertainty. Results show that engaging in all three markets increases VPP profits by 10.1% and reduces imbalances by 63.8%, highlighting the value of coordinated multi-market strategies under uncertainty and the relevance of near-real-time re-scheduling possibility for NP-RES. Similarly, Visser et al. [22] propose an operational bidding strategy that optimizes the participation of a PV power plant in electricity spot markets, considering day-ahead, intraday, and imbalance markets. Particularly, they develop a multistage stochastic optimization method using a scenario generation algorithm. The method is evaluated using both technical and economic metrics, demonstrating a significant improvement over a reference strategy. The study shows that extending market participation to intraday increases revenues by 22% and reduces imbalance by 50%.

While the optimization of utility-scale PV + BESS systems for participation in electricity markets has been widely studied, the optimization of residential PV + BESS systems also offers significant benefits, particularly in enhancing energy flexibility and improving economic performance [23]. Li [24] develops a genetic algorithm-based method for the optimal sizing of grid-connected PV + BESS systems in residential buildings, incorporating a time-of-use electricity tariff scheme to maximize cost savings. Results show that jointly optimizing PV and battery sizes significantly reduces electricity imports and costs, though optimal system size varies with household consumption profiles and price conditions. Duman et al. [25] propose an optimal PV + BESS sizing model for prosumers equipped with Home Energy Management Systems (HEMS), emphasizing the role of day-ahead load scheduling in increasing self-consumption. The developed MILP model incorporates optimal PV tilt angle selection and load scheduling of various controllable appliances, while accounting for battery degradation and vehicle-to-home functionality. It first minimizes daily electricity costs through demand response and self-consumption, then simulates one year of HEMS operation to evaluate the Net Present Value (NPV) of different PV + BESS configurations. The optimal design is identified as the configuration yielding the highest NPV. Similarly, Zhou et al. [26] explore the allocation of PV and battery capacity in a smart home environment, explicitly considering the impact of electricity pricing mechanisms on system performance and cost-effectiveness. Particularly, the study accounts for various electricity pricing mechanisms—time-of-use (TOU), real-time pricing (RTP), and stepwise power tariffs (SPT)—as well as PV subsidies and uncertainties in solar generation and seasonal load profiles. A hybrid optimization approach, combining a cataclysmic genetic algorithm with the DICOPT solver in GAMS [27], is employed to identify the optimal configuration. Through six case studies and sensitivity analyses, the results reveal that under SPT and sufficient subsidies, installing only PV is optimal, whereas under RTP and TOU, both PV and BESS are necessary to achieve maximum profits. Rezaeimozafar et al. [28] propose a two-stage stochastic optimization framework for operating behind-the-meter PV + BESS systems, aiming to minimize daily grid consumption and battery degradation while accounting for PV output uncertainty. A feed-forward neural network combined with a statistical error analysis technique is used to generate accurate PV output scenarios, which are then reduced using a backward scenario reduction method to mitigate dimensionality issues in the optimization model. The reduced scenarios are input into the two-stage model to estimate expected electricity costs and degradation. The approach is validated with real-world data from a household in Ireland, demonstrating strong performance under PV variability and results comparable to an ideal forecast scenario.

A key aspect when optimizing commercial PV programs is to assess the economic impact of imbalances. Guo et al. [29] propose an optimal real-time operational framework aimed at minimizing imbalance costs of a VPP. The study highlights that the effectiveness of the adopted strategy depends on the specific imbalance settlement mechanisms in place. Similarly, Marneris et al. [30] investigate the rationale behind strategic bidding and the necessity of accurate forecasting under different imbalance settlement schemes. Their study presents a comprehensive approach to determine the optimal bidding strategy for a RES aggregator operating in the day-ahead and ancillary services markets, where prices are taken as given. To ensure precision in the optimization process, the interactions between the day-ahead market, ancillary services, and real-time balancing mechanisms are explicitly modelled.

Imbalance costs, and price signals coming from different energy markets, should correctly reflect the impact of PV-induced imbalances on the working operations of the power system managed by the TSO. Pierro et al. present a structured analysis addressing the challenges of integrating high shares of PV generation into the Italian power system [31,32,33]. In their earliest work [31], they propose two complementary strategies to mitigate the impacts of rising PV penetration on imbalance volumes: improving PV and net load day-ahead forecasts, and converting conventional PV systems into “flexible PV” plants equipped with remote curtailment capabilities and cost-optimized battery storage. Both strategies are shown to significantly reduce imbalance volumes and costs, with the second capable of eliminating imbalance impacts entirely, all at a cost lower than current imbalance management practices. Building on this foundation, Pierro et al. in [32] outline a broader roadmap to achieve a 100% renewable electricity mix in Italy by transforming intermittent PV production to firm 24/365 power generation. In their most recent contribution [33], the authors focus on the system-wide implementation of flexible PV plants to provide ancillary services through proactive curtailment and battery support. They show that such plants can reduce the national imbalance by 36% compared to 2016 levels with only 6% curtailment of PV production. This can be achieved at or below current dispatching costs. Moreover, a geographically distributed fleet of flexible PV plants ensures optimal regulation performance, with the sizing strategy found to be robust against year-specific variations in load, generation, and market conditions.

Complementing this system-level perspective, Lazard offers a widely referenced benchmark for the cost of transforming PV generation into firm power through battery storage [34]. According to their last report, the Levelized Cost of Electricity (LCOE) for standalone PV in the CAISO market [35] is approximately $43/MWh. When battery-based firming is included, this cost rises to $141/MWh. The firming cost is defined as the additional expense required to compensate for the non-firm portion of PV output using a firm resource—assumed in this case to be a 4-hour lithium-ion battery with a power rating equal to 50% of the installed PV capacity. This definition highlights the context-specific nature of firming and distinguishes it from broader concepts such as dispatchability. Consequently, firming cost estimates must be interpreted with caution, as they depend heavily on local assumptions about resource adequacy, system needs, and market design.

The optimization of bidding strategies for PV plants has been extensively explored in the literature, with various studies addressing participation to different stages of the electricity market. However, very few studies have specifically analyzed the benefits of PV systems participation in IDM focusing on economic profitability and imbalance management. Additionally, while numerous optimization models incorporate the minimization of imbalance quantities and, consequently, imbalance costs, they rarely integrate this aspect with the optimization of BESS sizing. However, these two dimensions are strongly interrelated, as the extent of imbalances is inherently influenced by storage capacity, given that BESS can mitigate forecasting errors. This comes at the expense of investment and operational costs associated with storage infrastructure which must be carefully evaluated. Moreover, no study has assessed the difference between utility-scale and distributed PV configurations when it comes to market participation and imbalance management, with possible different impacts on the optimal BESS size. The technical and economic performance of distributed PV systems has been widely examined in scientific literature. However, direct comparisons with utility-scale PV power plants, which currently represent the predominant form of PV generation in the power system [36], remain limited.

1.3. Contribution and Novelty

This research investigates the optimal participation of PV + BESS systems within the Italian electricity market. Specifically, it focuses on the optimization of DAM and IDM bidding strategy through a multistage stochastic MILP, followed by a real-time operation phase, in which the flexibility of the battery is leveraged to minimize imbalances. The optimization framework is implemented in Python, using Pyomo library for mathematical modelling [37] and Gurobi as the solver [38]. To generate PV production scenarios, a Monte Carlo method is applied in combination with k-means clustering.

Two different system configurations are analyzed in this study:

• a distributed system of fixed rooftop PV systems located at the Leonardo campus of Politecnico di Milano [39];
• a simulated utility-scale plant situated at the same location.

The dataset for estimating daily PV power production profiles is sourced from the EU Photovoltaic Geographical Information System (PVGIS) [40]. To ensure comparability between the distributed and the utility-scale configurations, a fixed installed peak capacity of 1 MW is considered. By clustering 12 groups of days (corresponding to four seasons and three different irradiance conditions per season), the annual performance of the two configurations is assessed. The BESS size is varied from 0 to 5 MWh to assess its impact on both system’s economic performance and the resulting imbalance volumes.

This work contributes to the existing literature by providing a comparative analysis of two distinct PV system configurations—utility-scale and distributed—using a consistent methodological framework. While these configurations are not mutually exclusive, the study highlights their respective strengths and limitations in terms of profitability and dispatchability. The analysis evaluates system economic performance in electricity spot markets as a function of BESS capacity, offering insights into optimal BESS sizing and bidding strategies for each configuration. Furthermore, the study explores the trade-off between economic optimization and dispatchability requirements, assessing how an imposed imbalance threshold of 5% influence battery sizing and key economic metrics such as LCOE and NPV. This comprehensive approach provides a deeper understanding of the value of BESS integration in PV systems under varying operational constraints.

Although this study is conducted within the Italian electricity market and focuses on Milan’s PV generation, findings can be easily extended to other regulatory contexts and alternative geographical locations. The remainder of this paper is structured as follows: Chapter 2 provides an overview of the Italian electricity market, with a particular emphasis on the TIDE regulatory framework. Chapter 3 details the methodological approach, from scenario generation to the mathematical model used for optimization. Chapter 4 presents the case study and discusses the results obtained. Finally, Chapter 5 summarizes the key findings and outlines potential directions for future research.

2. Italian Electricity Market

This research considers the Italian context within the framework of the new TIDE regulatory initiative, ensuring alignment with the evolving European energy landscape. The section begins by outlining the core principles and objectives of TIDE, providing the necessary foundation for understanding its implications on market operations. It then offers an overview of the structure of the Italian electricity spot markets. The Ancillary Services Market (ASM) is excluded from the analysis, as it lies beyond the scope of this study.

2.1. TIDE

TIDE [7], introduced in January 2025, represents a significant step toward aligning the Italian electricity market with the EU regulatory landscape. Developed in compliance with EU Regulation 2017/2195 [41] and other European directives, TIDE establishes new principles and operational frameworks aimed at enhancing competition and liquidity in electricity markets. Particularly, TIDE aims to support the evolving needs of the electricity system as the share of non-dispatchable renewable energy sources continues to grow. Its primary objective is to implement an economically efficient dispatch model in which all network resources—including consumption units—can, in principle, serve a dual role: their primary function as energy producers or consumers, managed by Balance Responsible Parties (BRPs), and an additional function of providing flexibility services, managed by Balance Service Providers (BSPs). The model promotes competition among all types of resources, centralized or distributed, through organized market platforms that select the most efficient providers of flexibility services. This selection adheres to the principle of technological neutrality, meaning services are procured regardless of the providing technology as long as technical and cost criteria are met.

A key regulatory development introduced by TIDE is the formal recognition of Virtual Aggregated Units (Unità Virtuali Abilitate, UVA) as market participants, enabling aggregated DERs to compete alongside conventional power plants. This allows for greater flexibility in the participation of decentralized assets, ensuring that even small-scale resources can contribute to system balancing and market liquidity.

In summary, TIDE fosters a more inclusive and competitive electricity market by integrating RES, storage, and distributed systems and simultaneously enhancing system flexibility. This transformation includes a restructuring of market mechanisms, particularly within the spot markets, which are central to this study. The following sections outline the Italian electricity market structure, focusing on DAM, IDM and imbalance settlement within this new regulatory framework.

2.2. Day-Ahead Market and Intra-Day Market

The first phase of the Italian electricity spot market is the DAM, where market participants must submit their energy bids by 12:00 on the previous day (D-1), with market-clearing results published at 12:55 on D-1. Following the DAM, the IDM provides market participants with the opportunity to adjust their commercial positions by submitting additional supply offers or purchase requests. The IDM consists of three scheduled auctions (CRIDA) and a continuous trading session (IDM-XBID) [42].

The three intraday auction sessions are structured as follows:

• CRIDA1: Opens at 12:55 on D-1 and closes at 15:30 on D-1.
• CRIDA2: Opens at 12:55 on D-1 and closes at 22:00 on D-1.
• CRIDA3: Opens at 12:55 on D-1 and closes at 10:00 on D.

Prior to the introduction of TIDE, the resolution of spot market transactions was one hour, but starting in 2025, this transitioned to a 15-minute resolution to with the ASM. Consequently, a market time unit (MTU) of 15 minutes is adopted in this research.

In terms of spatial granularity, the Italian electricity market is divided into seven market zones: North (NORD), North-Central (CNOR), South-Central (CSUD), South (SUD), Calabria (CALA), Sicily (SICI), and Sardinia (SARD) [43].

In the DAM, operators submit bids that include: the MTU of delivery, the sale or purchase quantities and the respective price. Offers are accepted based on the auction results after the market closes. All supply offers accepted in the DAM are valued at the marginal clearing price of their respective zones. This price is determined for each MTU by the intersection of demand and supply curves and varies across zones when transit limits are saturated. With the introduction of TIDE, all purchasing offers are valued at the PUN Index, calculated by GME ex-post as the average of zonal prices, weighted by purchased quantities in each market zone.

In the CRIDA, demand bids and supply offers are selected on the basis of the same criterion as described for the DAM, employing a system marginal pricing method. In contrast, in the IDM-XBID, bids are matched immediately whenever a sell offer is priced lower than—or equal to—a corresponding buy bid.

This study considers a PV + BESS system that participates in the DAM and adjusts its commercial output in CRIDA2, based on improved PV forecast accuracy between 12:55 (D-1), the gate closure time (GTC) of the DAM, and 22:00 (D-1), the GTC of CRIDA2.

2.3. Imbalance Settlement

The nodal imbalance volume for a single unit is calculated for each Imbalance Settlement Period (ISP)—the time interval over which a BRP imbalance volume is financially settled—as the difference between the actual energy exchanged with the grid and the scheduled energy resulting from market clearing. This deviation, measured in MWh, indicates whether the unit has a positive (overproduction or underconsumption) or negative (underproduction or overconsumption) imbalance. In addition to nodal imbalances, macrozonal imbalances are also considered, which represent the net imbalance volume across all units within the same macrozone, which is defined as the aggregation of one or more pricing area.

Currently, in Italy is employed a static scheme with two macrozones: a North Macrozone corresponding to the ’NORD’ market zone and a South macrozone corresponding to the six remaining market zones. Imbalance volumes are calculated for each 15-minute ISP and are economically settled using the Single Pricing methodology. This pricing mechanism, established by ARERA (the Italian Regulatory Authority for Energy, Networks and Environment) through Deliberation 523/2021/R/eel [44], ensures alignment with European regulatory requirements [41]. Under the single pricing mechanism, the imbalance price for each ISP is calculated based on the overall macrozonal imbalance volume. Therefore, a uniform imbalance price is applied to all market participants in a given macrozone, regardless of whether they contribute positively or negatively to the overall imbalance. Particularly, the macrozonal imbalance price is determined as the volume-weighted average price of balancing energy activated by the TSO in response to system needs. Consequently, in a macrozone with positive imbalance, the imbalance price tends to be lower than the DAM price, whereas in a macrozone with negative imbalance, it tends to be higher than the DAM price. This design benefits participants unbalancing in the opposite direction of the macrozone imbalance volume: those who help mitigate it are rewarded with better prices than DAM prices, while those who worsen it face penalties. However, it must be noted that, since the imbalance price depends on real-time system conditions and balancing actions, it is not known during real-time operations. Instead, Terna, the Italian TSO, publishes the official imbalance prices on the day following market execution.

3. Methodology

This section outlines the modeling approach used to simulate the PV + BESS system operations and its interactions with electricity markets. The process begins with the generation of PV production scenarios using a Monte Carlo random sampling method, widely recognized in scientific literature for its flexibility in handling uncertainty and ease of implementation [45]. To manage computational complexity, a scenario reduction technique is also applied, specifically the k-means clustering algorithm. After generating PV scenarios for both DAM and IDM, the next step is to realistically model the knowledge of electricity prices, with particular attention to how imbalance prices are represented at different stages of the optimization. The final step involves the optimization of the PV + BESS system’s market participation through a multistage stochastic MILP. Given the inherent uncertainties in PV generation, the stochastic approach is considered the most suitable solution. The decision to adopt a multistage model is based on the progressive refinement of available information throughout the day, with improved PV production forecasts between DAM and IDM and the resolution of uncertainty during real time operation.

3.1. PV Scenario Generation

3.1.1. Day-Ahead Market PV Scenarios

Generating PV scenarios is essential to accurately capture the uncertainty associated with solar radiation during the considered day. As previously mentioned, the adopted methodology combines the Monte Carlo method with k-means clustering, following an approach similar to that described in reference [18]. To ensure a robust dataset, the study utilizes three years of historical PV production data with a 15-minute resolution.

The first step of the methodology involves classifying days by season. Within each seasonal group, k-means clustering is applied to group days based on solar radiation conditions. This clustering technique partitions a dataset of n observations into k clusters by assigning each observation to the cluster with the nearest mean, which serves as the cluster’s representative. The k-means clustering technique has been selected because it is particularly effective for extracting meaningful patterns from data, enabling the identification of trends in time-series behavior. By grouping similar data points, it reduces computational complexity while preserving the key statistical characteristics of the dataset [46].

To determine the optimal number of clusters for each season, the Elbow Method is applied. This technique, widely used in the literature, helps identifying the appropriate number of clusters by analyzing the within-cluster sum of squares (WCSS), which quantifies the total squared distance between each point in a cluster and its centroid. The Elbow Method involves plotting the WCSS on the y-axis against the number of clusters on the x-axis. The “elbow point”—where the decrease in WCSS starts to plateau—indicates the optimal number of clusters.

To maintain consistency and ensure computational efficiency, a common number of clusters is adopted across all four seasons. This choice is justified by the similarity in optimal cluster values obtained for each season and the benefits of uniformity in the subsequent analysis. As a result, the profiles are organized into 4×N matrices, where 4 represents the seasons and N denotes the number of clusters.

Once the clusters have been identified, the following steps are carried out for each cluster:

• Calculation of the average profile of the cluster being analyzed.
• Derivation of error profiles as the difference between the various historical profiles and the average profile, with a 15-minutes resolution.
• Application of the Monte Carlo method through a random addition of errors to the average profile. The selected errors must respect the time constraint (an error obtained in a given quarter hour can only be added to the PV production value of the average profile in the same quarter hour) but can belong to different error profiles. The process continues iteratively until the convergence criterion is satisfied.

Scenario generation is terminated when the uncertainty, defined as in eq. (1), reaches a value lower than the imposed threshold:

$${\xi }_{P^{PV}}\left(i\right)={\displaystyle \frac{1}{E_i \left[P^{PV}\right]}}\ast \sqrt{{\displaystyle \frac{\sum _{k=1}^i{\left(P^{PV}\left(k\right)-E_i \left[P^{PV}\right]\right)}^2}{i}}},$$

(1)

where i is the number of scenarios generated so far, $P^{PV}\left(k\right)$ represents the PV production profile for the k-th scenario, while $E_i \left[P^{PV}\right]$ denotes the average PV production profile over the first i scenarios.

Once the iterations are completed and all scenarios have been generated, a scenario reduction method is applied to limit computational efforts. The optimal number of reduced scenarios is determined using the elbow method, which is reapplied to assess the trade-off between the number of scenarios and the retained variance. After establishing a common number of reduced scenarios for each cluster, the k-means clustering method is employed to group similar scenarios based on daily PV power production. This process results, for each cluster in each season, in a final set of representative scenarios, each assigned a probability, ensuring that the reduced set accurately reflects the statistical properties of the original distribution.

3.1.2. Intra-Day Market Improved PV Scenarios

The main difference between DAM and IDM optimization lies in the improvement in PV forecast accuracy that occurs between the DAM GTC (12:55 on D-1) and the CRIDA2 GTC (22:00 on D-1). However, due to the lack of historical forecast data, this improvement is synthetically simulated. The steps applied, starting from the real profiles, remain the same, except for an additional step introduced between the Monte Carlo method and scenario reduction. Specifically, a subset of the newly generated profiles is selected based on a predefined filtering criterion. The goal is to retain only those scenarios that match more closely the actual PV production for the day being simulated—representative of the specific combination of season and irradiance conditions under analysis—while discarding those that deviate significantly. To achieve this, the normalized Root Mean Square Error (nRMSE) is used as the filtering metric, and an appropriate threshold is chosen to ensure that only the most realistic scenarios are retained. Equation (2) shows the calculation of the nRMSE:

$$nRMSE={\displaystyle \frac{\sqrt{\frac{1}{N}\ast \sum _{q=1}^N{\left(P_q^{gen}-P_q^{real}\right)}^2}}{\frac{1}{N}\ast \sum _{q=1}^N{P}_q^{real}}},$$

(2)

$N$ represents the number of time intervals (96 quarter-hour periods), $P_q^{gen}$ is the power at interval q of the generated scenario, and $P_q^{real}$ is the real power production.

This parameter was chosen as it effectively balances differences in both absolute value and shape, making it well-suited for assessing forecast accuracy. Its widespread use in the literature as a standard metric for evaluating forecasting model performance further supports its application in this study. For instance, [47] applies nRMSE to assess the accuracy of solar radiation forecasts, particularly in comparing different forecasting methodologies, while [48] analyzes various studies that frequently adopt nRMSE for performance evaluation, highlighting its relevance across different forecasting applications.

3.2. Mathematical Model

Once the scenarios have been generated, they serve as inputs for the mathematical optimization models that simulate the participation of the PV + BESS system in the electricity market. Figure 1 illustrates the sequence of the different optimization phases, along with the respective inputs and outputs. The process begins with the optimization of the DAM commercial program using the generated PV scenarios. This is followed by a refinement step in the IDM, leveraging more accurate PV forecasts. It is important to note that in both stages, the actual imbalance prices are not used. Instead, penalties are applied to both positive and negative imbalances to discourage intentional deviations, reflecting the operator’s uncertainty regarding future imbalance prices. Finally, based on the commercial program defined at IDM closure, the real-time operation of the plant is simulated. The objective in this phase is to minimize imbalances. At this stage, the simulation allows the calculation of daily profits using actual imbalance volumes and prices.

3.2.1. Day-Ahead Market

The initial PV scenarios, produced using the procedure outlined in section 3.1.1, along with DAM prices and imbalance penalties, are used as input for the MILP model that optimizes the bidding strategy in the DAM.

Equations (3)-(6) represent the objective function of the first stage of the mathematical model. As shown, the optimal solution corresponds to maximizing the DAM expected profit, while accounting for imbalance costs and battery costs.

$$\max \left(expected profit-imbalance cost-battery cost\right)$$

(3)

$$expected profit=\sum _{q=1}^{q=96}{€ }_q^{DAM}\ast {p }_q^{DAM}\ast ∆$$

(4)

$$imbalance cost=\sum _{i=1}^{i=n_{sc}}{\pi }_i\ast ( \sum _{q=1}^{q=96} {\lambda }_q^{im{b}_{negative}}\ast {p }_{i, q}^{im{b}_{positive}}\ast ∆ -\sum _{q=1}^{q=96} {\lambda }_q^{im{b}_{negative}}\ast {p }_{i,q}^{im{b}_{positive}} \ast ∆)$$

(5)

$$battery cost=\sum _{i=1}^{i=n_{sc}}{\pi }_i\ast {€}_{mean}^{DAM}\ast {(SOC}_{initial}-{soc}_{i, q=96})\ast {BESS}_{capacity}$$

(6)

The expected profit is calculated as the revenue derived from selling bids submitted in the DAM. Unlike the other two equations, this is not a weighted average across different possible outcomes, since the energy offered in the DAM remains the same across all PV scenarios. In contrast, imbalance quantities and BESS operations vary based on the assumed PV power production profile. The aim of the optimization is to define a strategy to be submitted to the DAM, representing the best possible compromise for all potential PV daily outcomes before uncertainty is resolved. The term associated with the DAM bid is multiplied by the time interval length (0.25 h) to maintain unit consistency. DAM prices are modeled under a price-taker assumption, given the high market liquidity, and are considered perfectly known in this study. This approach is commonly adopted in the literature for simplifying market modeling and optimizing bidding strategies [49]. Given that the system operates at a 15-minute resolution, whereas prices are still at an hourly resolution, a linear interpolation is applied to obtain intermediate values. This choice anticipates the changes expected under TIDE, which will standardize all markets to a single MTU of 15 minutes.

As described in section 2.2, Italy follows a single pricing system for imbalance settlement, meaning that imbalance prices depend on the real-time system position for each quarter-hour. However, if perfect foresight were assumed, this pricing mechanism could lead the optimizer to deliberately create imbalances in order to exploit potential advantages over the DAM price. Since such behavior is unrealistic, given that operators can hardly know imbalance prices in advance, penalties are introduced for both positive and negative imbalances during the DAM and IDM optimization phases. The calculation of the imbalance penalties follows these steps:

• Collection of annual DAM prices and imbalance prices.
• Calculation of the difference between imbalance prices and DAM prices over each quarter of an hour of the year.
• Computation of the annual average imbalance price deviations for each quarter-hour, separately for negative and positive system imbalances.
• Finally, each time a day is simulated, imbalance penalties are computed for each quarter-hour by adding the average imbalance price deviations to the actual daily DAM price profile.

Finally, the battery cost is a fictitious compensation term, preventing artificial profits from the energy initially stored in the battery. The calculation assumes that the cost of restoring the initial SOC is based on the DAM annual average price.

To ensure consistency between the physical power flows processed by the PV + BESS system and the commercial dispatch schedule, the following constraint must be imposed:

$${P }_{i,q}^{PV}-{p }_{i,q}^{charge}+{p }_{i,q}^{discharge}={p }_q^{DAM}-{p }_{i,q}^{im{b}_{negative}}\ast y_{i,q}^{im{b}_{negative}}+{p }_{i,q}^{im{b}_{positive}}\ast y_{i,q}^{im{b}_{positive}}.$$

(7)

In equation (7), the DAM bids remain constant across all scenarios, while the other variables are scenario-dependent. The two variables $y_{i,q}^{im{b}_{negative}}$ e $y_{i,q}^{im{b}_{positive}}$are binary variables that take a value of 1 depending on the type of imbalance. It is important to note that, in this study, the battery is assumed to be charged exclusively by PV power production, with no electricity being withdrawn from the grid.

BESS operations are governed by the following five constraints:

$${soc}_{i,q} \ge { SOC}_{min},$$

(8)

$${soc}_{i,q} \le {SOC}_{max},$$

(9)

$${soc}_{i,q}={soc}_{i,q-1}+∆\ast {\displaystyle \frac{ {\eta }_{charge}\ast {p }_{i,q}^{charge}\ast y_{i,q}^{charge}}{BES{S}_{capacity}}}-∆\ast {\displaystyle \frac{ {p }_{i,q}^{discharge}\ast y_{i,q}^{discharge}}{{\eta }_{discharge}\ast BES{S}_{capacity}}},$$

(10)

$${p }_{i,q}^{charge} \le BES{S}_{Pmax},$$

(11)

$${p }_{i,q}^{discharge} \le BES{S}_{Pmax}.$$

(12)

The SOC of the battery can neither exceed 100% nor fall below the 0% threshold in any scenario for any quarter-hour. Excess or insufficient power is converted into positive and negative imbalances, respectively, thereby introducing an economic penalty. Equation (10) defines the SOC update for each quarter-hour based on the BESS processed power. Two binary variables $y_{i,q}^{charge}$ and $y_{i,q}^{discharge}$ are defined to represent if the BESS is charging or discharging. The final two constraints ensure that power flows remain within the maximum allowable limits imposed by the technical specifications of the BESS.

3.2.2. Intra-Day Market

The second stage of the stochastic MILP model consists of solving the IDM. The primary distinction between the two phases is the improvement in PV forecasts. The newly refined forecast scenarios are incorporated as inputs into the model, alongside IDM prices, imbalance penalties, and the bids already submitted in the DAM.

The methodology remains largely unchanged. The objective function follows the same structure of eq. (3), with the only variation concerning the expected profit component, as shown in eq. (13):

$$expected profit=\sum _{q=1}^{q=96}{€ }_q^{DAM}\ast {P }_q^{DAM}\ast ∆+\sum _{q=1}^{q=96}{€ }_q^{IDM}\ast {p }_q^{IDM}\ast ∆.$$

(13)

In this stage, the additional revenue from IDM transactions is incorporated, complementing the earnings already secured in the DAM. The IDM prices are assumed to be known in advance for the same reasons outlined in the first stage. Given that this study focuses exclusively on participation in the CRIDA2 auction, the corresponding prices are used. To achieve adjustment to the DAM submitted commercial program, IDM bids can take both positive and negative values. A positive IDM bid reflects an additional purchase or an increase in energy injection to compensate for an underestimated DAM position, whereas a negative IDM bid corresponds to a reduction in energy injection or an increase in energy withdrawal, correcting an overestimated DAM position. These adjustments must always remain within the nominal power limits of the system, ensuring feasibility. Additionally, the sum of DAM and IDM bids must never be negative, preserving operational consistency.

Finally, the power balance equation is modified as follows:

$${P }_{i,q}^{PV}-{p }_{i,q}^{charge}+{p }_{i,q}^{discharge}={P }_q^{DAM}+{p }_q^{IDM}-{p }_{i,q}^{im{b}_{negative}}+{p }_{i,q}^{im{b}_{positive}}.$$

(14)

3.2.3. Real-Time Operation

The last stage of the mathematical model does not involve an optimization process but rather an operational strategy for real-time management. At this stage, the actual PV power production for each quarter-hour is revealed, fully eliminating uncertainty. Unlike the previous phases, where forecasts and stochastic scenarios influenced decision-making, the real-time strategy relies solely on known values, with no probabilistic components to consider. At each quarter-hour, the model minimizes imbalance quantities as defined in the objective function expressed in eq. (15):

$$\underset{q}{\min }\left(p_q^{im{b}_{negative}}+p_q^{im{b}_{positive}}\right).$$

(15)

Rather than optimizing for the entire day at once, the model operates sequentially, updating at each time step q. The battery SOC is carried forward from one interval to the next, reflecting a realistic operational framework in which PV production is not known in advance but uncertainty is resolved progressively. To reduce imbalances, the model utilizes the battery’s flexibility. However, due to its physical constraints, long or short imbalances may still occur when the battery reaches either its maximum or minimum SOC.

4. Case Studies and Results

4.1. Analyzed Case Studies

PV production data is sourced from PVGIS [40], a widely used tool that offers significant flexibility in various aspects:

• It allows users to specify the exact geographical location of the PV system; in this study, Piazza Leonardo da Vinci in Milan is selected.
• The time period can vary between 2005 and 2024. In this study, the range 2021–2023 is chosen to ensure a comprehensive representation of production profiles.
• The system type can be customized by selecting fixed panels or tracking systems, allowing further customization of tilt and orientation.
• The installed peak power and system losses can be defined.

The distributed PV power plant analyzed in this study consists of eight fixed rooftop installations located across the Politecnico di Milano Leonardo Campus, with an overall peak power of 1MW and the characteristics summarized in Table 2. The systems employ three main mounting configurations, reflecting the architectural diversity and structural constraints of the selected rooftops. The majority of the systems are rooftop-coplanar, where PV modules are installed flush with the slope and orientation of the existing roof surfaces. Two installations use a ballasted mounting system on flat roofs; in this configuration, modules are secured by weight rather than mechanical anchoring, enabling flexible placement without compromising the roof structure [50]. Finally, one installation is mounted on a barrel-shaped roof, approximated using a simplified model with three tilt sections to represent the curved geometry:

• 50% of the panels positioned at 0° inclination,
• 25% at 15° inclination,
• 25% at -15° inclination (equivalent tilt but opposite azimuth).

Additionally, Figure 2 illustrates the geographical distribution of these rooftop PV systems providing an overview of their spatial layout within the university campus. However, it is important to note that, due to the close proximity of the buildings, PVGIS does not account for localized differences in solar irradiance across the campus. As a result, all rooftop systems are treated as if they are located at the same geographical point, and temporal variations in cloud cover across individual buildings are not captured. Despite this limitation, each system retains its unique tilt and azimuth values, which result in distinct PV production profiles even under the same irradiance conditions. The losses were estimated by observing the actual performances of the systems installed across the university.

Once the system parameters are defined, PVGIS is used to simulate daily PV production profiles for the period 2021–2023. The output from the eight rooftop systems is then aggregated to generate the overall 1 MW production profile representing the distributed PV configuration.

The utility-scale plant considered in this study is a hypothetical ground-mounted system, located at the same geographical site as the rooftop installations and designed with the same peak power of 1 MW to enable a consistent comparison. Unlike the fixed rooftop systems, this configuration employs a horizontal single-axis tracking system, which allows the PV modules to rotate along a north–south axis. This tracking mechanism follows the sun’s movement from east to west throughout the day, increasing solar exposure and thereby enhancing power production compared to fixed systems [51]. The system is assumed to be installed on flat ground, typical of utility-scale deployments, and system losses are set to 14%, in line with the default value in PVGIS and the assumptions made in ref. [21].

The original PVGIS dataset has an hourly resolution. However, given the 15-minutes resolution required for this study, a linear interpolation method was applied. As a result, two matrices were constructed, each with dimensions 1095 × 96 (corresponding to three years of daily production profiles, with a 15-minute resolution). These matrices represent the complete production datasets for both the rooftop and utility-scale systems.

Once the two matrices have been obtained, they are divided into four seasonal subsets. For each season, the Elbow method is applied to determine the optimal number of clusters for further partitioning the sub-matrices. Figure 3 provides a graphical example of this process for Autumn.

The analysis identifies three optimal clusters per season, each corresponding to different total radiation levels, influenced by seasonal variations and climatic conditions. Particularly, clusters were defined as Sunny (highest production, clear-sky conditions), Variable (intermediate production, mixed weather), and Cloudy (lowest production, cloudy). This classification led to 4 x 3 = 12 clusters per plant configuration. From each of these 12 clusters, one representative day from 2023 is selected to simulate the market participation of the PV + BESS system. Rather than simulating all days in the year, the analysis is performed only on these 12 representative days, each assumed to typify the operational and production conditions of its corresponding cluster. This methodological choice was driven by both computational and time constraints, allowing for a feasible simulation process while still preserving the seasonal and weather-related variability in solar production. The selected representative days are shown in Table 3.

To reconstruct a full-year performance profile from these limited simulations, an occurrence probability is assigned to each cluster based on how frequently it occurs within the three analyzed years. The results obtained from each simulated representative day are then weighted by this probability and scaled to 365 days, effectively allowing the model to approximate annual performance and market participation outcomes from a limited set of representative days. To verify the effectiveness of the method adopted, the Capacity Factor (CF) is calculated as 18.6% for the utility-scale plant and 14.2% for the distributed system, values that are consistent with real-world performance expectations for similar PV installations.

4.2. PV Scenario Generation and Market Data

The first phase of simulating market participation involves generating PV scenarios for DAM using the Monte Carlo method, followed by scenario reduction through k-means clustering, as described in Section 3.1. This process is carried out separately for each selected representative day. It is important to clarify that the scenario generation method from Section 3.1 is applied using all PV production profiles within the same cluster—meaning all days from the 2021–2023 period that fall under the same category (Sunny, Variable, or Cloudy) and season as the representative day.

To ensure convergence in scenario generation, an uncertainty threshold of 1% is imposed. If this threshold is not reached, the process is automatically terminated after generating 2000 scenarios. These parameters are set to strike a balance between accuracy and computational effort, in line with the values used in similar works [18].

After completing the scenario generation phase, the number of profiles is reduced using the Elbow method to identify representative cases. Unlike the seasonal clustering, the WCSS curves do not display a clear elbow point, making it challenging to determine the optimal number of scenarios. Nevertheless, six clusters are selected, as this number provides the best trade-off between capturing the key variability in the data and maintaining computational efficiency, with reference to the approach adopted in [18]. As a result, six representative PV production profiles are selected for each of the 12 representative days and for both plant configurations. These profiles represent the forecasted possible PV power generation scenarios used in the DAM simulation.

To generate the scenarios for the IDM simulation, as detailed in section 3.1.2., a filtering criterion is introduced. A nRMSE threshold is imposed, varying for each cluster and plant configuration. This differentiation is necessary because a fixed threshold would not be equally effective across all conditions. A threshold that is highly selective for a cloudy day might fail to filter any Monte Carlo-generated scenarios for a sunny day. Conversely, a threshold suitable for a sunny day would eliminate almost all scenarios generated for variable or cloudy clusters. Similarly, the two plant configurations require distinct nRMSE thresholds for the same cluster. To ensure uniformity in the filtering process across all days, the nRMSE threshold is set to retain between 1% and 5% of the generated scenarios. The final step in scenario processing involves further reducing the filtered scenarios using k-means clustering, obtaining six representative scenarios, matching the number used for the DAM, to serve as inputs for the IDM optimization stage.

Figure 4 provides an example of the various processing steps applied to four representative days. The thicker blue line shows the actual PV production on the representative day for the given cluster and is therefore identical for both the DAM and IDM stages. What differs are the six reduced PV scenarios used as inputs for the two stages of the stochastic optimization: the IDM scenarios, produced after the filtering process, are much closer to the actual production profile, while the DAM scenarios show greater deviations since they result from non-filtered Monte Carlo profiles.

Finally, it is necessary to collect real market price data to provide all the required inputs for the multistage stochastic optimization. The real DAM and IDM prices for the analyzed days were collected from the Italian Energy Market Operator (GME) [52]. The imbalance prices were instead retrieved from Terna’s official database [53].

4.3. Techno-Economic Analysis of PV + BESS Market Participation

The objective of this study is to determine the optimal battery size from both a profitability and a dispatchability perspective. To achieve this, it is crucial to obtain results across a wide range of storage capacities, allowing for a comprehensive evaluation. To avoid excessive computational effort, the Energy-to-Power Ratio (EPR) is fixed at 2. Consequently, storage capacity is varied between 0 and 5 MWh, in steps of 100kWh. This results in a total of 4 (seasons) × 3 (clusters) × 51 (battery sizes) = 612 market participation simulation, repeated for both plant configurations. At each iteration, the three mathematical models described in Section 3.2. are solved and optimized. Finally, the total profit is computed as follows:

$${profit}^{total}={profit}^{DAM}+{profit}^{IDM}-{cost}^{imb}-{cost}^{fina{l}_{SOC}},$$

(16)

$${profit}^{DAM}=\sum _{q=1}^{96}{€}_q^{DAM}\ast p_q^{DAM}\ast ∆,$$

(17)

$${profit}^{IDM}=\sum _{q=1}^{96}{€}_q^{IDM}\ast p_q^{IDM}\ast ∆,$$

(18)

$${cost}^{imb}=\sum _{q=1}^{96}{€}_q^{imb}\ast P_q^{im{b}_{negative}}-{€}_q^{imb}\ast P_q^{im{b}_{positive}},$$

(19)

$${cost}^{fina{l}_{SOC}}={€}_{mean}^{DAM}\ast {(SOC}_{initial}-{SOC}_{q=96})\ast {BESS}_{capacity},$$

(20)

where $P_q^{im{b}_{negative}}$ and $P_q^{im{b}_{positive}}$are the residual imbalances following the system real-time operations.

The total profit is initially calculated for each of the 12 representative days across the different battery sizes. To estimate the annual total profit for each BESS capacity, these daily profits are combined using a weighted average. The weights correspond to the occurrence probabilities of each cluster, which represent how frequently each type of day (Sunny, Variable, or Cloudy) in each season (Winter, Spring, Summer, Autumn) occurs during the 2021–2023 period used to generate the PV scenarios. The occurrence percentages of each cluster, for both the utility-scale plant and the distributed PV system, are presented in Table 4. It is important to note that these percentages differ slightly between the two configurations because the clustering process was performed separately based on their respective PV production profiles. However, since both systems are modeled at the same geographical location, the percentages are quite similar, reflecting consistent radiation level patterns.

In the following sections, we present the results of our analysis assessing the economic impact and improved dispatchability achieved through BESS installation for both configurations. The analysis consists of three steps. First, we determine the minimum BESS size required to reduce annual imbalance volumes—relative to the PV plant’s total production—to below 5%. This threshold is based on [32], which reports that imbalances in Italy accounted for 5.6% of national demand in 2016. To isolate the effect on dispatchability, we run simulations using constant DAM and IDM prices across both stages, so that the optimization focuses solely on minimizing imbalances, excluding any gains from energy arbitrage. This allows us to estimate the net cost of making the PV system dispatchable. In the second step, recognizing that a BESS cannot be financially justified if used solely to reduce imbalances, we assess the system economic performance when the storage is optimally used also for energy arbitrage in DAM and IDM. This allows us to identify the optimal BESS size and associated revenues for both the utility-scale and distributed configurations. Finally, we combine the two perspectives to evaluate the optimal BESS size when a 5% imbalance constraint is enforced. The economic analysis is based on two indicators: the LCOE, which is used to capture the system cost increase due to BESS installation; and the NPV, which estimates the economic returns over the system’s lifetime.

The LCOE is computed as in eq. (21) over a 30-year horizon aligned with the PV plant’s expected lifetime:

$$LCOE={\displaystyle \frac{{CAPEX}_{PV}+\sum _{t=1}^{30}\frac{{OPEX}_{PV}}{{\left(1+i\right)}^t}+\sum _{t=0}^2\frac{{CAPEX}_{BESS}}{{\left(1+i\right)}^{10t}}+\sum _{t=1}^{30}\frac{{OPEX}_{BESS}}{{\left(1+i\right)}^t}}{\sum _{t=1}^{30}\frac{{PVPower}_{sold}}{{\left(1+i\right)}^t}}}.$$

(21)

PV CAPEX and OPEX are taken from [54]: the Commercial PV category is used for the distributed system (CAPEX = €1300k/MW, OPEX = €15k/MW), while the Utility Scale category is used for the utility-scale system with tracking (CAPEX = €950k/MW, OPEX = €17.5k/MW). BESS is assumed to have a 10-year lifetime, requiring two replacements over 30 years and ending without any residual value. BESS CAPEX is modeled as in [55], with an energy-related cost of €250k/MWh (battery banks) and a power-related cost of €80k/MW (inverter and grid integration). OPEX is estimated at €5k/MWh/year. The LCOE denominator includes only the energy exchanged with the grid—i.e., net of BESS round-trip losses in PV + BESS configurations.

Based on these parameters and assuming a 5% discount rate—consistent with values used in similar studies [56,57,58]—the NPV can be calculated over a 10-year horizon as follows:

$$NPV={-CAPEX}_{BESS}+\sum _{t=1 }^{10}{\displaystyle \frac{{\Delta }_{Revenues} \left[\frac{€}{year}\right]-OPE{X}_{BESS}}{{(1+i)}^t}} .$$

(22)

This NPV formulation evaluates whether the BESS investment adds economic value to the system. It accounts only for the BESS-related CAPEX and OPEX, and the incremental annual revenues generated with BESS compared to the baseline case without storage.

4.3.1. Economic Impact of Dispatchability Constraints

Assuming that only the PV plant participates in the DAM and IDM and using the PV production scenarios described in Section 4.2., the annual imbalance volume for the utility-scale configuration with horizontal single-axis tracking amounts to approximately 13% of total energy production, while for the distributed configuration it is around 11%. If a BESS is added solely to reduce imbalance volumes below the 5% threshold, the required storage capacity is 200 kWh for the utility-scale system and 100 kWh for the distributed one, reducing imbalance levels to 4.7% and 3.8%, respectively.

Figure 5 shows the LCOE for both configurations in the base case (without BESS) and in the firmed cases described above. The cost of firming is approximately €6/MWh for the utility-scale system, resulting in a 12.6% increase in LCOE, and around €3/MWh for the distributed system, with a 3.8% increase.

However, it is important to highlight that installing such battery sizes solely to minimize imbalances does not generate sufficient revenues to recover the BESS investment costs. Specifically, calculating the NPV using the formulation in Eq. (22), under current BESS costs, yields –55k€ for the utility-scale system and –42k€ for the distributed one, highlighting the need to explore additional uses for the BESS, such as its potential for arbitrage in the DAM and IDM.

4.3.2. Unconstrained Optimal BESS Sizing

This section analyzes the optimal operation of the PV+BESS system in electricity markets under varying BESS size and without imposing any constraints on imbalance volumes. Figure 6 presents the annual profit trend and its derivative as battery capacity increases.

As expected, the profit generated by the PV + BESS system increases with battery size. However, the rate of increase is not uniform: for smaller storage capacities, profit grows significantly, whereas for larger capacities, the marginal benefit decreases. This decreasing revenue suggests that increasing BESS capacity beyond a certain threshold no longer provides proportional economic benefits.

To assess the economic feasibility of BESS installation we compute the NPV as in eq. (22). The results indicate that, for any battery size, the NPV remains negative. This suggest that, under 2023 electricity market conditions and assumed battery costs, for a participation limited to DAM and IDM, a battery does not generate sufficient revenue to fully recover its costs. From a strictly economic perspective, the most financially advantageous solution is to operate without a BESS, for both system configurations.

Therefore, we identify the maximum battery bank cost at which at least one PV+BESS configuration remains economically viable, defined as achieving a positive NPV. This threshold is approximately half the current cost, around 125 k€/MWh. Figure 7 shows the total 10-year revenue for each battery size after accounting for this reduced BESS cost. The optimal battery sizes that maximize net profit fall within the 1–2 MWh range for both configurations, though differences in profitability within this range are minimal and difficult to distinguish visually. Specifically, the highest net profit is achieved with a 1.4 MWh BESS for the utility-scale system and a 1.3 MWh BESS for the distributed system. For comparison, under current BESS costs, these same configurations would result in NPVs of –239 k€ and –221 k€, respectively.

The shaded areas in Figure 7 highlight the configurations where the PV+BESS system outperforms the PV-only baseline under the assumed BESS cost scenario. These areas largely overlap, with the utility-scale system exhibiting a slightly broader profitable range. In both cases, installing more than 2 MWh of storage does not provide economic benefits.

4.3.3. Optimal BESS Sizing Under Dispatchability Constraint

This study also aims to investigate how the dispatchability of the system is influenced by increasing battery size. To achieve this, annual normalized imbalance volumes under optimal market participation are presented for each considered BESS capacity for both the utility-scale plant and the distributed system.

Figure 8 illustrates the normalized imbalance volumes as BESS capacity increases. The graph shows that the utility-scale plant exhibits higher normalized imbalanced volumes than the distributed system for most of the BESS size range.

To investigate the reasons behind the imbalance volume differences between the two configurations, we analyze the days with the largest imbalance discrepancies. For this analysis, we select a BESS capacity of 1.4 MWh, identified in the previous section as an optimal size. With this fixed capacity, we collect the seven components of the system’s energy balance at 15-minute resolution across 12 representative days: solar production, charging power, discharging power, DAM bids, IDM bids, and both short and long imbalance volumes. Figure 9 and Figure 10 display the corresponding diagrams for the utility-scale and the distributed configurations operating on the representative autumn day with variable conditions.

It can be seen that the horizontal tracking system of the utility-scale plant produces sharper profiles with two peaks in power production throughout the day. In contrast, the distributed system generally exhibits smoother profiles, with a smaller volume of generated energy. This is further illustrated in Figure 11, which shows the percentage difference in annual PV production between the two systems as function of time. It is clear that the largest production difference occurs at the beginning and end of the day, with the distributed system generating significantly less energy during these periods. This suggests that the utility-scale plant benefits from a longer daily production window, enabling a wider range of production patterns and introducing greater forecast uncertainty. As a result, there is a higher likelihood of committing to a commercial schedule that leads to unavoidable imbalances.

Referring back at Figure 8, we observe that the minimum BESS sizes required to meet the imposed 5% imbalance threshold, under optimal BESS market participation, are 1.1 MWh for the distributed system and 1.7 MWh for the utility-scale system. This indicates that the dispatchability constraint does not limit the economic optimum for the distributed configuration, as its optimal BESS size (1.3 MWh) exceeds the minimum required (1.1 MWh). In contrast, for the utility-scale system, the optimal size (1.4 MWh) falls below the 1.7 MWh needed for dispatchability, meaning that enforcing this constraint reduces profitability.

As a final step, we aim to determine the cost of making a PV system dispatchable defined as achieving an annual imbalance volume below 5% of total annual production under the condition that the installed BESS must also engage in energy arbitrage to recover its investment costs. Therefore, Figure 12 compares the LCOE of the PV-only system with that of the PV+BESS systems for both the utility-scale and distributed configurations, using the minimum BESS sizes required to keep the relative imbalance below the 5% threshold (respectively 1.7MWh and 1.1MWh) and employing current BESS costs.

It is clear that the LCOE is significantly higher for BESS-integrated solutions in both configurations. For the utility-scale system, the LCOE increases from 46€/MWh to 95€/MWh, while for the distributed system, it rises from 76€/MWh to 119€/MWh. Therefore, if the BESS is also used to maximize system profit through spot market arbitrage, the relative increase in LCOE amounts to 104% for the utility-scale system and 57% for the distributed one.

The LCOE increase reflects the cost of firming which the two systems must bear to minimize imbalances at a negligible level when performing energy arbitrage. If the additional profits generated by a larger BESS compensate for this cost increase, the system achieves an optimum in both profitability and dispatchability. For the PV + BESS configurations reported in Figure 12, the NPV under current BESS costs is negative in both cases and worse than in the pure firming scenario presented in Section 3.1: –295 k€ for the utility-scale system compared to –55 k€, and –191 k€ for the distributed system compared to –42 k€. However, assuming reduced BESS costs, as discussed in Section 3.2, could make larger BESS capacities more economically viable.

It is also important to note that this analysis considers only revenues from participation in the day-ahead and intraday markets. Additional revenue streams—such as ancillary services—could further improve the economic performance, particularly if the battery is underutilized in the current configuration. To evaluate this, we compute the battery’s Capacity Factor (CF) as follows:

$$CF \left[\%\right]={\displaystyle \frac{\sum _{d=1}^{12}{\pi }_d\sum _{q=1}^{96}{P}_{d,q}^{charge}\ast ∆+P_{d,q}^{discharge}\ast ∆ }{{BESS}_{P max}\ast 24h}}\ast 100 .$$

(23)

Here, d represents the 12 selected typical days. A weighted average of the daily contributions is used to determine the total energy exchanged by the battery over an average day. The denominator accounts for the maximum theoretically exchangeable energy in a single day. For battery sizes ensuring imbalance percentages around 5%, the calculated CF values are 18.6% for the utility-scale plant and 14.2% for the distributed system. These values confirm that the battery is not fully utilized across all operational hours, indicating strong potential for additional service provision.

5. Discussion

This study develops a multistage stochastic optimization framework to assess the participation of a 1 MW PV system coupled with a BESS in the Italian electricity spot markets, accounting for both DAM and IDM operations. Two configurations—utility-scale and distributed—are analyzed to determine the optimal battery size for each.

To explore the technical and economic interplay between storage sizing and system performance, the study investigates the trade-off between dispatchability and profitability. When considering PV-only systems, the LCOE is calculated at 46.39 €/MWh for the utility-scale system and 75.83 €/MWh for the distributed one. Introducing a BESS used solely to meet the 5% annual imbalance threshold—without optimizing for market revenues—requires 200 kWh of storage for the utility-scale system and 100 kWh for the distributed system. This increases the LCOE to 52.24 €/MWh and 78.77 €/MWh, respectively, which remain significantly below the firmed PV LCOE estimated by Lazard in the CAISO context (141 $/MWh) [34]. For a better comparison, it should be noted that, for the same PV peak power, Lazard’s study assumes a 2MWh BESS. However, this strategy proves economically inefficient with respect to PV stand-alone, as minimizing imbalances alone does not yield sufficient revenue to offset current storage costs, resulting in negative NPVs of –55 k€ and –42 k€, respectively.

An alternative approach is to leverage the BESS for energy arbitrage in both DAM and IDM, while also minimizing imbalances in real time. Yet, under current market conditions and battery costs, BESS operation results in negative NPVs across all tested sizes—making PV-only systems the most economically viable option. Nevertheless, future reductions in battery CAPEX, as well as access to additional revenue streams (e.g., frequency regulation services), could significantly improve the financial viability of PV + BESS systems.

To simulate such a future scenario, we assume a 50% reduction in battery bank investment costs—from 250 k€/MWh to 125 k€/MWh—representing the break-even threshold at which at least one BESS size yields a positive NPV for both configurations. In this scenario, the optimal BESS size (without dispatchability constraints) becomes 1.4 MWh for utility-scale and 1.3 MWh for distributed systems, yielding slightly positive NPVs. For reference, under current BESS costs, the same configurations would result in NPVs of –239 k€ and –221 k€.

Finally, we reintroduce the 5% imbalance constraint while maintaining energy arbitrage strategies. For the utility-scale system, this raises the minimum required BESS size to 1.7 MWh, resulting in a return to negative NPV, even under reduced battery costs. In contrast, the distributed system with a 1.3 MWh BESS remains within the dispatchability threshold, as the minimum required size to meet firm operation is 1.1 MWh. Therefore, no resizing is needed.

Calculating the LCOE under current BESS costs for these minimum dispatchable configurations yields 94.53 €/MWh for the utility-scale system with a 1.7 MWh BESS, and 118.94 €/MWh for the distributed system with a 1.1 MWh BESS. These higher LCOE values reflect the trade-off between system firmness and storage utilization. Using the BESS solely for imbalance minimization would require smaller storage capacities, resulting in lower LCOE but also negative NPVs. Conversely, leveraging the BESS for energy arbitrage—assuming it generates enough revenue to offset its cost—leads to higher LCOE values, but potentially enables positive NPVs.

These findings highlight that distributed systems—due to smoother production profiles—benefit from improved forecast accuracy when using a Monte Carlo scenario generation approach. This allows them to achieve dispatchability with smaller storage capacities. In contrast, utility-scale systems face more variability and require larger batteries to meet the same dispatchability targets, making investment recovery more difficult.

Several directions for future research emerge from this work. First, incorporating the ancillary services market, and particularly the provision of frequency regulation services, would provide a more complete assessment of BESS revenue streams. Second, battery degradation modeling should be introduced to evaluate long-term performance and cost recovery. Third, considering a variable BESS EPR could provide deeper insights into the balance between storage capacity and power availability. Finally, scenario generation could be refined by using historical PV forecasts and actual plant data, allowing for more realistic DAM-IDM transitions and better modeling of forecast uncertainty across time horizons.