Hull Clustering with Blended Representative Periods for Energy System Optimization Models

TL;DR

This paper introduces hull clustering with blended representative periods, a novel method that improves the accuracy and efficiency of energy system planning models by better capturing demand and renewable variability.

Contribution

It proposes a new clustering approach using extreme points and blended RPs to enhance traditional representative period methods in energy system optimization.

Findings

Outperforms traditional RP techniques in accuracy and computational efficiency.

Reduces the number of RPs needed for accurate modeling.

Demonstrated effectiveness on European energy data.

Abstract

The growing integration of renewable energy sources into power systems requires planning models to account for not only demand variability but also fluctuations in renewable availability during operational periods. Capturing this temporal detail over long planning horizons can be computationally demanding or even intractable. A common approach to address this challenge is to approximate the problem using a reduced set of selected time periods, known as representative periods (RPs). However, using too few RPs can significantly degrade solution quality. In this paper, we propose the method of hull clustering with blended RPs to enhance traditional clustering-based RP approaches in two key ways. First, instead of selecting typical cluster centers (e.g., centroids or medoids) as RPs, our method is based on extreme points, which are more likely to be constraint-binding. Second, it represents base periods as weighted combinations of RPs (e.g., convex or conic blends), approximating the full time horizon more accurately and with fewer RPs. Through two case studies based on data from the European network operators, we demonstrate that hull clustering with blended RPs outperforms traditional RP techniques in both regret and computational efficiency.

Equations

$$\mathcal{A}=\mathcal{G}\bigsqcup \mathcal{S}\bigsqcup \mathcal{C},\text{where}\mathcal{S}={\mathcal{S}}^{\text{ST}}\bigsqcup {\mathcal{S}}^{\text{seas}}.$$

$$\text{if }\rho (d)=r\text{ then }{W}_{d,r}=1,\text{ otherwise }{W}_{d,r}=0.$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle}{\mathrm{subject~to}}$}][l]{\underset{\displaystyle}{\mathrm{minimize}}}\quad c^{\text{inv}}+c^{\text{op}}$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle\phantom{}}{\mathrm{subject~to}}$}][c]{{\mathrm{subject~to}}}$$

$$c^{\text{inv}}=\text{smashoperator}[]\sum _{a^{\text{inv}}\in {\mathcal{A}}^{\text{inv}}}{I}_{a^{\text{inv}}}{P}_{a^{\text{inv}}}{i}_{a^{\text{inv}}}$$

$$\displaystyle c^{\text{op}}=W^{\text{op}}\sum_{r\in\mathcal{R}}W_{\!r}\sum_{h\in\mathcal{H}}\Bigl(\sum_{g\in\mathcal{G}}V_{g}p^{\text{out}}_{g,r,h}+\makebox[0.0pt][l]{$\tfrac{1}{\tau}\smashoperator[]{\sum_{s^{\text{se}}\in\mathcal{S}^{\text{seas}}}^{}}\bigl(V^{\text{sp}}_{s^{\text{se}}}e^{\text{sp}}_{s^{\text{se}},r,h}+V^{\text{br}}_{s^{\text{se}}}e^{\text{br}}_{s^{\text{se}},r,h}\bigr)\Bigr)$}$$

$$\text{smashoperator}[]\sum _{a\in {\mathcal{A}}_{n,x}^{\text{in}}}{p}_{a,r,h}^{\text{out}}-\text{smashoperator}[]\sum _{a\in {\mathcal{A}}_{n,x}^{\text{out}}}{p}_{a,r,h}^{\text{in}}-\text{smashoperator}[]\sum _{l\in {\mathcal{L}}_{n,x}^{\text{out}}}{f}_{l,r,h}+\text{smashoperator}[]\sum _{l\in {\mathcal{L}}_{n,x}^{\text{in}}}{f}_{l,r,h}=D_{n,x,r,h}{D}_{n,x}^{\text{max}}$$

$${\sigma }_{s^{\text{ST}},r,1}^{\text{intra}}-{\sigma }_{s^{\text{ST}},r}^{\text{intra,0}}=({\eta }_{s^{\text{ST}}}^{\text{in}}{p}_{s^{\text{ST}},r,1}^{\text{in}}-\frac{p_{s^{\text{ST}},r,1}^{\text{out}}}{{\eta }_{s^{\text{ST}}}^{\text{out}}})\tau$$

$${\sigma }_{s^{\text{ST}},r,h}^{\text{intra}}-{\sigma }_{s^{\text{ST}},r,h-1}^{\text{intra}}=({\eta }_{s^{\text{ST}}}^{\text{in}}{p}_{s^{\text{ST}},r,h}^{\text{in}}-\frac{p_{s^{\text{ST}},r,h}^{\text{out}}}{{\eta }_{s^{\text{ST}}}^{\text{out}}})\tau$$

$$\displaystyle\sigma^{\text{intra}}_{s^{\text{se}},r,1}-\sigma^{\text{intra,\,0}}_{s^{\text{se}},r}=\makebox[0.0pt][l]{$\bigl(\eta^{\text{in}}_{s}p^{\text{in}}_{s^{\text{se}},r,1}-\frac{p^{\text{out}}_{s^{\text{se}},r,1}}{\eta^{\text{out}}_{s^{\text{se}}}}\bigr)\tau-e^{\text{sp}}_{s^{\text{se}},r,1}+e^{\text{br}}_{s^{\text{se}},r,1}+E_{s^{\text{se}},r,1}E^{\text{max}}_{s^{\text{se}}}$}$$

$$\forall s^{\text{se}},r,$$

$$\displaystyle\sigma^{\text{intra}}_{s^{\text{se}},r,h}-\sigma^{\text{intra}}_{s^{\text{se}},r,h-1}=\makebox[0.0pt][l]{$\bigl(\eta^{\text{in}}_{s^{\text{se}}}p^{\text{in}}_{s^{\text{se}},r,h}-\frac{p^{\text{out}}_{s^{\text{se}},r,h}}{\eta^{\text{out}}_{s^{\text{se}}}}\bigr)\tau-e^{\text{sp}}_{s^{\text{se}},r,h}+e^{\text{br}}_{s^{\text{se}},r,h}+E_{s^{\text{se}},r,h}E^{\text{max}}_{s^{\text{se}}}$}$$

$$\forall s^{\text{se}},r,h\ne 1,$$

