# A Scenario-Spatial Decomposition Approach With a Performance Guarantee for the Combined Bidding of Cascaded Hydropower and Renewables

Luca Santosuosso, Simon Camal, Arthur Lett, Guillaume Bontron, Jalal Kazempour, Georges Kariniotakis

<figure><img src="image_1.png" /></figure>

## To cite this version:

Luca Santosuosso, Simon Camal, Arthur Lett, Guillaume Bontron, Jalal Kazempour, et al. A Scenario-Spatial Decomposition Approach With a Performance Guarantee for the Combined Bidding of Cascaded Hydropower and Renewables. IEEE Transactions on Power Systems, 2025, pp.1-12. 10.1109/TPWRS.2025.3638005. hal-04738929v2

HAL Id: hal-04738929

https://minesparis-psl.hal.science/hal-04738929v2

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Scenario-Spatial Decomposition Approach With a Performance Guarantee for the Combined Bidding of Cascaded Hydropower and Renewables

Luca Santosuosso, *Graduate Student Member, IEEE*, Simon Camal, *Member, IEEE*, Arthur Lett, Guillaume Bontron, Jalal Kazempour, *Senoir Member, IEEE*, and Georges Kariniotakis, *Senoir Member, IEEE*

**Abstract**—This study develops a scalable co-optimization strategy for the joint bidding of cascaded hydropower, wind, and solar energy units, treated as a unified entity in the day-ahead market. Although hydropower flexibility can manage the stochasticity of renewable energy, the underlying bidding problem is complex due to intricate coupling constraints and nonlinear dynamics. A decomposition in both scenario and spatial dimensions is proposed, enabling the use of distributed optimization. The proposed distributed algorithm is eventually a heuristic due to non-convexities arising from the system's physical dynamics. To ensure a performance guarantee, trustworthy upper and lower bounds on the global optimum are derived, and a mathematical proof is provided to demonstrate their existence and validity. This approach reduces the average runtime by up to 35% compared to alternative distributed methods and by 57% compared to the centralized optimization. Moreover, it consistently delivers solutions, whereas both centralized and alternative distributed approaches fail as the size of the optimization problem grows.

**Index Terms**—Trustworthy bounds, distributed optimization, stochastic programming, cascaded hydropower, bidding strategy.

## NOMENCLATURE

### Indices and Sets

* $\mathcal{N}$ Set of hydropower plants, indexed by $n$

* $S_n^{\text{oc}}$ Set of inflow segments for the piecewise linear approximation of asset $n$'s operational curve, indexed by $i$

* $T$ Set of time steps, indexed by $t$

* $\Omega$ Set of scenarios, indexed by $\omega$

* $k$ Index for iteration of the distributed algorithm

### Physical Constants

* $g$ Acceleration due to gravity (m/s$^2$)

* $\omega$ Water density (kg/m$^3$)

### Modeling Parameters

* $\bar{H}_n$ Maximum hydraulic head of asset $n$ (m)

* $\underline{H}_n$ Minimum hydraulic head of asset $n$ (m)

* $P_\omega$ Probability associated with scenario $\omega$ (-)

* $\bar{P}_n^H$ Capacity of asset $n$ (MW)

* $\bar{P}_{\omega,t}^S$ Solar farm output in scenario $\omega$ at time $t$ (MW)

* $\bar{P}_{\omega,t}^W$ Wind farm output in scenario $\omega$ at time $t$ (MW)

L. Santosuosso, S. Camal and G. Kariniotakis are with the PERSEE Center of MINES Paris PSL University, 06904 Sophia-Antipolis, France (e-mail: luca.santosuosso@minesparis.psl.eu; simon.camal@minesparis.psl.eu; georges.kariniotakis@minesparis.psl.eu).

J. Kazempour is with the Technical University of Denmark, 2800 Kgs. Lyngby, Denmark (e-mail: jalal@dtu.dk).

G. Bontron and A. Lett are with the Compagnie Nationale du Rhône, 69004 Lyon, France (e-mail: g.bontron@cnr.tn.fr; a.lett@cnr.tn.fr).

* $\bar{Q}_n^{\text{br}}$ Minimum barrage discharge of asset $n$ (m$^3$/s)

* $\bar{Q}_{n,\omega,t}^{\text{ext}}$ Forecast of external inflow (from river and its tributaries) to asset $n$ in scenario $\omega$ at time $t$ (m$^3$/s)

* $\bar{Q}_{n,i}^{\text{in}}$ i-th inflow segment of asset $n$ (m$^3$/s)

* $\bar{Q}_n^{\text{tr}}$ Maximum turbines discharge of asset $n$ (m$^3$/s)

* $\bar{Q}_n^{\text{tr}}$ Minimum turbines discharge of asset $n$ (m$^3$/s)

* $S_{n,\omega,t}$ n-th asset surface area in scenario $\omega$ at time $t$ (m$^2$)

* $Z_n^{\text{fl},0}$ Initial forebay water level of asset $n$ (m)

* $\bar{Z}_{n,i}^{\text{fl}}$ i-th maximum forebay segment of asset $n$ (m)

* $\bar{Z}_{n,i}^{\text{fl}}$ i-th minimum forebay segment of asset $n$ (m)

* $Z_n^{\text{tr}}$ Tailrace water level of asset $n$ (m)

* $\Delta_T$ Sampling period (s)

* $\Delta q_n^{\text{tr}}$ Turbines outflow ramp limit of asset $n$ (m$^3$/s)

* $\eta_n^{\text{tr}}$ Turbine efficiency of asset $n$ (-)

* $\hat{\pi}_{\omega,t}^{\text{E}}$ Forecast of day-ahead energy price in scenario $\omega$ at time $t$ (€/MWh)

* $\hat{\pi}_{\omega,t}^{\text{E}\uparrow}$ Forecast of positive imbalance penalty price in scenario $\omega$ at time $t$ (€/MWh)

* $\hat{\pi}_{\omega,t}^{\text{E}\downarrow}$ Forecast of negative imbalance penalty price in scenario $\omega$ at time $t$ (€/MWh)

* $\tau_{a,b}^{\text{br}}$ Barrage water travel time from asset $a$ to $b$ (s)

* $\tau_{a,b}^{\text{tr}}$ Turbined water travel time from asset $a$ to $b$ (s)

### Notation for the Proposed Distributed Algorithm

* $\bar{K}$ Maximum number of iterations allowed (-)

* $J$ Objective function of the bidding problem (€)

* $J^{\text{LB}}$ Lower bound on the problem's global optimum (€)

* $J^{\text{UB}}$ Upper bound on the problem's global optimum (€)

* $J^{\text{gap}}$ Gap between upper and lower bounds (%)

* $\epsilon^{\text{gap}}$ Tolerance parameter of termination criterion (%)

* $\epsilon^{\text{UB}}$ Tolerance parameter of upper bound update (%)

* $\rho$ Penalty parameter (-)

### Variables

* $b_{n,\omega,t}^{\text{br}}$ Binary variable associated with the barrage safety constraint of asset $n$ in scenario $\omega$ at time $t$ (-)

* $b_{n,\omega,t,i}^{\text{oc}}$ Binary variable of the i-th operational curve segment of asset $n$ in scenario $\omega$ at time $t$ (-)

* $e_{\omega,t}$ Energy offer in scenario $\omega$ at time $t$ (MWh)

* $h_{n,\omega,t}$ Hydraulic head of asset $n$ in scenario $\omega$ at time $t$ (m)

* $\bar{P}_{n,\omega,t}^{\text{H}}$ Output of asset $n$ in scenario $\omega$ at time $t$ (MW)

* $q_{n,\omega,t}^{\text{br}}$ Outflow of barrage $n$ in scenario $\omega$ at time $t$ (m$^3$/s)

* $q_{n,\omega,t}^{\text{in}}$ Inflow of asset $n$ in scenario $\omega$ at time $t$ (m$^3$/s)

* $q_{n,\omega,t}^{\text{out}}$ Outflow of asset $n$ in scenario $\omega$ at time $t$ (m$^3$/s)

* $q_{n,\omega,t}^{\text{out}}$ Outflow of turbines $n$ in scenario $\omega$ at time $t$ (m$^3$/s)

* $z_{n,\omega,t}^{\text{fl}}$ Forebay level of asset $n$ in scenario $\omega$ at time $t$ (m)
$\delta_{\omega,t}^{E\uparrow}$ Positive imbalance in scenario $\omega$ at time $t$ (MWh)

$\delta_{\omega,t}^{E\downarrow}$ Negative imbalance in scenario $\omega$ at time $t$ (MWh)

$\lambda_{n,\omega}$ Vector of dual variables corresponding to the consen-
sus constraints of asset $n$ in scenario $\omega$

# I. INTRODUCTION

## A. Literature Review

As the energy sector shifts from conventional generators to variable renewable energy sources (vRES) such as wind and solar, optimizing the bidding strategy and operation of these stochastic resources by means of controllable units has gained significant attention [1]. Hydropower plays a critical role in this transition, providing nearly half of the world's low-carbon electricity and offering an unparalleled 1500 TWh of storage capacity, dwarfing global battery storage by a factor of 2200 in 2021 [2]. Among hydropower systems, cascaded hydropower plants (CHPP) are effective in enhancing vRES market value [3] and supporting grid operations [4]. However, co-optimizing the market-based operation of CHPP and vRES is challenging due to uncertainty sources, complex spatiotemporal constraints, and complicated nonlinear dynamics.

Approaches for managing uncertainty in CHPP-vRES systems can be broadly categorized into robust optimization [5], interval optimization [6], chance-constrained programming [7], and scenario-based stochastic programming [8]. A comprehensive analysis of these methods within the broader context of power system optimization under uncertainty is detailed in [9]. Among these, scenario-based stochastic programming is often preferred for its flexibility in modeling uncertainty by incorporating multiple potential outcomes, or *scenarios*, each weighted by its associated probability. Recent studies have shown the effectiveness of this approach in managing multiple sources of uncertainty inherent in the operation of CHPP, with applications including the optimization of market participation strategies [10], maintenance scheduling [11], and day-ahead operational planning [12]. Scenario-based stochastic programming is particularly well-suited for handling unbounded uncertainty sources such as prices [13], inflow [14], and demand [15], which cannot typically be confined a priori to a specific set or bounded range, as required in robust and interval optimization approaches, respectively. Moreover, the scenario-based stochastic programming modeling approach yields solutions tailored to expected outcomes [16], in contrast to approaches like interval optimization, which typically focus on hedging against worst-case realizations [6]. However, it can be computationally inefficient, particularly with large sets of scenarios. In bidding problems, scenario-based approaches allow for the use of price-quantity bidding curves, enabling energy bid adjustments based on predicted prices and generation levels [17]. This is particularly beneficial for vRES paired with flexible resources such as hydropower units [18].