$${\sigma }_{s^{\text{se}},1}^{\text{inter}}-{\sigma }_{s^{\text{se}}}^{\text{inter,0}}=\text{smashoperator}[]\sum _{r\in \mathcal{R}}{W}_{1,r}({\sigma }_{s^{\text{se}},r,H}^{\text{intra}}-{\sigma }_{s^{\text{se}},r}^{\text{intra,0}})$$

$${\sigma }_{s^{\text{se}},d}^{\text{inter}}-{\sigma }_{s^{\text{se}},d-1}^{\text{inter}}=\text{smashoperator}[]\sum _{r\in \mathcal{R}}{W}_{d,r}({\sigma }_{s^{\text{se}},r,H}^{\text{intra}}-{\sigma }_{s^{\text{se}},r}^{\text{intra,0}})$$

$${\sigma }_{s^{\text{se}}}^{\text{inter,0}}=S_{s^{\text{se}}}^{\text{0}}={\sigma }_{s^{\text{se}},D}^{\text{inter}}=\text{smashoperator}[]\sum _{r\in \mathcal{R}}{W}_{D,r}{\sigma }_{s^{\text{se}},r,H}^{\text{intra}}$$

$${\sigma }_{s^{\text{ST}},r,H}^{\text{intra}}={\sigma }_{s^{\text{ST}},r}^{\text{intra,0}}$$

$${\eta }_c^{\text{in}}{p}_{c,r,h}^{\text{in}}=\frac{p_{c,r,h}^{\text{out}}}{{\eta }_c^{\text{out}}}$$

$${\pi }_{a^{\text{inv}}}=P_{a^{\text{inv}}}(U_{a^{\text{inv}}}^{\text{0}}+i_{a^{\text{inv}}})$$

$${\pi }_{a^{\text{not}}}=P_{a^{\text{not}}}{U}_{a^{\text{not}}}^{\text{0}}$$

$$p_{a,r,h}^{\text{out}}\le A_{a,r,h}{\pi }_a$$

$$p_{s,r,h}^{\text{in}}\le {\pi }_s$$

$$-R_g{\pi }_g\tau \le p_{g,r,h}^{\text{out}}-p_{g,r,h-1}^{\text{out}}\le R_g{\pi }_g\tau$$

$$-R_g{\pi }_g\tau \le \text{smashoperator}[]\sum _{r\in \mathcal{R}}(W_{d,r}{p}_{g,r,1}^{\text{out}}-W_{d-1,r}{p}_{g,r,H}^{\text{out}})\le R_g{\pi }_g\tau$$

$${\sigma }_{s,r,h}^{\text{intra}}\le S_s$$

$$S_{s^{\text{se}},d}^{\text{min}}{S}_{s^{\text{se}}}\le {\sigma }_{s^{\text{se}},d}^{\text{inter}}\le S_{s^{\text{se}},d}^{\text{max}}{S}_{s^{\text{se}}}$$

$$-L_l^{\text{imp}}\le f_{l,r,h}\le L_l^{\text{exp}}$$

$$\sum _{r\in \mathcal{R}}{W}_{d,r}=1,W_{d,r}\ge 0,\forall r\in \mathcal{R},d\in \mathcal{D}.$$

$$\sum _{r\in \mathcal{R}}{W}_{d,r}\le 1,W_{d,r}\ge 0,\forall r\in \mathcal{R},d\in \mathcal{D}.$$

$$\sum _i{a}_i{y}_{i,d}$$

$$=\sum _{r\in \mathcal{R}}{W}_{d,r}\underset{\le b}{\underbrace{\sum _i{a}_i{y}_{i,r}}}\le b\underset{\le 1}{\underbrace{\sum _{r\in \mathcal{R}}{W}_{d,r}}}\le b.$$

$$W_{d,r}\ge 0,\forall r\in \mathcal{R},d\in \mathcal{D}.$$

$${\mathbb{W}}^{\text{Dirac}}\subset {\mathbb{W}}^{\text{conv}}\subset {\mathbb{W}}^{\le 1}\subset \mathbb{W},$$

$$\left.\begin{aligned} \mathbf{D}&=[D_{n,x,d,h}]_{(n,x,h)\in\mathcal{N}\times\mathcal{X}\times\mathcal{H}}^{d\in\mathcal{D}},\\ \mathbf{A}&=[A_{g,d,h}]_{(g,h)\in\mathcal{G}\times\mathcal{H}}^{d\in\mathcal{D}},\\ \mathbf{E}&=[E_{s,d,h}]_{(s,h)\in\mathcal{S}\times\mathcal{H}}^{d\in\mathcal{D}},\\ \end{aligned}\quad\right\}\,\longrightarrow\quad\begin{aligned} \mathbf{C}=\left[\begin{array}[]{c}\mathbf{D}\\ \hline\cr\mathbf{A}\\ \hline\cr\mathbf{E}\end{array}\right].\end{aligned}$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle}{\mathrm{subject~to}}$}][l]{\underset{\displaystyle}{\mathrm{minimize}}}\quad\tfrac{1}{2}\lVert\mathbf{R}\mathbf{w}_{\!d}-\mathbf{c}_{d}\rVert_{2}^{2}$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle\phantom{}}{\mathrm{subject~to}}$}][c]{{\mathrm{subject~to}}}\quad$$

$$g({\mathbf{w}}_d)={\nabla }_{{\mathbf{w}}_d}\frac{1}{2}\Vert \mathbf{R}{\mathbf{w}}_d-{\mathbf{c}}_d{\Vert }_2^2={\mathbf{R}}^{\top }(\mathbf{R}{\mathbf{w}}_d-{\mathbf{c}}_d).$$

$${\Delta }^{\text{cur}}$$

$${\mathbf{w}}_d$$

Full text

1 Introduction

Energy System Optimization Models (ESOMs) are widely used to plan resilient and cost-effective energy systems. Common examples in the power sector include generation expansion planning (GEP) [1], integrated resource planning [2], mid-term operational planning [3], and transmission expansion planning [4]. As renewable energy sources like wind and solar gain traction, the variability and uncertainty in generation have increased—and will continue to rise with further integration [5]. Accurately capturing this variability requires high-resolution temporal modeling of demand and production. However, incorporating such fine-grained data over long planning horizons greatly increases the computational burden, making ESOMs difficult—or even intractable—to solve for large-scale systems [6].

1.1 Motivation

To address this challenge, researchers and practitioners often apply reduction techniques to temporal data. Instead of using all time series, they work with smaller subsets of time-varying data such as demand and renewable availability profiles. When periods are sufficiently similar, this can significantly reduce runtime without introducing substantial errors. Common approaches include time blocks and time slices [7], as well as representative periods (RPs) [8]. Recently, the latter has gained popularity, as it keeps the high temporal resolution while preserving the chronological structure within RPs, allowing for the modeling of seasonal storage, inter-period ramping, and other temporal constraints [9, 10].

The concept of RPs refers to specific time intervals, such as days or weeks, selected to capture the variability and characteristics of the ESOM over an extended period, such as a year or even decades. Selecting RPs involves a trade-off between model accuracy and computational efficiency. Therefore, identifying critical periods that accurately capture system behavior is an important research topic, as these periods influence how well the system is represented while still enabling accurate analysis [11].

1.2 Literature Review

Clustering algorithms are commonly utilized to obtain RPs based on time series of demand and renewable availability profiles [12]. Traditional clustering methods, such as $k$-medoids [9] and hierarchical clustering [13], are among the most common and straightforward approaches in the literature. However, these methods may not adequately capture extreme events or rare occurrences crucial for energy systems investment and operational planning decisions [14].

To address this limitation, researchers have proposed assigning greater weight to extreme values in the dataset [15]. This ensures that boundary conditions and rare events are captured in the reduced set of RPs. These weights are typically defined a priori, either through heuristic rules or expert knowledge of the system. Other approaches are data-driven. For example, [16] introduces a framework that selects RPs based on investment decisions rather than input data. However, this method requires solving multiple reduced model instances, which can be prohibitive—particularly for large-scale ESOMs with seasonal storage. Alternatively, hybrid clustering techniques, such as combining $k$-means with hierarchical clustering, have also been explored to improve the selection of RPs [17].

The concept of RP mixing has emerged as a means to recover information from the non-represented periods as a linear combination of the RPs. This concept was studied in [18] to better capture the diversity of original periods to reflect a broader range of conditions in the RP version of the ESOM. Nevertheless, the proposal in [18] overlooked the impact of the extreme conditions in the representative selection and proposed a mixed-integer programming optimization to find the RPs, which can be computationally intensive for large-scale problems. Additionally, that approach does not take into account the fact that the RPs will be mixed with each other in the expansion planning problem.

In conclusion, challenges persist in developing temporal aggregation methods that can effectively reduce computational requirements without significant loss of accuracy or that depend either on specific predefined rules or previous solutions of the reduced versions of the ESOM. The need for approaches that can capture both typical and extreme system conditions computationally efficiently is increasingly recognized.

1.3 Contributions

In response to the challenges shown in the previous section, this paper proposes the hull clustering with blended representative periods that aims to capture extreme conditions in the input data with fewer representatives than the counterparts, efficiently achieving high levels of accuracy, while maintaining a high fidelity when representing the non-extreme periods. The fundamental elements of the hull clustering with blended RPs are as follows.

Selection of Extreme Points: Instead of using cluster centroids or medoids, blended RPs select representatives from the data hull, focusing on extreme data points. This ensures that critical scenarios, which are often constraint-binding, are adequately represented, leading to optimization outcomes with lower regret. 2. 2.

Conic Sums for Blending: By using conic sums of representative periods, blended RPs create weighted combinations that more accurately reflect the data in the base periods. This approach reduces the number of RPs needed while preserving the integrity of the original data distribution, thus enhancing computational efficiency without compromising accuracy.

The positive effects of the contributions of this paper include:

Improved Solution Quality: The use of extreme points ensures that the selected RPs better capture the variability of the data, leading to more accurate solutions. 2. 2.

Reduction in Computational Burden: By reducing the number of RPs needed, the method lowers computational requirements, making longer planning horizons and more detailed models feasible.

However, the adoption of blended RPs introduces complexities, as the need for extreme points challenges conventional clustering techniques and demands the development of specialized algorithms to implement the hull clustering effectively.

In addition, we conduct two case studies, with data from the European network. The first one is a stylized version of a generation expansion planning problem, considering ramping but not storage. The second is a Power-to-X dispatch case study based on the global ambition scenario of the TYNDP for the year 2050 [19], which includes seasonal storage for hydro reservoir technologies. We demonstrate that our method achieves the same solution quality as conventional methods, using substantially fewer RPs and therefore runtime.

2 Model Formulation

Although the hull clustering with the blended RP method can be applied to a wide range of ESOMs, this paper focuses on two versions: a generation expansion planning (GEP) and a sector-coupling power and hydrogen dispatch problem (Power-to-X, P2X). Both versions include typical complexities in energy systems, such as transmission capacities, energy conversion, short-term and seasonal storage technologies, and ramping constraints, to illustrate the proposed method’s core ideas. In this section, we present the general formulation of an ESOM that can be adapted for both GEP and P2X optimization models. The GEP model aims to determine the optimal investment in new generation capacity while minimizing the sum of investment and operational costs, subject to changing demand, ramping, and resource availability. The P2X dispatch model aims to determine the optimal dispatch of existing capacities while minimizing operational costs, subject to available resource capacities, balance, conversion, and storage constraints. Given the complexity and temporal resolution of these problems, reducing the time series data through RPs is a common approach. Here, we outline the ESOM formulation and how RP methods are used to make large-scale problems computationally tractable.

2.1 Energy System Optimization Model

In energy system optimization, we consider a set of assets $a\in\mathcal{A}$. Each asset is located at a node $n\in\mathcal{N}$ in a network, which is connected by transmission lines $l\in\mathcal{L}$ transporting energy carriers $x\in\mathcal{X}$ (e.g., electricity or gas). We use $\mathcal{A}^{\text{in}}_{n,x}$ and $\mathcal{A}^{\text{out}}_{n,x}$ for assets that respectively consume and output carrier $x$ at node $n$. These assets are categorized into three functional types:

generation assets $g\in\mathcal{G}$, which produce energy (e.g., wind turbines, solar panels); 2. 2.

storage assets $s\in\mathcal{S}$, which store excess energy for later use (e.g., batteries, pumped hydro), divided into short-term $s^{\text{ST}}\in\mathcal{S}^{\text{ST}}$ and seasonal $s^{\text{se}}\in\mathcal{S}^{\text{seas}}$ storage; and 3. 3.

conversion assets $c\in\mathcal{C}$, which transform energy between carriers (e.g., fuel cells).

The asset types are disjoint:

[TABLE]