Besides uncertainty management, a key challenge in CHPP operation is balancing modeling accuracy with tractability. Comprehensive physical or operational models result in complex mixed-integer nonlinear programming (MINLP) problems [19], for which there is still no efficient off-the-shelf solver. Significant complexity arises from the coupling among units, and the complicated nonlinear dynamics involved [20].

To address this, the original model is often simplified using three main methods [21], namely aggregation, optimization/heuristic, and decomposition methods. Aggregation methods have been extensively applied to both short-term [22] and long-term [23] scheduling of CHPP, but are best suited to systems with similar reservoir and inflow characteristics, as the individual plant constraints are overlooked. Researchers have also investigated ways to simplify the MINLP model while accounting for the individual dynamics of each plant. Solving these optimization models directly with commercial solvers is often very challenging. A common method extensively used in the current literature is outer approximation [24], which decomposes the MINLP problem into a mixed-integer linear programming (MILP) master problem and a set of nonlinear programming (NLP) subproblems, although this often results in a suboptimal performance. Alternatively, the MINLP problem can be simplified in various manners, reducing to linear programming (LP) [25], MILP [26], quadratic programming [27], mixed-integer quadratic programming (MIQP) [28] or NLP [29] problems. MILP approximations are increasingly popular, offering a balance between modeling accuracy and tractability [30]. To solve these models, methods such as linear decision rules [31] and dynamic programming [29] have been applied. While rule-based approaches may oversimplify the problem, dynamic programming can be computationally expensive due to dimensionality issues.

To overcome these limitations, researchers have explored decomposition techniques. Two main strategies identified in the literature are scenario decomposition and spatial decomposition [32]. Scenario decomposition methods, such as Benders decomposition [24] or the consensus alternating direction method of multipliers (ADMM) [28], split the stochastic problem into multiple deterministic subproblems, typically one per scenario. Spatial decomposition, often achieved using Lagrangian relaxation techniques like dual decomposition [33], ADMM [34], or the auxiliary problem principle (APP) [35], breaks the problem down by units or groups of units. Both strategies can reduce computational burden compared to the centralized (non-decomposed) counterpart, however they often require assumptions of convexity to guarantee converging to the global solution. Even for MILP approximations, applying these methods to CHPP scheduling results in heuristic approaches [21].

## B. Motivation and Contributions

Co-optimized bidding strategy for CHPP and vRES is a complex stochastic MINLP problem. While MILP approximations have become popular for their ability to balance modeling accuracy versus computational tractability, they still incur significant computational burdens. Decomposition approaches have been effective in alleviating such a burden, however existing works focus on either spatial or scenario decomposition, addressing only one dimension of the problem's complexity, and most importantly, do not offer any solution guarantee in terms of optimality. To date, no decomposition method for CHPP-vRES bidding effectively scales both spatial and scenario dimensions while providing a performance guarantee.
<table><caption>TABLE I
COMPARISON OF THIS STUDY WITH THE RELEVANT RELATED LITERATURE</caption><thead><tr><th>Ref.</th><th>Computational strategy</th><th>Uncertainty modeled?</th><th>Problem</th><th>Solution method</th><th>Decomposition</th><th>Trustworthy</th></tr></thead><tbody><tr><td>[22]</td><td>Aggregation</td><td>×</td><td>MILP</td><td>Branch and bound</td><td>×</td><td>✓</td></tr><tr><td>[21]</td><td>Aggregation</td><td>✓(stochastic programming)</td><td>NLP</td><td>NLP solver</td><td>×</td><td>✓</td></tr><tr><td>[29]</td><td>Optimization</td><td>×</td><td>NLP</td><td>Dynamic programming</td><td>×</td><td>✓</td></tr><tr><td>[5]</td><td>Optimization</td><td>✓(robust optimization)</td><td>LP</td><td>LP solver</td><td>×</td><td>✓</td></tr><tr><td>[26]</td><td>Optimization</td><td>✓(interval programming)</td><td>MILP</td><td>Branch and bound</td><td>×</td><td>✓</td></tr><tr><td>[14]</td><td>Optimization</td><td>✓(stochastic programming)</td><td>MILP</td><td>Branch and bound</td><td>×</td><td>✓</td></tr><tr><td>[15]</td><td>Optimization</td><td>✓(stochastic programming)</td><td>MILP</td><td>Branch and bound</td><td>×</td><td>✓</td></tr><tr><td>[12]</td><td>Optimization</td><td>✓(stochastic programming)</td><td>MILP</td><td>Branch and bound</td><td>×</td><td>✓</td></tr><tr><td>[10]</td><td>Optimization</td><td>✓(stochastic programming)</td><td>MILP</td><td>Dynamic programming</td><td>×</td><td>✓</td></tr><tr><td>[31]</td><td>Optimization</td><td>✓(stochastic programming)</td><td>LP</td><td>Linear decision rules</td><td>×</td><td>✓</td></tr><tr><td>[35]</td><td>Decomposition</td><td>×</td><td>MILP</td><td>APP</td><td>Spatial</td><td>×</td></tr><tr><td>[34]</td><td>Decomposition</td><td>×</td><td>MIQP</td><td>ADMM</td><td>Spatial</td><td>×</td></tr><tr><td>[33]</td><td>Decomposition</td><td>✓(stochastic programming)</td><td>MILP</td><td>Dual decomposition</td><td>Spatial</td><td>×</td></tr><tr><td>[11]</td><td>Decomposition</td><td>✓(stochastic programming)</td><td>MILP</td><td>Benders decomposition</td><td>Scenario</td><td>×</td></tr><tr><td>[24]</td><td>Decomposition</td><td>✓(stochastic programming)</td><td>MINLP</td><td>Benders decomposition</td><td>Scenario</td><td>×</td></tr><tr><td>[28]</td><td>Decomposition</td><td>✓(stochastic programming)</td><td>MIQP</td><td>Consensus ADMM</td><td>Scenario</td><td>×</td></tr><tr><td>This paper</td><td>Decomposition</td><td>✓(stochastic programming)</td><td>MILP</td><td>Consensus ADMM</td><td>Spatial & scenario</td><td>✓</td></tr></tbody></table>

This paper develops a scalable co-optimization method for the joint bidding of cascaded hydropower, wind and solar energy units in the day-ahead market, accounting for uncertainties in water inflows, prices, and renewable power generation. A MILP formulation is provided, which is then decomposed using a consensus ADMM approach. Since a straightforward application of ADMM to a mixed-integer problem inherently results in a heuristic method, the framework is extended to incorporate trustworthy upper and lower bounds on the distributed solution. This enhancement not only accelerates the algorithm but also facilitates the assessment of solution quality, ensuring that the solution is $\epsilon$-suboptimal, meaning its objective value is guaranteed to be at most $\epsilon$ (e.g., 0.01%) worse than the global optimal solution. In other words, the proposed distributed algorithm offers a performance guarantee, as it provides solutions that operators can confidently rely on. Finally, drawing on recent advancements in consensus ADMM [36], the existence and validity of the bounds is demonstrated. Table I compares this study with the relevant related literature.

The key contributions of this paper are summarized as follows. First, a scenario-spatial decomposition of the centralized stochastic CHPP-vRES bidding strategy is proposed. This involves decomposing the problem with respect to three sets of coupling constraints: consistency constraints across price, inflow, and vRES scenarios (scenario coupling); power balancing among the aggregated energy units (spatial coupling); and hydraulic links within the hydropower cascade (spatial coupling). This full-scale decomposition distinguishes the proposed approach from previous studies that either address the centralized problem directly, e.g., [4], [5], [7], [31], or decouple only one set of constraints, e.g., [24], [28], [33]–[35]. Second, a consensus ADMM-based distributed algorithm is developed to derive trustworthy upper and lower bounds on the global optimum of the bidding problem. As a distinct feature, this method offers a performance guarantee, even for MILP problems, by using bounds to quantify the “distance” between the distributed solution and the global optimum. In contrast, existing decomposition approaches typically yield heuristic solutions when applied to mixed-integer models, e.g., [24], [28], [33]–[35]. Finally, the paper offers a mathematical

proof demonstrating the existence and validity of the bounds.

The proposed approach is validated through a case study simulating the real-world system managed by the French aggregator Compagnie Nationale du Rhône [37]. The model integrates real operational curves to account for the individual safety constraints of each hydropower plant. The non-convex nature of these curves introduces significant complexity into the optimization problem. Furthermore, dynamic price-quantity bidding curves are derived to enhance the adaptability of the bidding strategy to price forecasts, providing greater flexibility compared to fixed hourly bids, which are characterized by a uniform quantity offered for all price levels. Simulation results assess the modeling accuracy, the aggregator's performance in the market, and the scalability of the distributed method. An ex-post out-of-sample validation is performed to determine the minimum number of scenarios necessary to adequately represent the underlying uncertainties.

The remainder of the article is organized as follows: Section II details the problem and its decomposition, Section III discusses the results, and Section IV concludes the study.

## II. METHODOLOGY

This section formulates the problem and outlines the proposed distributed method with trustworthy bounds. Let $n \in \mathbb{N}$ denote a hydropower plant within the cascaded system. The problem is formulated in discrete time over a finite set of sampling times $t \in T$, with a sampling period of $\Delta_T$. The scenario $\omega \in \Omega$ denotes a potential realization of uncertainty. The uncertain parameters are denoted by the hat symbol ($\hat{\cdot}$). The symbol $:==$ denotes the definition of variables, vectors, functions, and sets, while $|\cdot|$ denotes the cardinality of a set.

### A. Problem Statement

The goal is to optimize the participation of a price-taking energy producer in the day-ahead market, known as the bidding problem [40]. The system mimics a real-world virtual power plant (VPP) operated by Compagnie Nationale du Rhône, comprising 12 CHPP (totaling 2.2 GW of installed capacity), 6 wind farms (278.4 MW), and 4 solar farms (60.8 MW), as
<figure><img src="image_3.png" /><figcaption>Fig. 1. Diagram of the VPP under consideration, comprising CHPP and vRES.</figcaption></figure>

shown in Fig. 1. The study considers three sources of uncertainty—water inflow, vRES output, and market prices—which are modeled using a finite set of scenarios.