Some assets $a^{\text{inv}}\in\mathcal{A}^{\text{inv}}\subseteq A$ are designated as investable, meaning their capacities can be expanded as part of the planning process, subject to investment costs and budget constraints. Non-investable assets $a^{\text{not}}\in\mathcal{A}^{\text{not}}=A\setminus\mathcal{A}^{\text{inv}}$, by contrast, represent existing infrastructure. In addition, we define $\mathcal{L}^{\text{out}}_{n,x}$ and $\mathcal{L}^{\text{in}}_{n,x}$ to indicate the subset of outgoing and incoming lines to the node and carrier, respectively.

To capture the system’s temporal dynamics, optimization is performed over a horizon $\lvert\mathcal{T}\rvert=\lvert\mathcal{D}\rvert\cdot\lvert\mathcal{H}\rvert$, which is discretized into a set of $\lvert\mathcal{D}\rvert$ time periods $d\in\mathcal{D}$ (e.g., days) of equal length $\lvert\mathcal{H}\rvert$, each consisting of time steps $h\in\mathcal{H}$ (e.g., hours).

The temporal dimension of the model can be reduced by replacing base periods $d\in\mathcal{D}$ with a much smaller set of RPs $r\in\mathcal{R}$, where $\lvert\mathcal{R}\rvert\ll\lvert\mathcal{D}\rvert$. The RPs are selected so that they reflect typical combinations of the relevant conditions, reducing the computational burden while preserving the essential temporal patterns needed for accurate decision-making.

RPs are typically identified using clustering methods. While $k$-means is commonly used, alternative methods such as $k$-medoids and hierarchical clustering may offer better performance depending on the structure of the data. For instance, $k$-medoids selects actual observations, making it more robust to outliers. After clustering, each base period $d$ is assigned an RP $d\mapsto\rho(d)$ by a clustering function $\rho:\mathcal{D}\to\mathcal{R}$. As a result, they are replaced with their RPs $r=\rho(d)$. We introduce weights $W_{\!d,r}$ as indicators of this assignment:

[TABLE]

Additionally, RPs $r$ are assigned total weights $W_{\!r}$ representing the number of periods that are replaced by this representative: After clustering, each base period $d$ is assigned a representative $r$ by a clustering function $\rho:\mathcal{D}\to\mathcal{R}$. In the model, each base period $d$ is then replaced by its representative $r=\rho(d)$. RPs $r$ are assigned weights $W_{\!r}=\sum_{d\in\mathcal{D}}W_{\!d,r}$ representing the number of periods that are replaced by this representative.

Using this approach, we propose the following ESOM that integrates investment and operational decisions. All input parameters are denoted by capital letters and decision variables by lower case letters. Except for the flows $f_{l,r,h}$, the decision variables are constrained to be non-negative, but these constraints are omitted for brevity.

[TABLE]

The model (2) minimizes the total cost (2a) comprising two main components. First, the total investment cost $c^{\text{inv}}$ in (2b) is computed as the sum of investment costs across all investable assets $a\in\mathcal{A}^{\text{inv}}$. For each asset, the cost is determined by its capacity per unit $P_{a}$, the annualized investment cost per MW $I_{a}$, and the number of units selected for investment $i_{a}$. Second, the total operational cost $c^{\text{op}}$ in (2c) accounts for generation and reliability-related costs across all RPs $r\in\mathcal{R}$ and time steps $h\in\mathcal{H}$. These include costs associated with generation output $p^{\text{out}}_{g,r,h}$, as well as penalties for spillage $e^{\text{sp}}_{s^{\text{se}},r,h}$ and borrowing $e^{\text{br}}_{s^{\text{se}},r,h}$ of water. Spillage and borrowing are slack variables for water reservoir level constraints; for non-reservoir technologies, these are set to zero. To ensure consistency with annualized investment cost, the operational cost is also annualized using a weight $W^{\text{op}}=\lvert\mathcal{Y}\rvert/\lvert\mathcal{T}\rvert$, where $\lvert\mathcal{Y}\rvert$ and $\lvert\mathcal{T}\rvert$ are the numbers of time steps in a year and in the model, and $\tau$ is the duration of the time step.

The next group of constraints enforces intra-period energy system balances. The node balance constraint (2d) ensures that the local demand $D_{n,x,r,h}D^{\text{max}}_{n,x}$ is satisfied by the sum of locally generated, consumed, and net imported power (imports minus exports). Here $D^{\text{max}}_{n,x}$ represents the peak demands, and $D_{n,x,r,h}\in[0,1]$ are the time-varying demand profiles. The intra-period storage balance constraints (2e)–(2h) track the state of charge $\sigma^{\text{intra}}_{s,r,h}$ of storage assets across time, accounting for the previous state of charge, initial storage level $\sigma^{\text{intra,\,0}}_{s,r}$, charging and discharging flows adjusted by efficiencies $\eta^{\text{in}}_{s}$ and $\eta^{\text{out}}_{s}$, spillage $e^{\text{sp}}_{s,r,h}$ and borrowing $e^{\text{br}}_{s,r,h}$, and external inflows $E_{s,r,h}E^{\text{max}}_{s,r,h}$, where $E_{s,r,h}\in[0,1]$, commonly representing hydro inflows for reservoirs. Some technologies may not include either internal $\eta^{\text{in}}_{s}p^{\text{in}}_{s,r,h}\tau$ or external inflows, in which case either $p^{\text{in}}_{s,r,h}$ or $E^{\text{max}}_{s,r,h}$ are set to zero.

The inter-period storage balance constraints (2i)–(2j) recover the chronological information for the seasonal storage when using representative periods by linking the inter- and intra-period storage levels through the use of the weights $W_{\!d,r}$ (note that the weights are not necessarily integer).

The cyclic constraints (2k) link the last inter-period storage level $\sigma^{\text{inter}}_{s^{\text{se}},D}$ to the initial one $\sigma^{\text{inter,\,0}}_{s^{\text{se}}}$ and set both to a pre-defined value $S^{\text{0}}_{s^{\text{se}}}$. This is done to avoid creating a solution where the storage is fully charged initially and discharges over the time horizon, creating free energy in the system. Intra-period cycling constraints (2l) serve the same role but for short-term storage. Additionally, in the model formulation so far, the intra-period storage variables $\sigma^{\text{intra}}_{s,r,h}$ are defined up to a constant: increasing all or them by the same amount does not change the inter-period variables $\sigma^{\text{inter}}_{s,d}$. To address this issue, we “tether” them to each other via the last part of (2k).