The energy volumes offered to the market are denoted by $e_{\omega,t}$, and the day-ahead energy price forecast by $\hat{\pi}^E$. A dual pricing settlement scheme is employed, with separate prices for positive imbalances, $\delta^{E\uparrow}$ (energy injection falling short of the offer), and negative imbalances, $\delta^{E\downarrow}$ (energy injection exceeding the offer). The forecasted penalty prices for these imbalances are denoted by $\hat{\pi}^{E\uparrow}$ and $\hat{\pi}^{E\downarrow}$, respectively. Bids are submitted to the market operator as price-quantity pairs, forming a set of price-quantity bidding curves [17]. The optimization problem considers all scenario prices as potential bidding points, with the energy produced for each hour and price determining the bid. The use of price scenarios as potential bidding points has been extensively documented in the literature on optimal bidding strategies for VPPs [38], particularly within the framework of dual pricing settlement schemes [39]. For a comprehensive comparison of the price-quantity bidding strategy with alternative price-only and quantity-only strategies, the reader is referred to [18]. Consistent with common rules in many European markets [40], all supply bids within an hour must be non-decreasing as the price increases.

## B. Centralized Optimization Problem

The VPP in Fig. 1 is modeled with three key assumptions: (i) vRES are uncontrollable, facilitating the assessment of the cascade's ability to manage their variability; (ii) tailrace water levels remain constant due to the minimal short-term impact of run-of-the-river plant discharges; and (iii) all turbines within each hydropower plant are of the same type. In accordance

<figure><img src="image_4.png" /><figcaption>Fig. 2. Example of an operational curve from Compagnie Nationale du Rhône.</figcaption></figure>

with assumption (iii), the formulation presumes that the turbines are homogeneous within each individual hydropower plant, while allowing for heterogeneous turbines across different plants within the cascaded system. This assumption is adopted to streamline the subsequent discussion and notation. Nevertheless, the proposed formulation can be readily extended to accommodate heterogeneous turbines within a single plant, for example by following the modeling approach presented in [41].

The dynamics of the forebay reservoir water levels, denoted by $z^{\text{fbl}}$, are computed as a function of the reservoir's surface area $S$, inflow $q^{\text{in}}$, outflow $q^{\text{out}}$, and the initial forebay water level $Z^{\text{fbl},0}$. For each asset pair, a specific point between the two reservoirs, termed the control point, is selected, and the forebay water level is derived from the volume balance at this control point, leading to the following state-space model:

$$z_{n,\omega,t}^{\text{fbl}} = z_{n,\omega,t-1}^{\text{fbl}} + \frac{(q_{n,\omega,t}^{\text{in}} - q_{n,\omega,t}^{\text{out}}) \Delta_T}{S_{n,\omega,t}}, \quad \forall n, \forall \omega, \forall t, \quad (1a)$$

$$z_{n,\omega,0}^{\text{fbl}} = z_{n,\omega,0|T|}^{\text{fbl},0} = Z_{n}^{\text{fbl},0}, \quad \forall n, \forall \omega. \quad (1b)$$

The inflow into each reservoir is composed of the water discharged from the upstream hydropower system through the barrage and turbines, denoted by $q^{\text{br}}$ and $q^{\text{tr}}$, respectively, and the external inflow, denoted by $\hat{Q}^{\text{ext}}$. Constraints (2) describe the inflow into the n-th reservoir under scenario $\omega$ at time $t$:

$$q_{n,\omega,t}^{\text{in}} = q_{n-1,\omega,t-\tau_{n-1,n}^{\text{br}}}^{\text{br}} + q_{n-1,\omega,t-\tau_{n-1,n}^{\text{tr}}}^{\text{tr}} + \hat{Q}_{n,\omega,t}^{\text{ext}}, \quad \forall n, \forall \omega, \forall t, \quad (2)$$

where $\tau_{n-1,n}^{\text{tr}}$ and $\tau_{n-1,n}^{\text{br}}$ represent the water propagation time from the plant $n-1$ to the control point of plant $n$ via turbines and barrage, respectively. Similarly, the outflow is given by

$$q_{n,\omega,t}^{\text{out}} = q_{n,\omega,t}^{\text{br}} + q_{n,\omega,t}^{\text{tr}}, \quad \forall n, \forall \omega, \forall t. \quad (3)$$

To ensure stability and safety, CHPP operate under strict operational curves. Fig. 2 shows an example of an operational curve for a hydropower plant. These curves are developed based on approved guidelines from relevant authorities and take into account various constraints, including navigation, irrigation, nuclear safety, agriculture, and tourism. During low-flow periods, water levels are kept high for navigation, while during floods, levels are controlled to mimic pre-barrage conditions. Reservoirs provide storage capacity only during normal inflow conditions, i.e., during energetic periods.
The operational curves are approximated using the piece-wise linear constraints in (4). This involves dividing the inflow range into segments $Q^{in}$, each linked to fixed forebay water levels $\underline{Z}^{fbl}$ and $\bar{Z}^{fbl}$. A set of inflow segments $S_n^{oc}$ is computed for each hydropower plant $n$:

$$q_{n,\omega,t}^{in} \geq b_{n,\omega,t,i}^{oc} Q_{n,i}^{in}, \quad \forall i \in S_n^{oc}, \forall n, \forall \omega, \forall t, \quad (4a)$$

$$\begin{array}{l} q_{n,\omega,t}^{in} < b_{n,\omega,t,i}^{oc} Q_{n,i+1}^{in} + (1 - b_{n,\omega,t,i}^{oc}) M^{oc}, \\ \qquad \forall i \in S_n^{oc}, \forall n, \forall \omega, \forall t. \end{array} \quad (4b)$$

Here, $M^{oc}$ is a large positive constant, and the binary variables $b^{oc}$ are constrained by

$$\sum_{i \in S_n^{oc}} b_{n,\omega,t,i}^{oc} = 1, \quad \forall n, \forall \omega, \forall t. \quad (5)$$

Constraints (6) enforce the forebay water level limits:

$$\sum_{i \in S_n^{oc}} \underline{Z}_{n,i}^{fbl} b_{n,\omega,t,i}^{oc} \leq z_{n,\omega,t}^{fbl} \leq \sum_{i \in S_n^{oc}} \bar{Z}_{n,i}^{fbl} b_{n,\omega,t,i}^{oc}, \\ \forall n, \forall \omega, \forall t. \quad (6)$$

The turbines are constrained by the ramp limit:

$$q_{n,\omega,t}^{tr} - q_{n,\omega,t-1}^{tr} \leq \Delta q_{n}^{tr}, \quad \forall n, \forall \omega, \forall t, \quad (7a)$$

$$q_{n,\omega,t-1}^{tr} - q_{n,\omega,t}^{tr} \leq \Delta q_{n}^{tr}, \quad \forall n, \forall \omega, \forall t. \quad (7b)$$

For safety reasons, the opening of the barrage is allowed exclusively when the reservoir is at full capacity:

$$b_{n,\omega,t}^{br} \underline{Q}_{n}^{br} \leq q_{n,\omega,t}^{br} \leq b_{n,\omega,t}^{br} M^{br}, \quad \forall n, \forall \omega, \forall t. \quad (8)$$

Here, $\underline{Q}^{br}$ is the minimum barrage discharge, $M^{br}$ is a large positive constant, and the binary variables $b^{br}$ satisfy

$$b_{n,\omega,t}^{br} \leq 1 - \frac{\sum_{i \in S_n^{oc}} \bar{Z}_{n,i}^{fbl} b_{n,\omega,t,i}^{oc} - z_{n,\omega,t}^{fbl}}{\max\left(\left\{\bar{Z}_{n,i}\right\}_{i \in S_n^{oc}}\right)}, \quad \forall n, \forall \omega, \forall t. \quad (9)$$

Similarly, when the turbines operate, their discharge is constrained by

$$\underline{Q}_{n}^{tr} \leq q_{n,\omega,t}^{tr} \leq \bar{Q}_{n}^{tr}, \quad \forall n, \forall \omega, \forall t. \quad (10)$$

The hydropower output, $p^H$, can be calculated as a function of the turbines outflow, net hydraulic head $h$, and turbine efficiency $\eta^{tr}$, as described by

$$p_{n,\omega,t}^{H} = w g \eta_{n}^{tr} q_{n,\omega,t}^{tr} h_{n,\omega,t} 10^{-6}, \quad \forall n, \forall \omega, \forall t, \quad (11)$$

where $g$ is the acceleration due to gravity and $w$ is the density of water. The factor $10^{-6}$ is used to convert the power output from watts to megawatts. The head is defined by

$$h_{n,\omega,t} = z_{n,\omega,t}^{fbl} - Z_{n}^{tr}, \quad \forall n, \forall \omega, \forall t. \quad (12)$$

Here, $Z^{tr}$ denotes the constant tailrace reservoir water level. The hydropower generation function exhibits nonlinearity due to its bilinear dependence on $q^{tr}$ and $h$. To address this, the

<figure><img src="image_5.png" /><figcaption>Fig. 3. Visualization of the McCormick approximation in (13), assuming a turbine efficiency of 95%.</figcaption></figure>

McCormick approximation [42] in (13) is employed to replace the bilinear term in (11) with a convex envelope:

$$p_{n,\omega,t}^{H} \geq C_n \left( \underline{Q}_{n}^{tr} h_{n,\omega,t} + \underline{H}_{n} q_{n,\omega,t}^{tr} - \underline{Q}_{n}^{tr} \underline{H}_{n} \right), \quad (13a)$$

$$p_{n,\omega,t}^{H} \geq C_n \left( \bar{Q}_{n}^{tr} h_{n,\omega,t} + \bar{H}_{n} q_{n,\omega,t}^{tr} - \bar{Q}_{n}^{tr} \bar{H}_{n} \right), \quad (13b)$$

$$p_{n,\omega,t}^{H} \leq C_n \left( \underline{Q}_{n}^{tr} h_{n,\omega,t} + \underline{H}_{n} q_{n,\omega,t}^{tr} - \underline{Q}_{n}^{tr} \underline{H}_{n} \right), \quad (13c)$$

$$p_{n,\omega,t}^{H} \leq C_n \left( \bar{Q}_{n}^{tr} h_{n,\omega,t} + \bar{H}_{n} q_{n,\omega,t}^{tr} - \bar{Q}_{n}^{tr} \bar{H}_{n} \right). \quad (13d)$$

Here, $C_n = 10^{-6} w g \eta_n^{tr}$, while $\underline{H}$ and $\bar{H}$ denote the maximum and minimum hydraulic head, respectively. Fig. 3 illustrates the McCormick approximation in (13). This linearization is a standard method for constructing convex envelopes around bilinear functions [43], with broad applications in hydropower modeling [15]. It is considered one of the most accurate and tight approximations available for bilinear terms [44]. A detailed error analysis is provided in [45], showing that the largest deviations occur at the intersections of the bounding hyperplanes defined by (13). For instance, assuming a turbines discharge in the range [0, 1000] m³/s, maximum head of 5 m, 95% turbine efficiency, water density of 1000 kg/m³, and gravitational acceleration of 9.81 m/s², this approximation yields a maximum error of ~3% of the generation capacity.

The range of admissible values for $p^H$ is enforced by

$$0 \leq p_{n,\omega,t}^{H} \leq \bar{P}_{n,\omega,t}, \quad \forall n, \forall \omega, \forall t. \quad (14)$$

Constraints (15) define the VPP's power balance as

$$\left( \sum_{n \in \mathcal{N}} p_{n,\omega,t}^{H} + \hat{P}_{\omega,t}^{W} + \hat{P}_{\omega,t}^{S} \right) \Delta_{T} - e_{\omega,t} - \delta_{\omega,t}^{E\downarrow} + \delta_{\omega,t}^{E\uparrow} = 0, \quad \forall \omega, \forall t, \quad (15)$$

where $\hat{P}^W$ and $\hat{P}^S$ denote the total wind and solar power generation of the VPP, respectively.

The aggregator minimizes the objective function defined as

$$J := \sum_{\omega \in \Omega} \sum_{t \in T} P_{\omega} \left( \hat{\pi}_{\omega,t}^{E\uparrow} \delta_{\omega,t}^{E\uparrow} - \hat{\pi}_{\omega,t}^{E\downarrow} \delta_{\omega,t}^{E\downarrow} - \hat{\pi}_{\omega,t}^{E} e_{\omega,t} \right), \quad (16)$$
representing the cost incurred by the VPP in the day-ahead market. Minimizing this cost is equivalent to maximizing the VPP's profit. Here, $P_\omega$ is the probability of scenario $\omega$. Constraints (17) ensure price-quantity bidding curves taking a non-decreasing form:

$$ (\hat{\pi}_{\omega,t}^{E} - \hat{\pi}_{\omega',t}^{E})(e_{\omega,t} - e_{\omega',t}) \geq 0, \forall \omega' \in \Omega \setminus \{\omega\}, \\ \forall \omega, \forall t. \quad (17) $$

Finally, the set of decision variables $z$ is defined as

$$ z := \left\{ e_{\omega,t}, \delta_{\omega,t}^{\mathrm{E}\uparrow}, \delta_{\omega,t}^{\mathrm{E}\downarrow}, q_{n,\omega,t}^{\mathrm{tr}}, q_{n,\omega,t}^{\mathrm{br}}, z_{n,\omega,t}^{\mathrm{fbl}}, p_{n,\omega,t}^{\mathrm{H}}, b_{n,\omega,t}^{\mathrm{br}}, b_{n,\omega,t}^{\mathrm{oc}} \right\}_{i \in S_n^{\mathrm{oc}}, n \in N, \omega \in \Omega, t \in T}. \quad (18) $$

Then, the MILP given in (19) defines the bidding problem:

$$ \min_{z} J \quad (19a) $$

$$ \text{subject to:} \quad (1) - (10), (12), (14), (15), (17), \quad (19b) $$

$$ (13), \forall n, \forall \omega, \forall t. \quad (19c) $$

### C. From Centralized to Distributed Optimization

The centralized optimization problem (19) is characterized by spatial coupling, due to the hydraulic links among CHP in (1) and the electrical coupling in (15), as well as the scenario coupling in (17). All coupling constraints involve the global variables $q_{n,\omega}^{\mathrm{tr}} := [q_{n,\omega,1}^{\mathrm{tr}}, \dots, q_{n,\omega,|T|}^{\mathrm{tr}}]^{\mathrm{T}}, q_{n,\omega}^{\mathrm{br}} := [q_{n,\omega,1}^{\mathrm{br}}, \dots, q_{n,\omega,|T|}^{\mathrm{br}}]^{\mathrm{T}}$ and $p_{n,\omega}^{\mathrm{H}} := [p_{n,\omega,1}^{\mathrm{H}}, \dots, p_{n,\omega,|T|}^{\mathrm{H}}]^{\mathrm{T}}$. By duplicating these variables, the centralized problem (19) is reformulated as a consensus problem, which reads as

$$ \min_{z^{\mathrm{B}}, \{z_{n,\omega}^{\mathrm{H}}\}_{n \in N, \omega \in \Omega}} J \quad (20a) $$

$$ \text{subject to} \quad z^{\mathrm{B}} \in \Xi(\beta), \quad (20b) $$

$$ z_{n,\omega}^{\mathrm{H}} \in \Gamma_{n,\omega}(\theta_{n,\omega}), \forall n, \forall \omega, \quad (20c) $$

$$ z_{n,\omega}^{\mathrm{B,H}} = \tilde{z}_{n,\omega}^{\mathrm{B,H}} : \lambda_{n,\omega}, \forall n, \forall \omega. \quad (20d) $$

Here, $z^{\mathrm{B}}$ collects the decision variables associated with the power balancing of the VPP, as defined by

$$ z^{\mathrm{B}} := \left\{ e_{\omega,t}, \delta_{\omega,t}^{\mathrm{E}\uparrow}, \delta_{\omega,t}^{\mathrm{E}\downarrow}, \bar{p}_{n,\omega}^{\mathrm{H}} \right\}_{n \in N, \omega \in \Omega, t \in T}. \quad (21) $$

The decision variables associated with the dynamics of each hydropower plant are collected in $z_{n,\omega}^{\mathrm{H}}$, as defined by

$$ z_{n,\omega}^{\mathrm{H}} := \left\{ b_{n,\omega,t}^{\mathrm{br}}, b_{n,\omega,t}^{\mathrm{oc}}, h_{n,\omega,t}, z_{n,\omega,t}^{\mathrm{fbl}}, \tilde{q}_{n,\omega}^{\mathrm{tr}}, \tilde{q}_{n,\omega}^{\mathrm{tr}}, \tilde{q}_{n,\omega}^{\mathrm{br}}, \tilde{q}_{n,\omega}^{\mathrm{br}}, \tilde{p}_{n,\omega}^{\mathrm{H}} \right\}_{t \in T}. \quad (22) $$

The feasible sets of $z^{\mathrm{B}}$ and $z_{n,\omega}^{\mathrm{H}}$ are denoted by $\Xi(\beta)$ and $\Gamma_{n,\omega}(\theta_{n,\omega})$, respectively. The former is defined by (14), (15) and (17). The latter is defined by (1)-(10) and (12)-(14). Vectors $\beta$ and $\theta_{n,\omega}$ collect the parameters associated with the constraints defining the feasible sets $\Xi$ and $\Gamma_{n,\omega}$, respectively.

Moreover, the vector $z_{n,\omega}^{\mathrm{B,H}} \in \mathbb{R}^{6 \times |T|}$ collects the global variables of the problem, as defined by

$$ z_{n,\omega}^{\mathrm{B,H}} := [q_{n-1,\omega}^{\mathrm{tr}}, q_{n-1,\omega}^{\mathrm{br}}, q_{n,\omega}^{\mathrm{tr}}, q_{n,\omega}^{\mathrm{br}}, p_{n,\omega}^{\mathrm{H}}, p_{n,\omega}^{\mathrm{H}}]^{\mathrm{T}}. \quad (23) $$

<figure><img src="image_6.png" /><figcaption>Fig. 4. Example of decomposed structure of problem (19) when reformulated as the consensus problem (20). The edges in the bipartite graph represent the consistency constraint linking global variables with their local copies.</figcaption></figure>

The vector $\tilde{z}_{n,\omega}^{\mathrm{B,H}} \in \mathbb{R}^{6 \times |T|}$ collects the local copies of the global variables in $z_{n,\omega}^{\mathrm{B,H}}$, as defined by

$$ \tilde{z}_{n,\omega}^{\mathrm{B,H}} := [\tilde{q}_{n-1,\omega}^{\mathrm{tr}}, \tilde{q}_{n-1,\omega}^{\mathrm{br}}, \tilde{q}_{n,\omega}^{\mathrm{tr}}, \tilde{q}_{n,\omega}^{\mathrm{br}}, \tilde{p}_{n,\omega}^{\mathrm{H}}, \tilde{p}_{n,\omega}^{\mathrm{H}}]^{\mathrm{T}}. \quad (24) $$

Consistency between global variables and local copies is enforced through the consensus constraints in (20d). The corresponding vector of dual variables $\lambda_{n,\omega} \in \mathbb{R}^{6 \times |T|}$ is defined as

$$ \lambda_{n,\omega} := [\tilde{\lambda}_{n-1,\omega}^{\mathrm{tr}}, \tilde{\lambda}_{n-1,\omega}^{\mathrm{br}}, \tilde{\lambda}_{n,\omega}^{\mathrm{tr}}, \tilde{\lambda}_{n,\omega}^{\mathrm{br}}, \tilde{\lambda}_{n,\omega}^{\mathrm{p}}, \tilde{\lambda}_{n,\omega}^{\mathrm{p}}]^{\mathrm{T}}. \quad (25) $$

Clearly, problem (20) exhibits a decomposable structure. Subproblem (20b) pertains to the VPP power balancing and is referred to as the power balancing subproblem. Subproblems (20c) address the hydropower plants individual dynamics and are referred to as the local hydropower scheduling subproblems. An example of the resulting decomposed structure is shown in Fig. 4 for $|N| = 3$ and $|\Omega| = 1$. The computational advantage of this decomposition is that the original stochastic MILP problem, solved for all $n$ and $\omega$, is now broken down into $|N| \times |\Omega|$ deterministic MILP subproblems, and one stochastic LP subproblem, solved for all $n$ and $\omega$.