The conversion balance constraint (2m) ensures that energy inputs and outputs at conversion assets remain consistent with their conversion efficiencies $\eta^{\text{in}}_{c}$ and $\eta^{\text{out}}_{c}$.

The next set of constraints captures the physical limitations of the assets. Constraints (2n) and (2o) define the total number of available units $\pi_{a}$ as the sum of existing units $U^{\text{0}}_{a}$ and newly invested units $i_{a}$. Operational capacity constraints then limit each asset’s activity based on this availability. Specifically, constraints (2p) and (2q) ensure that the produced power $p^{\text{out}}_{a,r,h}$ and consumed power $p^{\text{in}}_{a,r,h}$ do not exceed the time-dependent availability $A_{a,r,h}$. This parameter equals 1 for all assets except renewable generators, for which it reflects hourly availability with $A_{g,r,h}\in[0,1]$. Ramping constraints (2r) and (2s) limit the change in generation between consecutive time steps by the ramping parameter $R_{g}$, modeling real-world operational inflexibility. Constraints (2t) and (2u) ensure that storage levels remain withing required upper and lower limits, while transmission constraints (2v) enforce that line flows $f_{l,r,h}$ stay within import $-L^{\text{\!imp}}_{l}$ and export $L^{\text{\!exp}}_{l}$ limits.

In this paper, we focus on this specific problem formulation for clarity of presentation. When the investment decisions $i_{a}$ are treated as decision variables and restricted to generation assets—i.e., $\mathcal{A}^{\text{inv}}\subseteq\mathcal{G}$—the model (2) reduces to the GEP problem. Alternatively, if all investment decisions $i_{a}$ are fixed as input parameters, the model can be used for operations planning such as P2X, focusing solely on the cost-efficient dispatch of existing generation assets. The formulation is flexible and can be used for storage investment, and with trivial adjustment extended to transmission expansion planning.

3 Methodology

We have identified that any representative period method consists of two key components: a clustering method for finding RPs, and the mapping of base-period data to RPs. In this paper, we introduce a generalized framework that allows for flexible choices in both components. This leads to a new variant, which we term the method of blended RPs. The resulting procedure is summarized in Algorithm 1.

3.1 Overview of Blending in Representative Periods

The first key innovation of the method of blended representative periods lies in generalizing how base periods are represented in terms of a reduced set of RPs. Traditional clustering-based approaches assign each base period to exactly one RP. In such a framework, the clustering function $\rho:\mathcal{D}\to\mathcal{R}$ defines a partition of the original periods: each base period $d\in\mathcal{D}$ is mapped to a single RP $r\in\mathcal{R}$. This assignment, however, can be limiting, as in many practical settings—such as renewable generation—the base periods exhibit patterns that are mixtures of different underlying conditions. For instance, a mildly sunny, breezy afternoon may not be well captured by any single RP but may lie somewhere between a typical sunny day and a typical windy day. In such cases, forcing a one-to-one mapping can introduce large approximation errors.

To address this, we introduce the idea of blended representative periods, where each base period is approximated by a weighted combination of multiple RPs. This blend is captured by a weight matrix111In matrix definitions, we use subscript and superscript to describe row and column indices respectively. $\mathbf{W}=[W_{d,r}]_{d\in\mathcal{D}}^{r\in\mathcal{R}},$ where each row specifies how the base period $d$ is expressed in terms of RPs $r$. Instead of mapping each base period $d$ to a single RP $r$, we can now represent the problem data $z$ as $y_{d}\approx\sum_{r\in\mathcal{R}}W_{\!d,r}\cdot y_{r},$ where $y_{d}$ stands for a period-dependent parameter or variable, such as demand $D_{n,x,d,h},$ or produced power $p^{\text{out}}_{a,d,h}$. This generalization allows us to interpolate between RPs and better preserve the diversity and nuance of the original time series.

3.2 Blending Types: Dirac, Convex, and (Sub-Unit) Conical

The choice of weight matrix $\mathbf{W}$ is central to how base periods are approximated in the blended RP method. We have identified and analyzed four types of blending schemes, based on the geometric structure imposed on the weights: Dirac, convex, sub-unit conical, and conical. These represent increasing levels of flexibility in how base periods can relate to the representative periods.

3.2.1 Dirac Weights (‘Hard’ Clustering)

This is the traditional approach discussed in Section 2.1. Each base period $d\in\mathcal{D}$ is assigned to exactly one RP $r\in\mathcal{R}$, and the weight matrix $\mathbf{W}$ is binary with a single nonzero entry per row as per (1). This structure corresponds to a Dirac measure in discrete space (i.e., a measure concentrated on a single point). While computationally efficient, Dirac blending cannot capture overlap between RPs and may introduce abrupt transitions.

3.2.2 Convex Weights

A more flexible alternative allows each base period to be written as a convex combination of RPs:

[TABLE]

Convex blending allows smooth interpolation between RPs; as such, it has already been used in the context of ESOMs [18].

3.2.3 Sub-Unit Conic Weights

This more expressive variant allows the weights to be non-negative and sum at most into 1:

[TABLE]

This ensures that inequality constraints remain valid without modification. To see why, consider an arbitrary linear inequality involving variables $y_{i,r}$ associated with RPs, $\sum_{i}a_{i}y_{i,r}\leq b,\quad\forall r\in\mathcal{R},$ where $a_{i}$ are arbitrary constraint coefficients and $b$ is non-negative. This implies that $\sum_{i}a_{i}y_{i,d}\leq b$ in the base problem holds as well:

[TABLE]

By the same argument, inequalities of the form $-b\leq\sum_{i}a_{i}y_{i,d}$ hold for non-negative constants $b$. Thus, sub-unit conic blending naturally preserves the feasibility of inequality constraints in the blended problem, without any additional reformulation.

3.2.4 Conic Weights

The most general linear blending we consider imposes only non-negativity:

[TABLE]

This allows each base period to be approximated by any non-negative linear combination of RPs, including combinations with total weight greater than one. Such conic blending is the most expressive but may be less interpretable or lead to infeasibility of the solution to the original problem with respect to upper bounds, as the reasoning used in (3) no longer holds.

As far as we are aware, the (sub-unit) conic weights have not been described in the literature in the context of RPs.

Geometrically, each type of blending weight corresponds to projecting base period data onto a different type of hull formed by the RPs. This results in projection errors, which are the differences between the original base period data and its reconstruction from a weighted combination of representative periods, and it decreases as the space of admissible weights becomes more expressive. Specifically, Dirac weights select a single RP corresponding to the vertices of the set. Convex weights project each base period into the convex hull of the RPs; data points inside the hull can be approximated without projection error. Sub-unit conic weights extend the hull by adding the null point into it, reducing projection errors even further. Finally, general conic weights project base periods into the full conic hull of the RPs.

Figure 1 illustrates this visually. It shows how the projection error shrinks as we move from approximating a base period with a single RP (Dirac) to the full conic hull. The more expressive the hull, the more accurately we can approximate the data that is being clustered. This hierarchy reflects a trade-off between expressiveness and constraint preservation.

We use the term blended weights as a general descriptor encompassing all four weight types introduced here: Dirac, convex, sub-unit conic, and conic weights. These sets of weights satisfy the strict inclusions

[TABLE]

where $\mathbb{W}^{\,\cdot\,}$ is the space of all weight matrices $\mathbf{W}$ of a type $(\,\cdot\,)$ and $\mathbb{W}$ is the space of conic weights. Thus, each subsequent class further generalizes the conventional Dirac weights.

Some of these spaces have well-known geometric interpretations and standard names. For instance, the space of all conic weights corresponds to the positive orthant $\mathbb{R}_{+}^{n}$, $\mathbb{W}=\mathbb{R}_{+}^{n}$, the space of convex weights is the standard $n$-simplex $\bigtriangleup_{n}$, $\mathbb{W}^{\text{conv}}=\bigtriangleup_{n}=\bigl\{\mathbf{v}\in\mathbb{R}^{n}_{+}\bigm|\lVert\mathbf{v}\rVert_{1}=1\bigr\}$, and the space of Dirac weights is the standard $n$-dimensional basis $\mathbb{W}^{\text{Dirac}}=\{\mathbf{e}_{i}\}_{i=1,\dots,n}$.

Besides providing a unifying framework that generalizes existing RP methods, we propose to blend these periods using weights that are not necessarily convex. This allows our method to go beyond simple interpolation between RPs and instead express the base period data with greater flexibility and nuance.

3.3 Constructing the Clustering Matrix

The first step of Algorithm 1 is to represent each base period $d\in\mathcal{D}$ by a feature vector that captures its relevant time-varying characteristics. Collecting these vectors forms the clustering matrix $\mathbf{C}$. This matrix is required because clustering algorithms such as $k$-medoids operate in a feature space: to meaningfully group periods, each period must be described by a vector that reflects its time-varying data, rather than by its index alone. For ESOMs, these typically include demands $D_{n,x,d,h}$, availability profiles $A_{g,d,h}$ of renewable generation technologies, and time-varying energy storage inflows $E_{s,d,h}$.

Formally, we define $\mathbf{C}=[c_{k,d}]_{k\in\mathcal{K}}^{d\in\mathcal{D}}$ where each column $\mathbf{c}_{d}$ encodes the time series of all features for period $d$. To obtain these feature vectors, we first compute the data matrices for each feature type and then stack them vertically:

[TABLE]

To ensure that all features contribute comparably to the clustering process, we use normalized data; that is, the problem is formulated such that all values used in clustering lie between zero and one. Once $\mathbf{C}$ is constructed, it is passed to the clustering algorithm of choice to find the RPs.

3.4 Hull Clustering

Conventional clustering algorithms such as $k$-means and $k$-medoids are typically well-suited for Dirac-weight approximations, where each base period is assigned to a single RP. These methods minimize average point-to-cluster-center distances but do not consider whether the selected cluster centers (RPs) span a suitable hull for blending.

When using blended weights, it becomes essential that each base period can be reconstructed as a weighted combination of RPs, where the weights satisfy specific constraints. If the selected RPs do not span the correct hull, base periods will be projected outside of the feasible blending space, resulting in significant projection error.

3.4.1 Hull Types

Given a set of RPs, the choice of weight type determines the kind of hull within which the original periods must lie in order to be exactly reconstructible. Specifically:

• •

Convex hull: For convex weights, the original periods must lie within the convex hull of the representative periods. This is the smallest convex set containing them.

• •

Convex hull with null: For sub-unit conic weights, reconstruction is possible if the original periods lie within the convex hull of the representative periods and the origin (null). Intuitively, this corresponds to finding convex weights that may assign some mass to a null (zero) vector. After discarding the null, the remaining positive weights sum to less than one, hence sub-unit.

• •

Conic hull: For general conic weights, the reconstruction is valid as long as the original periods lie within the conic hull of the RPs.

We illustrate the differences between these hulls in Figure 2, where each successive hull type expands the space in which the original data points can be expressed using the RPs without introducing a projection error.

3.4.2 Greedy Hull Clustering

While the concept of using a convex or conic hull to minimize the projection error is natural, identifying such a hull in practice is nontrivial. Specifically, selecting a small subset of representative periods whose hull contains all the original periods is a computationally hard problem. This is related to well-known non-negative matrix factorization problem, which is NP-hard [20].

To address this challenge efficiently, we adapt the greedy conic hull algorithm [21, Algorithm 1] for all three hull types. The algorithm identifies points in the conic hull via a reduction to the convex hull with the gnomonic projection (i.e., a projection of rays onto points in a plane); omitting the projection yields a direct convex-hull method.

We additionally improve the algorithm by employing caching to avoid redundant distance computations during the greedy selection of the hull points. Specifically, once a point’s distance to the current convex hull has been computed, that value is stored in a cache. Since the convex hull can only grow as new points are added, the distance from any given point to the hull cannot increase over time. This property allows the algorithm to safely reuse previously computed distances unless there is strong evidence that a point’s position relative to the hull has changed significantly—namely, if its distance to the newly added point is smaller than its cached hull distance. By using this check, the algorithm reduces the number of expensive projection operations, thereby improving overall efficiency. Although this criterion is a heuristic, we nevertheless observed that in practice it offers a substantial speedup. The resulting method is presented in Algorithm 2.

Depending on the hull type, the full clustering procedure is:

• •

Convex hull: Algorithm 2 is applied directly to $\mathbf{C}$.

• •

Convex hull with null:

add a null vector $\mathbf{0}$ as a column to $\mathbf{C}$, 2. 2.

run Algorithm 2 with $\lvert\mathcal{R}\rvert+1$ representatives and initial set of RPs $\mathcal{R}=\{\mathbf{0}\}$, then 3. 3.