As illustrated in Fig. 4, the original global variables of the centralized problem (19), namely $q_{n,\omega}^{\mathrm{tr}}$, $q_{n,\omega}^{\mathrm{br}}$, and $p_{n,\omega}^{\mathrm{H}}$, are duplicated in the proposed decomposition. For each set of decoupled coupling constraints, two distinct sets of subproblems are introduced, each operating on its own set of duplicated variables, which are distinguished by distinct superscripts. For example, the global variable $p_{n,\omega}^{\mathrm{H}}$ is decomposed into local copies $\tilde{p}_{n,\omega}^{\mathrm{H}}$ for the power balancing subproblem, and distinct local copies $\tilde{p}_{n,\omega}^{\mathrm{H}}$ for the hydropower scheduling subproblems. The dual variables associated with the consensus constraints (20d), which are enforced for each pair of global and local variable copies, are denoted using the same superscript as their corresponding local variables. The resulting scenario-spatial decomposition enables a fully parallelizable solution of the original centralized problem (19) via the distributed optimization algorithm detailed in the following subsection.

Analyzing the participation of the VPP examined in this study in markets other than the day-ahead market is outside
the scope of this work. For a comprehensive evaluation of its profitability when simultaneously participating in both day-ahead and ancillary services markets, the reader is referred to [47]. Nonetheless, it is worth noting that the proposed scenario-spatial decomposition not only improves the scalability of the bidding strategy compared to traditional centralized approaches, but also enhances its flexibility and modularity. Specifically, the proposed decomposition approach decouples the local hydropower scheduling from the global power balancing of the VPP, resulting in two distinct sets of subproblems (see Fig. 4). In multi-market extensions, modifications to the bidding optimization model are typically confined to the global balancing subproblem and may include adjustments to the power balance constraints, updates to the objective function to reflect additional market revenues or penalties, and the inclusion of new constraints on the VPP's cumulative power generation [47]. The local hydropower scheduling subproblems in Fig. 4 remain unchanged. This modular structure supports both scalability and adaptability to the dynamic evolution of the VPP operations and objectives.

### D. Distributed Optimization with Trustworthy Bounds

The reformulation presented in (20) enables the application of various distributed optimization techniques. Among others, the consensus ADMM has demonstrated significant effectiveness [46]. Let $k$ denote the current iteration of the distributed algorithm, and $\rho^k$ denote the penalty parameter, which is updated dynamically, as detailed in [46]. To simplify the following formulation, define the matrix $\mathbf{A} \in \mathbb{R}^{6|T| \times 6|T|}$ as a diagonal matrix with a $5|T| \times 5|T|$ identity matrix in the top-left block and the last $|T|$ diagonal elements set to zero. By applying the consensus ADMM algorithm to (20), the iterative steps outlined below are derived.

**Step I.** Local primal variables update (26)-(27):

$$ z_{n,\omega}^{\mathbf{H}^{k+1}} = \underset{z_{n,\omega}^{\mathbf{H}} \in \Gamma_{n,\omega}}{\operatorname{argmin}} \left\{ \left( \mathbf{A} \boldsymbol{\lambda}_{n,\omega}^{k} \right)^{\top} \bar{\mathbf{z}}_{n,\omega}^{\mathbf{B},\mathbf{H}} + \frac{\rho}{2} \left\| \mathbf{A} \left( \bar{\mathbf{z}}_{n,\omega}^{\mathbf{B},\mathbf{H}} - z_{n,\omega}^{\mathbf{B},\mathbf{H}^{k}} \right) \right\|_{2}^{2} \right\}, \forall n, \forall \omega, \quad (26) $$

$$ z^{\mathbf{B}^{k+1}} = \underset{z^{\mathbf{B}} \in \Xi}{\operatorname{argmin}} \left\{ J + \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} \left( \left( \bar{\boldsymbol{\lambda}}_{n,\omega}^{\mathbf{P}} \right)^{\top} \bar{\boldsymbol{p}}_{n,\omega}^{\mathbf{H}} + \frac{\rho}{2} \left\| \bar{\boldsymbol{p}}_{n,\omega}^{\mathbf{H}} - \boldsymbol{p}_{n,\omega}^{\mathbf{H}} \right\|_{2}^{2} \right) \right\}. \quad (27) $$

**Step II.** Global primal variables update in (28):

$$ q_{n,\omega}^{\mathbf{tr}} \stackrel{k+1}{=} \frac{1}{2} \left( \bar{q}_{n,\omega}^{\mathbf{tr}} \stackrel{k+1}{=} \tilde{q}_{n,\omega}^{\mathbf{tr}} \right), \forall n, \forall \omega, \quad (28a) $$

$$ q_{n,\omega}^{\mathbf{br}} \stackrel{k+1}{=} \frac{1}{2} \left( \bar{q}_{n,\omega}^{\mathbf{br}} \stackrel{k+1}{=} \tilde{q}_{n,\omega}^{\mathbf{br}} \right), \forall n, \forall \omega, \quad (28b) $$

$$ p_{n,\omega}^{\mathbf{H}^{k+1}} = \frac{1}{2} \left( \bar{p}_{n,\omega}^{\mathbf{H}^{k+1}} + \tilde{p}_{n,\omega}^{\mathbf{H}^{k+1}} \right), \forall n, \forall \omega. \quad (28c) $$

**Step III.** Dual variables update in (29):

$$ \lambda_{n,\omega}^{k+1} = \lambda_{n,\omega}^{k} + \rho \left( \bar{z}_{n,\omega}^{\mathbf{B},\mathbf{H}^{k+1}} - z_{n,\omega}^{\mathbf{B},\mathbf{H}^{k+1}} \right), \forall n, \forall \omega. \quad (29) $$

The simulation results to be presented in the next section show that applying ADMM to solve (20) often yields high-quality solutions, highlighting its potential. However, the method's convergence depends on the convexity of the subproblems, which is not met here due to the presence of binary variables, making ADMM a promising yet heuristic approach for solving (20). To make it trustworthy, a method is introduced to derive effective upper and lower bounds on the optimal objective function value of (19), denoted by $J^*$, based on the distributed solution computed at each iteration of the consensus ADMM.

The initial lower bound, $J^{\mathrm{LB}^0}$, is obtained by solving the linear relaxation of (19). This relaxation expands the feasible space, inevitably yielding a new optimum that is less than or equal to $J^*$. Conversely, the initial upper bound, $J^{\mathrm{UB}^0}$, is established by constraining the feasible space, which can be achieved by eliminating the flexibility of CHPP. This is done by setting the upper and lower limits on the forebay water level of each hydropower plant to be identical, thus fixing the binary variables in (6). Consequently, the binary variables $b^{\mathrm{br}}$ in (8) are fixed to 1 for all plants, at every time step and scenario, allowing the constraints (9) to be neglected. Determining both $J^{\mathrm{LB}^0}$ and $J^{\mathrm{UB}^0}$ involves solving two LP problems before executing the distributed algorithm, with minimal computational impact. The proposed method then iteratively updates $J^{\mathrm{LB}^0}$ and $J^{\mathrm{UB}^0}$ as the ADMM algorithm proceeds, following two distinct strategies.

For the upper bound, each iteration $k$ yields an updated value, $J^{\mathrm{UB}^k}$, through a quasi-projection step. Let $z^k$ denote the primal variable value obtained via ADMM. Drawing inspiration from the outer approximation strategy [24], $z^k$ is projected onto the LP formed by fixing the binary variables in (19) to the values of $z^k$. This projection, denoted by $z^{\mathrm{prj}^k}$, remains feasible within (19) since it is derived from a projection onto a new feasible space formed by constraining (19)—specifically, fixing the binary variables. Thus, $J^* \le J^{\mathrm{UB}^k} := J(z^{\mathrm{prj}^k})$.

Proposition 1 describes the process of using dual variables associated with the consensus constraints (20) for the iterative updating of the lower bound $J^{\mathrm{LB}^k}$. It extends the theoretical findings from [36]—established for consensus ADMM in scenario decomposition—to the case where the problem is decomposed across both spatial and scenario dimensions. Notably, although stated for $J$ in (16), the result of Proposition 1 readily extends to any alternative convex objective function.

**Proposition 1.** Let $\lambda := \{\lambda_{n,\omega}\}_{n \in \mathcal{N}, \omega \in \Omega}$ satisfy the following conditions (component-wise):

$$ \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} P_{\omega} \left( \bar{\lambda}_{n,\omega}^{\mathbf{P}} + \tilde{\lambda}_{n,\omega}^{\mathbf{P}} \right) = 0, \quad (30) $$

$$ \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} P_{\omega} \left( \bar{\lambda}_{n,\omega}^{\mathrm{tr}} + \tilde{\lambda}_{n,\omega}^{\mathrm{tr}} \right) = 0, \quad (31) $$

$$ \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} P_{\omega} \left( \bar{\lambda}_{n,\omega}^{\mathrm{br}} + \tilde{\lambda}_{n,\omega}^{\mathrm{br}} \right) = 0. \quad (32) $$

Let

$$ D^{\mathrm{B}}(\lambda) := \min_{z^{\mathrm{B}} \in \Xi} \left\{ J + \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} P_{\omega} \left( \bar{\lambda}_{n,\omega}^{\mathbf{P}} \right)^{\top} \boldsymbol{p}_{n,\omega}^{\mathbf{H}} \right\}, \quad (33) $$
and

$$D_{n,\omega}^{\mathrm{H}}(\lambda_{n,\omega}) := \min_{z_{n,\omega}^{\mathrm{H}} \in \Gamma_{n,\omega}} \left\{ (\mathbf{A}\lambda_{n,\omega})^{\top} z_{n,\omega}^{\mathrm{B},\mathrm{H}} \right\}, \quad \forall n, \forall \omega. \quad (34)$$

Then

$$D(\lambda) := D^{\mathrm{B}}(\lambda) + \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} P_{\omega} D_{n,\omega}^{\mathrm{H}}(\lambda_{n,\omega}) \leq J^{*}. \quad (35)$$

*Proof.* Let $\{b^{\mathrm{br}^{*}}, b^{\mathrm{oc}^{*}}e^{*}, \delta^{\mathrm{E}^{*}}, p^{\mathrm{H}^{*}}, q^{\mathrm{tr}^{*}}, q^{\mathrm{br}^{*}}, z^{\mathrm{fbl}^{*}}\}$ be an optimal solution of the original problem (19). Feasibility implies that this optimal solution is also feasible within the relaxed feasible spaces of the subproblems (20b) and (20c), i.e., it is feasible in $\Xi$ and $\Gamma_{n,\omega}$, $\forall n, \forall \omega$, respectively. Thus:

$$D^{\mathrm{B}}(\lambda) \leq J(e^{*}, \delta^{\mathrm{E}^{*}}) + \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} P_{\omega} \left( \bar{\lambda}_{n,\omega}^{\mathrm{P}} \right)^{\top} p_{n,\omega}^{\mathrm{H}^{*}}, \quad (36)$$

and

$$D^{\mathrm{H}}(\lambda_{n,\omega}) \leq (\bar{\lambda}_{n,\omega}^{\mathrm{P}})^{\top} p_{n,\omega}^{\mathrm{H}^{*}} + (\bar{\lambda}_{n-1,\omega}^{\mathrm{tr}})^{\top} q_{n-1,\omega}^{\mathrm{tr}^{*}} \\ + (\bar{\lambda}_{n-1,\omega}^{\mathrm{br}})^{\top} q_{n-1,\omega}^{\mathrm{br}^{*}} + (\bar{\lambda}_{n,\omega}^{\mathrm{tr}})^{\top} q_{n,\omega}^{\mathrm{tr}^{*}} \\ + (\bar{\lambda}_{n,\omega}^{\mathrm{br}})^{\top} q_{n,\omega}^{\mathrm{br}^{*}}, \quad \forall n, \forall \omega. \quad (37)$$

Then

$$D(\lambda) \leq \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} P_{\omega} \left( \bar{\lambda}_{n,\omega}^{\mathrm{P}} + \tilde{\lambda}_{n,\omega}^{\mathrm{P}} \right)^{\top} p_{n,\omega}^{\mathrm{H}^{*}} \\ + \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} P_{\omega} \left( \bar{\lambda}_{n,\omega}^{\mathrm{tr}} + \tilde{\lambda}_{n,\omega}^{\mathrm{tr}} \right)^{\top} q_{n,\omega}^{\mathrm{tr}^{*}} \\ + \sum_{n \in \mathcal{N}} \sum_{\omega \in \Omega} P_{\omega} \left( \bar{\lambda}_{n,\omega}^{\mathrm{br}} + \tilde{\lambda}_{n,\omega}^{\mathrm{br}} \right)^{\top} q_{n,\omega}^{\mathrm{br}^{*}} \\ + J(e^{*}, \delta^{\mathrm{E}^{*}}) \\ = J(e^{*}, \delta^{\mathrm{E}^{*}}) = J^{*}. \quad (38)$$

The second last equality follows from (30)-(32), which are inherently satisfied by the ADMM, as shown in [46]. ☐

The proposed algorithm is referred to as the **Consensus ADMM with Bounds (CADMMB)** and is summarized in Fig. 5. A tolerance parameter $\epsilon^{\mathrm{UB}}$ is used to guide whether the upper bound update step should be performed at each iteration. The termination criterion is given in (39):

$$J^{\mathrm{gap}}^{k} = 100 \frac{|J^{\mathrm{UB}}^{k} - J^{\mathrm{LB}}^{k}|}{|J^{\mathrm{UB}}^{k}|} \leq \epsilon^{\mathrm{gap}} \quad \text{or} \quad k \geq \bar{K}, \quad (39)$$

where $\epsilon^{\mathrm{gap}}$ and $\bar{K}$ are parameters of the algorithm.

A key advantage of this approach, in contrast to heuristic applications of distributed methods, lies in the availability of *trustworthy* bounds on the global optimum: the solution provided by the CADMMB is $\epsilon^{\mathrm{gap}}$-suboptimal, meaning its objective value is at most $\epsilon^{\mathrm{gap}}$ (e.g., 0.01%) worse than the value achieved by the optimal solution. In other words, the CADMMB offers a *performance guarantee* by providing solutions that operators can confidently evaluate and trust.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Illustration of the proposed CADMMB algorithm.</figcaption></figure>

<table><caption>TABLE II
CHARACTERIZATION OF THE UNCERTAINTIES.</caption><thead><tr><th rowspan="2">Uncertainty</th><th colspan="3">Average hourly values (standard deviation)</th></tr><tr><th>February</th><th>March</th><th>April</th></tr></thead><tbody><tr><td>π<sup>E</sup> (€/MWh)</td><td>51.1 (12.4)</td><td>35.4 (9.8)</td><td>34.7 (8.8)</td></tr><tr><td>π<sup>E↑</sup> (€/MWh)</td><td>58.1 (13.9)</td><td>38.6 (10.7)</td><td>38.6 (9.7)</td></tr><tr><td>π<sup>E↓</sup> (€/MWh)</td><td>49.1 (12.8)</td><td>32.3 (8.4)</td><td>28.5 (8.3)</td></tr><tr><td>Q<sup>ext</sup> (m<sup>3</sup>/s)</td><td>1367.8 (456.1)</td><td>1932.2 (549.3)</td><td>798.68 (219.0)</td></tr><tr><td>P<sup>S</sup> (MWh)</td><td>6.2 (10.7)</td><td>9.6 (14.1)</td><td>13.5 (17.4)</td></tr><tr><td>P<sup>W</sup> (MWh)</td><td>91.5 (61.9)</td><td>84.6 (57.1)</td><td>52.8 (29.5)</td></tr></tbody></table>

<table><caption>TABLE III
CHARACTERIZATION OF THE HYDROPOWER CASCADE.</caption><thead><tr><th>Asset name</th><th>Average annual energy output (GWh)</th><th>Head (m)</th><th>Capacity (MW)</th></tr></thead><tbody><tr><td>Pierre-Bénite</td><td>528</td><td>9.0</td><td>80</td></tr><tr><td>Vaugris</td><td>332</td><td>6.7</td><td>72</td></tr><tr><td>Péage de Roussillon</td><td>885</td><td>12.2</td><td>160</td></tr><tr><td>Saint Vallier</td><td>668</td><td>11.5</td><td>120</td></tr><tr><td>Bourg lès Valence</td><td>1082</td><td>11.7</td><td>180</td></tr><tr><td>Beauchastel</td><td>1211</td><td>11.8</td><td>192</td></tr><tr><td>Baix le Logis Neuf</td><td>1177</td><td>11.7</td><td>192</td></tr><tr><td>Montélimar</td><td>1575</td><td>16.5</td><td>275</td></tr><tr><td>Donzère-Mondragon</td><td>2032</td><td>22.5</td><td>349</td></tr><tr><td>Caderousse</td><td>843</td><td>8.6</td><td>187</td></tr><tr><td>Avignon</td><td>857</td><td>9.5</td><td>186</td></tr><tr><td>Vallabrègues</td><td>1269</td><td>11.3</td><td>210</td></tr></tbody></table>

### III. SIMULATION RESULTS AND DISCUSSION

This section presents simulation results over three months (February-April 2017), covering diverse inflow conditions. The uncertainties are characterized in Table II. Compagnie Nationale du Rhône's analysis shows that the uncertainties in vRES, inflows, and prices are minimally correlated and can be treated as independent [47]. Based on this finding, distinct sets of scenarios are generated for each uncertainty source, following the methodology utilized in the real-world operations of Compagnie Nationale du Rhône [47]. A scenario
<figure><img src="image_8.png" /><figcaption>Fig. 6. Out-of-sample validation results considering 2000 testing scenarios.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Example of convergence of the standard consensus ADMM algorithm for various initial penalty parameter values, when considering 5 scenarios.</figcaption></figure>

reduction process is then applied, wherein a random subset of the original set of scenarios is selected, with equal probability assigned to each scenario. To ensure accurate uncertainty representation, the selected scenarios must adequately capture the dispersion of the original set. Consequently, an out-of-sample validation procedure is employed to determine the minimum number of scenarios necessary to adequately represent the underlying uncertainty [48]. It is important to note that various methods exist for both scenario generation and reduction. A comprehensive discussion of these methods is outside the scope of this study but can be found in the existing literature dedicated to this topic [49]. However, it is worth emphasizing that no scenario generation method can generally guarantee a perfect a priori characterization of the uncertainty. Therefore, the significance of the analysis presented here lies in its out-of-sample validation results, which assess the reliability of the selected scenarios in accurately representing the underlying uncertainty, irrespective of the methods employed for scenario generation and reduction.

The hydropower cascade configuration is shown in Fig. 1 and detailed in Table III. Simulations use historical day-ahead market prices in France, with the bidding problem solved iteratively each day within a 4-hour optimization window. A piecewise linearization of the operational curves in Fig. 2 is employed using 40 equally spaced segments, i.e., $|S_n^{oc}| = 40$ for all $n \in \mathcal{N}$. In the CADMMB, both $\epsilon^{gap}$ and $\epsilon^{UB}$ are set to 0.01%, while $\bar{K} = 5000$. The algorithm runs on an Intel i7 processor using Gurobi 9.5.0 as the solver. The proposed CADMMB algorithm, illustrated in Fig. 5, is implemented using the Python's multiprocessing module [50], enabling fully parallel execution of the distributed subproblems. Each subproblem is assigned to an independent process with separate memory allocation. This design ensures true parallelism and reflects realistic operational conditions, as the derived subproblems correspond to distinct units in the hydropower cascade and are expected to run on separate computational devices with isolated memory.

## A. Out-of-Sample Validation

This subsection reports the out-of-sample validation results for the proposed stochastic bidding strategy in (19). This standard procedure assesses the performance of the stochastic optimization tool [48]. It involves solving problem (19) using a

<figure><img src="image_10.png" /><figcaption>Fig. 8. Example of convergence of the proposed CADMMB algorithm for various initial penalty parameter values, when considering 5 scenarios.</figcaption></figure>

set of training scenarios, then evaluating the objective function (16) by applying the optimal decisions from the training to a new testing set, consisting of 2000 scenarios in this case. The difference between the training and testing objective function values is called the *objective function error*.

Fig. 6 presents the observed out-of-sample validation results. For each selected number of training scenarios (5, 10, 15, 20, and 25, as indicated in the figure), 100 independent rounds of random sampling are conducted on the original set of generated scenarios. In each round, a distinct set of random samples is selected to train the proposed stochastic optimization tool, which then generates a corresponding training objective function value. The average daily difference between the training and testing objective function values is computed across all 100 rounds. This value is reported in the figure to assess the average out-of-sample performance of the proposed stochastic optimization tool for each selected number of training scenarios, accounting for the variability in outcomes arising from the random selection of training scenarios. The rationale for conducting multiple independent rounds of random scenario sampling for each specified number of training scenarios is to improve the robustness of the out-of-sample validation process. Fig. 6 reports both the observed mean and standard deviation (std) values.