drop the null vector $\mathbf{0}$ from both $\mathbf{C}$ and $\mathcal{R}$.

• •

Conic hull: use the gnomonic projection [21]:

compute the column mean $\bar{\mathbf{c}}$ of the data and normalize it to obtain a reference direction $\mathbf{q}$, 2. 2.

for each column $\mathbf{c}_{d}$ of the clustering matrix $\mathbf{C}$, scale it by $1/\langle\mathbf{c}_{d},\mathbf{q}\rangle$ so that all points lie on the hyperplane $\mathcal{H}=\bigl\{\mathbf{x}\bigm|\langle\mathbf{x},\mathbf{q}\rangle=1\bigr\}$, 3. 3.

apply Algorithm 2 to the scaled clustering matrix, 4. 4.

return the corresponding unscaled RPs.

The greedy procedure ensures that each added representative improves coverage of the dataset, and it can be terminated once a fixed number of representatives is reached.

At each iteration, the method aims to select the data point that most reduces the projection error of the remaining points onto the current hull. This approach allows us to build up the RP set incrementally while maintaining a good approximation of the original data within the chosen hull type.

Computing the distance from a data point $\mathbf{c}_{d}$ to the current convex hull $\operatorname{Hull}\{\mathbf{c}_{r}\mid r\in\mathcal{R}\}$, as required in line 2 of Algorithm 2, is a nontrivial task. We formulate this as a projection problem: given a point $\mathbf{c}_{d}$, we seek the closest point in the convex hull in terms of Euclidean distance $\operatorname{dist}(\mathbf{R}\mathbf{w}_{\!d},\mathbf{c}_{d})=\lVert\mathbf{R}\mathbf{w}_{\!d}-\mathbf{c}_{d}\rVert_{2}$. Since any point in the hull of the columns of $\mathbf{R}$ can be expressed as a convex combination of those columns using weights $\mathbf{w}_{\!d}$, that is, $\mathbf{R}\mathbf{w}_{\!d}$, this projection reduces to the following constrained optimization problem:

[TABLE]

Without the constraint (4b), the problem has a well-known solution $\mathbf{R}^{\dagger}\mathbf{c}_{d}$ where $\mathbf{R}^{\dagger}$ is the Moore–Penrose pseudoinverse of $\mathbf{R}$. This might result in negative weights, and even projecting back onto the subspace of acceptable weights $\operatorname{proj}_{\mathbb{W}^{\text{conv}}}(\mathbf{R}^{\dagger}\mathbf{c}_{d})$ yields a good feasible but not necessarily optimal solution.

To address this issue, we use projected gradient descent (PGD, described in Algorithm 3) to further fine-tune the initial weight guess. Note that replacing $\lVert\,\cdot\,\rVert_{2}$ with $\tfrac{1}{2}\|\cdot\|^{2}_{2}$ does not alter the location of the minimizer, since the function $f(x)=\tfrac{1}{2}x^{2}$ is strictly increasing on $\mathbb{R}_{+}$. Thus, minimizing the squared distance is equivalent to minimizing the distance itself, while offering a simpler subgradient $g$ of the objective:

[TABLE]

This subgradient is smooth, making the projection step more amenable to efficient optimization using PGD, which we employ to compute the distance to the hull in line 2 of Algorithm 2 as follows:

[TABLE]

We use $\operatorname{proj}_{\mathbb{W}^{\text{conv}}}(\mathbf{R}^{\dagger}\mathbf{c}_{d})$ as the initial weight guess and (5) as the subgradient function $g$. For the projection operator $\operatorname{proj}_{\mathbb{W}^{\text{conv}}}$, we utilize the fast algorithm [22, Figure 2] for projection onto the simplex $\bigtriangleup_{n}=\mathbb{W}^{\text{conv}}$. For the implementation details, the reader is referred to the original paper.

3.5 Finding the Blended Weights

After hull clustering is complete, we obtain a matrix $\mathbf{R}$ of the RP data. The next step is to find the weights $W_{\!d,r}$ that best approximate the original data for each period $d$.

This is the weight-fitting problem (4) but with the weights potentially restricted to a different domain $\mathbb{W}^{\,\cdot\,}$. It can be solved with Algorithm 3 using different projection operators $\operatorname{proj}_{\mathbb{W}^{\,\cdot\,}}$ depending on the weight type.

• •

Dirac weights $\mathbb{W}^{\text{Dirac}}$: project by setting the largest coordinate of $\mathbf{x}$ to 1 and all others to 0. In practice, cluster assignments directly provide this projection, so explicit projection is often unnecessary.

• •

Convex weights $\mathbb{W}^{\text{conv}}$: apply the aforementioned fast projection algorithm [22, Figure 2].

• •

Conic sub-unit weights $\mathbb{W}^{\leq 1}$: use the convex hull method with the null point trick as described previously.

• •

Conic weights $\mathbb{W}$: project by thresholding negative values to zero, $\operatorname{proj}_{\mathbb{W}}(\mathbf{x})=\mathbf{x}^{+}=\max(\mathbf{x},\mathbf{0}).$

For the initial weight guess, we use either $\operatorname{proj}_{\mathbb{W}^{\text{conv}}}(\mathbf{R}^{\dagger}\mathbf{c}_{d})$ or the Dirac weights returned by the clustering algorithms, whichever results in a smaller projection error.

4 Case Studies

We test the proposed methodology in two case studies.222The input data and implementation for both case studies are available online at https://github.com/greg-neustroev/blended-rep-periods. The first case study is a GEP problem based on information from the TYNDP 2022 [19]. This stylized case study allows for investments in renewable energy from greenfield setting. For the sake of simplicity, this case study includes ramping constraints but no storage and conversion technologies.

The second case study is also based on information from the TYNDP, but from the 2024 study; in this case, we include hydrogen demand for each node, which can be met through electrolyzers that convert electricity into hydrogen or through steam methane reforming. Unlike the GEP case study, this analysis excludes ramping constraints and investment decisions, with all capacities predetermined. It includes both short-term and seasonal storage options, such as batteries and pumped hydro storage, and models hydropower plants, reservoirs, and inflows in each country. This approach highlights the importance of seasonal behavior, particularly in the Nordic countries.

4.1 Experiments Setup

For each case study, we design experiments to evaluate the performance of different clustering methods and weight types.

• •

Weight types: we compare standard Dirac weights with the proposed conic, conic sub-unit, and convex weights.

• •

Clustering methods: in addition to the common approaches ($k$-means, $k$-medoids), we consider hull methods specific to each weight type (convex for Dirac).

• •

Number of representatives: we vary the number of RPs to analyze its impact on each method’s performance.

Each experiment is run with five different random seeds to account for various sources of randomness.

4.2 Regret as a Measure of Solution Quality

RP methods approximate full-resolution ESOMs—e.g., modeling a full year hourly—but it can be hard to judge how good these approximations are. Rather than just comparing objective values, we use relative regret, which measures the extra cost (in percent) incurred when decisions are made using the approximate model instead of the full one. In other words, regret shows how much worse the solution is because of the simplification, with lower regret values being better.

In the GEP case study, we compute the regret of the investment decisions, that is, the investment variables from the approximate model are fixed, and the full ESOM is rerun to calculate the actual total cost (i.e., investment and operation). These costs are compared against the costs of the full original ESOM.

In the P2X case, seasonal storage levels (states of charge at the end of every day) are fixed, and rerunning the full problem gives the operational cost; regret is again the difference from the benchmark.

4.3 Results for Generation Expansion Planning

We run the GEP case studies for 5, 10, 20, 40, and 80 RPs.

Figure 3 shows the regret of all of the experiments grouped by the weight type. In case of both $k$-means and $k$-medoids, the weight type has almost no effect on the performance. In each case, the hull method outperforms both conventional clustering methods. Moreover, with just 5 RPs, conic clustering achieved relative regret of 7.4%, a level not achieved by any of the $k$-means results (lowest relative regret 26.1%) and matched only by $k$-medoids with 80 RPs (6.4%), as seen in Figure 4.

Figure 5 shows the total time, that is, the time required to read, preprocess, and cluster the data, fit the weights, formulate the model and solve it. Conic hull with 5 RPs took $5.3\pm 0.4$ s, compared to $54.0\pm 0.6$ s for 80-RP $k$-means; it was, therefore, ten times faster and achieved 3.5 times lower regret. $k$-medoids with 80 RPs took $75.1\pm 1.9$ s, being 14 times slower than conic hull method to achieve similar regret levels.

Blended weight methods, however, scale worse than Dirac methods: at 80 RPs they were slower than solving the full problem directly, as weight fitting dominates runtime. In contrast, Dirac weights with twice as many RPs often matched blended performance, highlighting the trade-off between weight sophistication and the number of RPs.

In summary, in GEP hull clustering—especially with conic weights—provided superior accuracy and efficiency, dominating both $k$-means and $k$-medoids.

4.4 Results for Power-to-X Dispatch with Seasonal Storage

As the P2X case study is more complex, 5 RPs were insufficient; therefore, experiments began with 10 RPs.

Figure 6 reports relative regret for all methods. Unlike in the GEP case, weight type plays a decisive role: blended weights outperform Dirac weights, regardless of clustering method. The runtime–regret trade-off (Figure 8) highlights two parts of the Pareto front: convex-weight hull clustering, which achieves regret below 15% with a runtime of 17 s, and Dirac $k$-medoids, which runs faster but with more than 50% regret.

Blended weights also capture reservoir dynamics better. Figure 7 shows seasonal storage levels for Norway, the largest reservoir in the system. From 40 RPs onward, hull clustering closely tracks the full model, whereas $k$-means diverges. With just 10 or 20 RPs, convex hull clustering yields near-optimal solutions (13.3% and 16.1% regret) compared to much higher regret levels for $k$-means (64.5% and 42.7%) and $k$-medoids (54.8% and 47.6%) with the same numbers of RPs.

Overall, in P2X the choice of weight type is more important, with convex hull methods striking the best balance between regret, runtime, and reservoir fidelity.

5 Discussion

The results highlight several benefits of using RP methods with advanced weighting and clustering strategies. Hull clustering with blended weights consistently delivers lower regret, faster runtimes, and more faithful approximations. These improvements are especially valuable when solving large-scale, high-resolution ESOMs, where runtime and memory requirements would otherwise be prohibitive.

At the same time, the case studies reveal that performance is strongly problem-dependent. In GEP with ramping, the main challenge is capturing chronological variability; what matters most is selecting RPs that preserve extreme ramps. Clustering type therefore dominates performance: hull methods excel because they select diverse boundary points that best capture these extremes. By contrast, in P2X, the binding constraints are less about ramps and more about seasonal storage balances. Once RPs are chosen, weights determine how they are “stitched back together” into a synthetic year. Here, the choice of weights is critical. Even though clustering type is more important for GEP and weight type for P2X, the combination of hull clustering and blended weights performs well in both cases; while not necessarily optimal, using them together appears to be important, as they exploit different problem structures.

The experiments highlight potential pitfalls. The problem-dependence means that applying a method blindly across models may lead to poor results. Moreover, blended weights, while powerful, introduce additional computational overhead; at the same time, we do not expect this to be an issue for large-scale problems, since the time to solve the optimization problem will very likely dominate the clustering time. For very large numbers of RPs, weight fitting becomes the bottleneck, sometimes exceeding the runtime of the full problem. This makes parameter tuning, such as the learning rate and step length in the gradient descent algorithm used for weight fitting, critical to achieving convergence without excessive iterations, calling for further research.

6 Conclusion

Existing clustering-based approaches using representative periods (RPs) to solve large ESOMs can be captured in a general framework. Within this framework, we introduced hull clustering with blended RPs, a novel approach to selecting and combining extreme and typical time series to efficiently reduce temporal data in ESOMs. By focusing on extreme points and blending them, the method captures critical system conditions more accurately than conventional clustering techniques, while requiring fewer RPs. Through two case studies—a stylized GEP problem and a more comprehensive P2X dispatch problem—we demonstrated that the approach achieves high solution quality with substantially reduced computational burden. The results highlight its ability to preserve important behavior patterns, improving fidelity, and achieving lower regret. While implementing blended RPs introduces new challenges in the clustering methodology, the approach offers a promising avenue for scaling ESOMs to longer horizons and more detailed models. Future work could extend the method to a more general formulation, test it on more diverse case studies, and integrate it directly into existing ESOM tools.

Acknowledgments

This publication is part of the project NextGenOpt with project number ESI.2019.008, which is financed by the Dutch Research Council (NWO) and supported by eScienceCenter under project number NLeSC C 21.0226. In addition, this research received partial funding from the European Climate, Infrastructure and Environment Executive Agency under the European Union’s HORIZON Research and Innovation Actions under grant agreement № 101095998.