Two distinct sets of out-of-sample validation results are presented in the figure. The first is obtained using the original formulation of the bidding strategy in (19), modeled as a MILP. The second is derived from an LP relaxation of (19),
<figure><img src="image_11.png" /><figcaption>Fig. 9. Example of ADMM performance when considering 20 scenarios.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Example of CADMMB performance when considering 20 scenarios.</figcaption></figure>

in which the integer variables are relaxed to take continuous values within the range [0, 1]. The results indicate that 20 scenarios are sufficient to minimize the error for both the MILP problem (19) and its LP relaxation. Using more than 20 scenarios raises computational complexity without significantly improving accuracy.

## B. Stylized Illustrative Example

Before presenting the full case study, results from a stylized example with 5 scenarios are presented. Figs. 7 and 8 show the convergence rates of the consensus ADMM and the proposed CADMMB, respectively, under different initial values of the penalty parameter $\rho$. The results reveal that the standard ADMM converges in less than 1140 iterations at best but exhibits undesirable oscillatory behavior. In contrast, Fig. 8 shows that the CADMMB not only eliminates the oscillations but also achieves bounds convergence with 99.99% accuracy in 203 iterations at best. The initial lower bound guess is already close to the global optimum and remains stable throughout the iterations. The initial gap between bounds is relatively small (2.54%), highlighting the benefits of warm starting the algorithm with initial bounds, as shown in Fig. 5. While consensus ADMM shows a satisfactory performance, the actual accuracy of its solutions remains indeterminate without explicitly computing the global optimum through centralized optimization. In contrast, our trustworthy CADMMB not only provides bounds that quantify the “distance” from the global optimum but also enables precise control over the desired accuracy by tuning the $\epsilon$ gap parameter.

## C. Full Scale Case Study

Figs. 9 and 10 illustrate the performance of the standard consensus ADMM and the proposed CADMMB, respectively,

<figure><img src="image_13.png" /><figcaption>Fig. 11. Comparison between centralized MILP, consensus ADMM, and CADMMB based on average runtime per scheduling period and failure rates (% of total simulation days), as a function of the number of scenarios.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Example of price-quantity bidding curves derived via CADMMB. Peak price periods are highlighted in red.</figcaption></figure>

<table><caption>TABLE IV
COMPARISON BETWEEN EX-POST PROFIT WITH FIXED HOURLY BIDS AND PRICE-QUANTITY BIDDING CURVES. IN BRACKETS, THE RELATIVE DIFFERENCE WITH RESPECT TO FIXED HOURLY BIDS.</caption><thead><tr><th rowspan="2">Num. scenarios</th><th colspan="2">Profit per MWh of energy produced (€)</th></tr><tr><th>Fixed bids</th><th>Price-quantity bidding curves</th></tr></thead><tbody><tr><td>5</td><td>33.63</td><td>34.38 (+2.2%)</td></tr><tr><td>10</td><td>35.74</td><td>36.80 (+3.0%)</td></tr><tr><td>15</td><td>37.82</td><td>39.22 (+3.7%)</td></tr><tr><td>20</td><td>39.45</td><td>41.22 (+4.5%)</td></tr><tr><td>25</td><td>39.89</td><td>41.36 (+3.7%)</td></tr></tbody></table>

when applied to the full scale case study involving 20 scenarios. As shown in Fig. 9, the ADMM fails to address the complexity of the problem. In contrast, our CADMMB exhibits a robust performance, efficiently converging with 99.99% accuracy in about 550 iterations. The initial lower bound estimate in this case is significantly distant from the optimum, necessitating iterative updates to refine $J^{LB}$. Notably, the evolution of $J^{LB}$ and $J^{UB}$ throughout the iterations differs markedly. While $J^{LB}$ is continuously updated based on the progressive evolution of the dual variables, $J^{UB}$ is determined through a projection step into a new feasible space, which results in the large jump observed in Fig. 10.

Fig. 11 generalizes these findings by showing the average
<table><caption>TABLE V
COMPARISON OF EX-POST PROFITS OBSERVED ACROSS THREE DISTINCT INFLOW CONDITIONS: ENERGETIC, LOW-FLOW, AND FLOOD.</caption><thead><tr><th rowspan="2">Num. scenarios</th><th colspan="3">Profit per MWh of energy produced (€)</th></tr><tr><th>Energetic</th><th>Low-flow</th><th>Flood</th></tr></thead><tbody><tr><td>5</td><td>38.81</td><td>33.56</td><td>31.16</td></tr><tr><td>10</td><td>41.58</td><td>35.99</td><td>33.13</td></tr><tr><td>15</td><td>44.88</td><td>38.56</td><td>34.79</td></tr><tr><td>20</td><td>48.03</td><td>40.05</td><td>35.31</td></tr><tr><td>25</td><td>48.63</td><td>40.39</td><td>35.32</td></tr></tbody></table>

runtime for solving the daily bidding problems and the number of failures—defined as the number of days for which the algorithm fails to provide a solution within the 4-hour optimization window—for different algorithms, as a function of the number of scenarios. In the reference case (20 scenarios), CADMMB reduces the average runtime by approximately 35% compared to consensus ADMM and by 57% compared to a traditional centralized optimization. Moreover, while the number of intractable instances explodes for both the standard ADMM and centralized algorithms as the number of scenarios grows, our trustworthy CADMMB consistently provides a distributed solution with over 99.99% accuracy ($\epsilon_{gap} = 0.01\%$).

Fig. 12 illustrates an example of the price-quantity bidding curves generated by our CADMMB. Despite the provided parallelization of the problem shown in Fig. 4, the CADMMB algorithm effectively coordinates both vRES and hydropower resources to optimize energy offers during daily peak price periods. These findings are generalized in Table IV, where it is reported that employing our CADMMB with price-quantity bidding curves enhances the average ex-post market profit of the VPP by up to 4.5% compared to using a traditional centralized solution with fixed hourly bids.

Finally, Fig. 13 illustrates a typical daily operation of the last hydropower asset in the cascade using CADMMB, reflecting the overall behavior of the hydropower system. The figure presents three distinct inflow conditions: energetic, low-flow, and flood periods. Due to constraints imposed by hydropower operational curves (e.g., Fig. 2), the hydropower plant can fully respond to vRES energy drops by discharging the reservoir only during energetic periods when sufficient storage capacity is available. During low-flow periods, CADMMB still enables to address the most significant vRES energy drops, while during flood periods, the hydropower plants are constrained to operate as run-of-the-river assets with no storage capacity.

These findings are generalized in Table V, which reports the average ex-post market profit achieved by the proposed CADMMB under varying inflow conditions. The simulation hours are grouped into three inflow regimes: energetic, low-flow, and flood periods, with corresponding average profits reported. The highest profits occur during energetic periods, when ample storage enables full utilization of the hydropower flexibility. In contrast, the limited flexibility during low-flow and flood periods significantly reduces the VPP profitability.

## IV. CONCLUSION AND FUTURE WORK

This study introduces a stochastic ADMM-based distributed algorithm for the co-optimization of vRES and CHPP in

<figure><img src="image_15.png" /><figcaption>Fig. 13. Example of coordination between the hydropower cascade and vRES under various inflow conditions when the proposed CADMMB is employed.</figcaption></figure>

the day-ahead market. First, the original MILP problem is decomposed across spatial and scenario dimensions. A distributed strategy is then applied to derive trustworthy bounds on the global optimum, with mathematical proof provided to demonstrate their existence and validity. Unlike a conventional ADMM application, which can lead to heuristic methods for non-convex problems, the proposed CADMMB algorithm offers a performance guarantee. Testing in collaboration with the French aggregator Compagnie Nationale du Rhône shows the algorithm's effectiveness in managing complex aggregations. The CADMMB achieves a 35% reduction in average runtime compared to consensus ADMM and a 57% reduction compared to traditional centralized optimization. Additionally, CADMMB consistently delivers solutions with over 99.99% accuracy, whereas both centralized and consensus ADMM methods fail to solve the problem. Future research will aim to refine initial bound guess strategies, explore the effects of parameter tuning, extend the analysis to multi-market trading case, and integrate asynchronous communication to enhance the robustness of the proposed distributed algorithm.

## REFERENCES

[1] D. Koraki, and K. Strunz, "Wind and solar power integration in electricity markets and distribution networks through service-centric virtual power plants," *IEEE Trans. Power Syst.*, vol. 33, no. 1, pp. 473-485, Jan. 2018.

[2] International Energy Agency, "Hydropower special market report," 2021. [Online]. Available: https://www.iea.org/reports/hydropower-special-market-report. [Accessed: Oct. 9, 2024].

[3] J. Zhang, C. Cheng, S. Yu, and H. Su, "Chance-constrained co-optimization for day-ahead generation and reserve scheduling of cascade hydropower-variable renewable energy hybrid systems," *Appl. Energy*, vol. 324, pp. 119732, Jul. 2022.

[4] Z. Zhao, C. Cheng, S. Liao, Y. Li, and Q. Lü, "A MILP based framework for the hydro unit commitment considering irregular forbidden zone related constraints," *IEEE Trans. Power Syst.*, vol. 36, no. 3, pp. 1819-1832, May 2021.

[5] D. Apostolopoulou, Z. De Grève, and M. McCulloch, "Robust optimization for hydroelectric system operation under uncertainty," *IEEE Trans. Power Syst.*, vol. 33, no. 3, pp. 3337-3348, May 2018.

[6] B. Zhou, G. Geng, and Q. Jiang, "Hydro-thermal-wind coordination in day-ahead unit commitment," *IEEE Trans. Power Syst.*, vol. 31, no. 6, pp. 4626-4637, Nov. 2016.
[7] Y. Qiu, J. Lin, F. Liu, Y. Song, G. Chen, and L. Ding, "Stochastic online generation control of cascaded run-of-the-river hydropower for mitigating solar power volatility," *IEEE Trans. Power Syst.*, vol. 35, no. 6, pp. 4709-4722, Nov. 2020.

[8] S. E. Fleten, and T. K. Kristoffersen, "Stochastic programming for optimizing bidding strategies of a Nordic hydropower producer," *Eur. J. Oper. Res.*, vol. 181, no. 2, pp. 916-928, Sep. 2007.

[9] L. A. Roald, D. Pozo, A. Papavasiliou, D. K. Molzahn, J. Kazempour, and A. J. Conejo, "Power systems optimization under uncertainty: A review of methods and applications," *Elect. Power Syst. Res.*, vol. 214, pp. 108725, Jan. 2023.

[10] Z. Li, P. Yang, Y. Yang, G. Lu, and Y. Tang, "Hydropower preventive maintenance scheduling in a deregulated market: A multi-stage stochastic programming approach," *Electr. Power Syst. Res.*, vol. 221, pp. 109379, Aug. 2023.

[11] Z. Zhong, N. Fan, and L. Wu, "Multistage stochastic optimization for mid-term integrated generation and maintenance scheduling of cascaded hydroelectric system with renewable energy uncertainty," *Eur. J. Oper. Res.*, vol. 318, no. 1, pp. 179-199, Oct. 2024.

[12] Y. Liu, X. Chen, N. Fan, Z. Zhao, and L. Wu, "Stochastic day-ahead operation of cascaded hydropower systems with Bayesian neural network-based scenario generation: A Portland general electric system study," *Int. J. Electr. Power Energy Syst.*, vol. 153, pp. 109327, Nov. 2023.

[13] A. Helseth, M. Fodstad, and B. Mo, "Optimal medium-term hydropower scheduling considering energy and reserve capacity markets," *IEEE Trans. Sustain. Energy*, vol. 7, no. 3, pp. 934-942, Jul. 2016.

[14] C. Wang et al., "Resilience-based coordinated scheduling of cascaded hydro power with sequential heavy precipitation," *IEEE Trans. Sustain. Energy*, vol. 14, no. 2, pp. 1299-1311, Apr. 2023.

[15] Y. Shi, H. Wang, C. Li, M. Negnevitsky, and X. Wang, "Stochastic optimization of system configurations and operation of hybrid cascade hydro-wind-photovoltaic with battery for uncertain medium- and long-term load growth," *Appl. Energy*, vol. 364, pp. 123127, Jun. 2024.

[16] H. Zhang, Z. Lu, W. Hu, Y. Wang, L. Dong, and J. Zhang, "Coordinated optimal operation of hydro-wind-solar integrated systems," *Appl. Energy*, vol. 242, pp. 883-896, May 2019.

[17] N. Mazzi, J. Kazempour, and P. Pinson, "Price-taker offering strategy in electricity pay-as-bid markets," *IEEE Trans. Power Syst.*, vol. 33, no. 2, pp. 2175-2183, Mar. 2018.

[18] M. Karasavvidis, D. Papadaskalopoulos, and G. Strbac, "Optimal offering of energy storage in electricity markets with loop blocks," *IEEE Trans. Power Syst.*, vol. 39, no. 1, pp. 299-315, Jan. 2024.

[19] Y. Yu, Y. Wu, and Q. Sheng, "Optimal scheduling strategy of cascade hydropower plants under the joint market of day-ahead energy and frequency regulation," *IEEE Access*, vol. 9, pp. 87749-87762, Jun. 2021.

[20] S. Séguin, P. Côté, and C. Audet, "Self-scheduling short-term unit commitment and loading problem," *IEEE Trans. Power Syst.*, vol. 31, no. 1, pp. 133-142, Jan. 2016.

[21] E. Shayesteh, M. Amelin, and L. Söder, "Multi-station equivalents for short-term hydropower scheduling," *IEEE Trans. Power Syst.*, vol. 31, no. 6, pp. 4616-4625, Nov. 2016.

[22] L. S. M. Guedes, P. de Mendonça Maia, A. C. Lisboa, D. A. G. Vieira, and R. R. Saldanha, "A unit commitment algorithm and a compact MILP model for short-term hydro-power generation scheduling," *IEEE Trans. Power Syst.*, vol. 32, no. 5, pp. 3381-3390, Sep. 2017.

[23] A. Helseth, and B. Mo, "Hydropower aggregation by spatial decomposition—An SDDP approach," *IEEE Trans. Sustain. Energy*, vol. 14, no. 1, pp. 381-392, Jan. 2023.

[24] C. J. López-Salgado, O. Añó, and D. M. Ojeda-Esteybar, "Stochastic unit commitment and optimal allocation of reserves: A hybrid decomposition approach," *IEEE Trans. Power Syst.*, vol. 33, no. 5, pp. 5542-5552, Sep. 2018.

[25] D. Apostolopoulou, and M. McCulloch, "Optimal short-term operation of a cascaded hydro-solar hybrid system: A case study in Kenya," *IEEE Trans. Sustain. Energy*, vol. 10, no. 4, pp. 1878-1889, Oct. 2019.

[26] N. Lu, G. Wang, C. Su, Z. Ren, X. Peng, and Q. Sui, "Medium- and long-term interval optimal scheduling of cascade hydropower-photovoltaic complementary systems considering multiple uncertainties," *Appl. Energy*, vol. 353, pp. 122085, Jan. 2024.

[27] E. C. Finardi, E. L. Silva, and C. Sagastizabal, "Solving the unit commitment problem of hydro plants via lagrangian relaxation and sequential quadratic programming," *Computational Appl. Math*, vol. 24, pp. 317-341, Dec. 2005.

[28] S. Takriti, J. R. Birge, and E. Long, "A stochastic model for the unit commitment problem," *IEEE Trans. Power Syst.*, vol. 11, no. 3, pp. 1497-1508, Aug. 1996.

[29] J. Huang, X. Wu, Z. Zheng, Y. Huang, and W. Li, "Multi-objective optimal operation of combined cascade reservoir and hydrogen system," *IEEE Trans. Ind Appl.*, vol. 58, no. 2, pp. 2836-2847, Apr. 2022.

[30] B. Tong, Q. Zhai, and X. Guan, "An MILP based formulation for short-term hydro generation scheduling with analysis of the linearization effects on solution feasibility," *IEEE Trans. Power Syst.*, vol. 28, no. 4, pp. 3588-3599, Nov. 2013.

[31] R. Egging, S.-E. Fleten, I. Grønvik, A. Hadziomerovic, and N. Ingvoldstad, "Linear decision rules for hydropower scheduling under uncertainty," *IEEE Trans. Power Syst.*, vol. 32, no. 1, pp. 103-113, Jan. 2017.

[32] M. R. Scuzziato, E. C. Finardi, and A. Frangioni, "Comparing spatial and scenario decomposition for stochastic hydrothermal unit commitment problems," *IEEE Trans. Sustain. Energy*, vol. 9, no. 3, pp. 1307-1317, Jul. 2018.

[33] A. Helseth, S. Jaehnert, and A. L. Diniz, "Convex relaxations of the short-term hydrothermal scheduling problem," *IEEE Trans. Power Syst.*, vol. 36, no. 4, pp. 3293-3304, Jul. 2021.

[34] Y. Ma et al., "Decentralized monthly generation scheduling of cascade hydropower plants in multiple time scale markets," *Int. J. Elect. Power Energy Syst.*, vol. 135, pp. 107420, Feb. 2022.

[35] L. Santosuosso, S. Camal, A. Lett, G. Bontron, and G. Kariniotakis, "Distributed economic model predictive control for the joint energy dispatch of wind farms and run-of-the-river hydropower plants," *Elect. Power Syst. Res.*, vol. 235, pp. 110805, Oct. 2024.

[36] D. Gade, G. Hackebeil, S. M. Ryan, J. P. Watson, R. J. B. Wets, and D. L. Woodruff, "Obtaining lower bounds from the progressive hedging algorithm for stochastic mixed-integer programs," *Math. Program.*, vol. 157, pp. 47-67, Apr. 2016.

[37] V. Piron, G. Bontron, and M. Pochat, "Operating a hydropower cascade to optimize energy management," *Hydropower & Dams*, vol. 22, no. 5, pp. 2-6, Sep. 2015.

[38] A. Baringo, and L. Baringo, "A stochastic adaptive robust optimization approach for the offering strategy of a virtual power plant," *IEEE Trans. Power Syst.*, vol. 32, no. 5, pp. 3492-3504, Sept. 2017.

[39] N. Mazzi, A. Trivella, and J. M. Morales, "Enabling active/passive electricity trading in dual-price balancing markets," *IEEE Trans. Power Syst.*, vol. 34, no. 3, pp. 1980-1990, May 2019.

[40] G. Steeger, L. A. Barroso, and S. Rebennack, "Optimal bidding strategies for hydro-electric producers: A literature survey," *IEEE Trans. Power Syst.*, vol. 29, no. 4, pp. 1758-1766, Jul. 2014.

[41] J. S. Anagnostopoulos, and D. E. Papantonis, "Optimal sizing of a run-of-river small hydropower plant," *Energy Convers. Manag.*, vol. 48, no. 10, pp. 2663-2670, Oct. 2007.

[42] G. P. McCormick, "Computability of global solutions to factorable nonconvex programs: Part I—Convex underestimating problems," *Math. Program.*, vol. 10, no. 1, pp. 147-175, Dec. 1976.

[43] J. Najman, and A. Mitsos, "On tightness and anchoring of McCormick and other relaxations," *J. Global Optim.*, vol. 74, pp. 677-703, Dec. 2017.

[44] E. Blom, and L. Söder, "Single-level reduction of the hydropower area equivalent bilevel problem for fast computation," *Renew. Energy*, vol. 225, pp. 120229, May 2024.

[45] B. Flamm, A. Eichler, J. Warrington, and J. Lygeros, "Two-stage dual dynamic programming with application to nonlinear hydro scheduling," *IEEE Trans. Control Syst. Technol.*, vol. 29, no. 1, pp. 96-107, Jan. 2021.

[46] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, "Distributed optimization and statistical learning via the alternating direction method of multipliers," *Found. Trends Mach. Learn.*, vol. 3, no. 1, pp. 1-122, Nov. 2011.

[47] L. Santosuosso, "Distributed stochastic optimization for operating complex virtual power plants: Leveraging cascaded run-of-the-river hydropower flexibility for renewable energy integration," Ph.D. dissertation, Université Paris Sciences et Lettres, Paris, France, 2025.

[48] A. Arrigo, C. Ordoudis, J. Kazempour, Z. De Grève, J.F. Toubeau, and F. Vallée, "Wasserstein distributionally robust chance-constrained optimization for energy and reserve dispatch: An exact and physically-bounded formulation," *Eur. J. Oper. Res.*, vol. 296, no. 1, pp. 304-322, Jan. 2022.

[49] H. Li et al., "A review of scenario analysis methods in planning and operation of modern power systems: Methodologies, applications, and challenges," *Electr. Power Syst. Res.*, vol. 205, pp. 107722, Apr. 2022.

[50] Python Software Foundation, "Process-based parallelism," 2025. [Online]. Available: https://docs.python.org/3/library/multiprocessing.html. [Accessed: May 26, 2025].